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})