I need to understand the mechanism of scipy.integrate.LSODA function.
I have written a test script that integrate a simple function. According to the LSODA webpage inputs of functions can be rhs function, t_min, initial y and t_max. On the other hand, when I run the code, I get nothing. What should I do?
import scipy.integrate as integ
import numpy as np
def func(t,y):
return t
y0=np.array([1])
t_min=1
t_max=10
N_max=100
t_min2=np.linspace(t_min,t_max,N_max)
first_step=0.01
solution=integ.LSODA(func,t_min,y0,t_max)
solution2=integ.odeint(func,y0,t_min2)
print(solution.t,solution.y,solution.nfev,'\n')
print(solution2)
Solution give
1 [ 1.] 0
[[ 1.00000000e+00]
[ 9.48773662e+00]
[ 9.00171421e+01]
[ 8.54058901e+02]
[ 8.10308559e+03]]
1.) You only initiate an instance of the LSODA solver class, no computation occurs, just initialization of the arrays with the initial data. To get an odeint-like interface, use solve_ivp with the option method='LSODA'.
2.) Without the option tfirst=True, the LSODA solver will solve y'(t)=t, while odeint will solve y'(t)=y(t)
To get comparable results, one should also equalize the tolerances, as the default tolerances can be different. One thus can call the methods like
print "LSODA"
solution=integ.solve_ivp(func,[t_min, t_max],y0,method='LSODA', atol=1e-4, rtol=1e-6)
print "odeint"
solution2=integ.odeint(func,y0,t_min2, tfirst=True, atol=1e-4, rtol=1e-6)
Even then you get no information on the internal steps of odeint, even if the FORTRAN code has an option for that, the python wrapper does not replicate it. You could add a print statement to the ODE function func so that you see at what places this function is actually called, this should average to 2 calls with close-by arguments per internal step.
This reports
LSODA
1.0 [1.]
1.00995134265 [1.00995134]
1.00995134265 [1.01005037]
1.01990268529 [1.02010074]
1.01990268529 [1.02019977]
10.0 [50.50009903]
10.0 [50.50009903]
odeint
1.0 [1.]
1.00109084546 [1.00109085]
1.00109084546 [1.00109204]
1.00218169092 [1.00218407]
1.00218169092 [1.00218526]
11.9106363102 [71.43162985]
where the reported steps in the output of LSODA are
[ 1. 1.00995134 1.01990269 10. ] [[ 1. 1.01005037 1.02019977 50.50009903]] 7
Of course, a high-order method will integrate the linear polynomial y'=t to the quadratic polynomial y(t)=0.5*(t^2+1) with essentially no error independent of the step size.