3

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?

Breno
  • 748
  • 8
  • 18
  • I am stuck on the same question for a different system. `solve_ivp` also gives me worse results than `odeint` even though I'm using the kwargs `method='LSODA',atol=1.49e-8,rtol=1.49e-8` in `solve_ivp` and even though I'm using the correct order of y and t. Did you come across any other insights along the way? – jollycat Jan 15 '21 at 09:16

1 Answers1

0

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

Breno
  • 748
  • 8
  • 18
  • In recent numpy version, the trick is to give `tfirst=True` to `odeint`, this way the interface of the `odeint` and `solve_ivp` are the same. `dqdt = lambda t,y: 2*y-t`, `res1 = odeint (dqdt, y0, t_eval, tfirst=True)` and `res2 = solve_ivp(dqdt, (t_eval[0],t_eval[-1]), y0, method='LSODA', t_eval=t_eval, atol=1.49e-8, rtol=1.49e-8)`. This snippet also make sure the returned evaluation times are the same for both. – e-malito Jul 19 '22 at 21:55