Is it possible to update funcion variables every time it's been called?

288 Views Asked by At

I'm trying to solve an ODE system with solve_ivp and i want to change the local variables of the function every time it's been called by the solver. In particular I wand to update the lagrange multipliers (lambdas_i) so that the next time step of solve_ivp, uses the values of the previous one. (The ''reconstruct'' function is from a python module that uses a method to reconstruct a size distribution from given moments) Is there a way to do this? I'll post the code below:

import time 
import numpy as np
import scipy.integrate as integrate
from numpy import sqrt, sin, cos, pi
import math 
import pylab as pp
from pymaxent import reconstruct


start = time.time()


'''Initialize variables'''

t=np.linspace(0, 60,31)
tspan=(0,60)

initial_m=[]
for i in range(4):
    def distr(L,i=i):
        return (L**i)*0.0399*np.exp(-((L-50)**2)/200)
   
    m, err=integrate.quad(distr, 0, np.inf)
    print('m(',i,')=',m)
    initial_m.append(m)
    

''' Solving ode system using Maximum Entropy, G(L)=1+0.002*L'''

def moments(t,y):
    m0 = y[0]
    m1 = y[1]
    m2 = y[2]
    m3 = y[3]
    Lmean=m1
    σ=(m2-Lmean**2)**(1/2)
    Lmin=Lmean-3*σ
    Lmax=Lmean+4*σ
    bnds=[Lmin,Lmax]
    L=np.linspace(Lmin,Lmax)
    
    sol, lambdas_i= reconstruct(mu=y ,bnds=bnds)
    print(lambdas_i)
    
    dm0dt=0
        
    def moment1(L):
        return(sol(L)+0.002*sol(L)*L)                
    dm1dt, err1=integrate.quad(moment1,Lmin,Lmax)
    
    
    def moment2(L):
        return(2*L*(sol(L)+0.002*sol(L)*L))
    
    dm2dt, err2=integrate.quad(moment2,Lmin,Lmax)
    
    def moment3(L):
        return(3*L**2*(sol(L)+0.002*sol(L)*L))
    
    dm3dt, err3=integrate.quad(moment3,Lmin,Lmax)
    
    return(dm0dt,dm1dt,dm2dt,dm3dt)

'''Χρήση της BDF, step by step'''

r=integrate.solve_ivp(moments,tspan,initial_m,method='BDF',jac=None,t_eval=t,rtol=10**(-3))

end = time.time()

print('Total time =',{end-start})
1

There are 1 best solutions below

0
On

Here is one way to achieve what you want. I won't be using your actual code, but using a simpler example problem, and you can use the same strategy to solve yours.

def seq():
  l = [1, 0]
  while True:
    p = l[1]
    n = l[0] + p
    l = [p, n]
    print(n)
    input()

Above is an example function that would print the next (or first) term of the fibonacci sequence. (In which each term is the sum of the two previous terms.) In this case, the input is solely used to pause between each iteration (as it is infinite). Now to transform this into a generator, to allow more flexibility, you could rewrite the function as:

def seq():
  l = [1, 0]
  while True:
    p = l[1]
    n = l[0] + p
    l = [p, n]
    yield n

Now, if you wanted to get the same result, you could:

for item in seq():
  print(item)
  input()

However, this is mostly useless. The point of the generator comes in when you want to gather the next number of the sequence, but at any point of the code. Not necessarily inside the loop that repeats itself until you're done. You could achieve this:

gen = seq()
next(gen) # returns 1
next(gen) # returns 1
next(gen) # returns 2
next(gen) # returns 3
next(gen) # returns 5

And so on...


Another way to solve this problem would be to use a global variable instead of a local variable.

l = [1, 0]

def seq():
  p = l[1]
  n = l[0] + p
  l[0] = p
  l[1] = n
  return n

The variable l is defined outside the function, not inside, so it will not be discarded when the function exits. (To run the code the same way as before:)

while True:
  print(seq())
  input()

Both of these should be implementable in your code.