Search code examples
pythonscipynumerical-integrationodeint

Scipy integrate: is solve_ivp "always worse" than odeint?

Analytical solution for an example ODE

I was testing the difference between these methods by solving the following initial value problem:

y'=2*y-t

You can solve this analytically by considering that y(t) is a linear combination of a homogeneous solution y_h(t)=c1*exp(2t) and a particular solution y_p(t)=t/2+1/4. The constant c1 is found by substituting y(t0)=y0. The analytical solution is then:

y(t)=(y0-t0/2-1/4)*exp(2*(t-t0))+t/2+1/4

Note that if y0=0.25 and t0=0, this is the same as y(t)=t/2+1/4. In this case y(1)=0.75.

Comparison between solve_ivp and odeint

First from scipy.integrate import solve_ivp, odeint.

By writing odeint(lambda y,t: 2*y-t,[0.25],[0,1]) we get the expected results y(0)=0.25 and y(1)=0.75.

But, by writing solve_ivp(lambda y,t: 2*y-t,t_span=[0,1],y0=[0.25],t_eval=[0,1]) we get the results y(0)=0.25 and y(1)=0.82775742.

As mentioned in this question, solve_ivp should have the 'LSODA' method and have its tolerances adjusted in order to fairly compare it to odeint. By reading the scipy odeint documentation we see that the tolerances are around 1.49e-8.

But solve_ivp(lambda y,t: 2*y-t,t_span=[0,1],y0=[0.25],t_eval=[0,1],method='LSODA',atol=1.49e-8,rtol=1.49e-8) still yields y(0)=0.25 and y(1)=0.82772876.

And if you try this for greater time spans, the results only get worse with solve_ivp, for this particular example.

Am I missing something?

Solution

  • I need to sleep...

    According to the odeint documentation, the function has to be of the type func(y, t) (although it does accept f(t,y) if you specify it).

    And according to the solve_ivp documentation, the function has to be of the type func(t, y).