Step by step time integrators in Python

Question

I am solving a first order initial value problem of the form:

dy/dt = f(t,y(t)), y(0)=y0

I would like to obtain y(n+1) from a given numerical scheme, like for example :

using explicit Euler's scheme, we have

y(i) = y(i-1) + f(t-1,y(t-1)) * dt

Example code:

# Test code to evaluate different time integrators for the following equation:
# y' = (1/2) y + 2sin(3t) ; y(0) = -24/37

def dy_dt(y,t):
    
    func = (1/2)*y + 2*np.sin(3*t)
    
    return func


import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint


tmin = 0
tmax = 50
delt= 1e-2

t = np.arange(tmin,tmax,delt)

total_steps = len(t)

y_explicit=np.zeros(total_steps)
#y_ODEint=np.zeros(total_steps)


y0 = -24/37

y_explicit[0]=y0
#y_ODEint[0]=y0


 # exact solution
 
y_exact = -(24/37)*np.cos(3*t)- (4/37)*np.sin(3*t) + (y0+24/37)*np.exp(0.5*t)


# Solution using ODEint Python

y_ODEint = odeint(dy_dt,y0,t)



for i in range(1,total_steps):
    
 # Explicit scheme  
 
 y_explicit[i] = y_explicit[i-1] + (dy_dt(y_explicit[i-1],t[i-1]))*delt   

 # Update using ODEint 
 
 # y_ODEint[i] = odeint(dy_dt,y_ODEint[i-1],[0,delt])[-1]
  


plt.figure()

plt.plot(t,y_exact)
plt.plot(t,y_explicit)
# plt.plot(t,y_ODEint)

The current issue I am having is that the functions like ODEint in python provide the entire y(t) as opposed to y(i). like in the line "y_ODEint = odeint(dy_dt,y0,t)"

See in the code, how I have coded the explicit scheme, which gives y(i) for every time step. I want to do the same with ODEint, i tried something but didn't work (all commented lines)

I want to obtain y(i) rather than all ys using ODEint. Is that possible ?

Answer 1

Your system is time variant so you cannot translate the time step from (t[i-1], t[i]) to (0, delt).

The step by step integration will is unstable for your differential equation though

Here is what I get

def dy_dt(y,t):
    
    func = (1/2)*y + 2*np.sin(3*t)
    
    return func

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint


tmin = 0
tmax = 40
delt= 1e-2

t = np.arange(tmin,tmax,delt)

total_steps = len(t)

y_explicit=np.zeros(total_steps)
#y_ODEint=np.zeros(total_steps)


y0 = -24/37

y_explicit[0]=y0

# exact solution

y_exact = -(24/37)*np.cos(3*t)- (4/37)*np.sin(3*t) + (y0+24/37)*np.exp(0.5*t)


# Solution using ODEint Python

y_ODEint = odeint(dy_dt,y0,t)
# To be filled step by step
y_ODEint_2 = np.zeros_like(y_ODEint)
y_ODEint_2[0] = y0

for i in range(1,len(y_ODEint_2)):
    # update your code to run with the correct time interval
    y_ODEint_2[i] = odeint(dy_dt,y_ODEint_2[i-1],[tmin+(i-1)*delt,tmin+i*delt])[-1]

plt.figure()

plt.plot(t,y_ODEint, label='single run')
plt.plot(t,y_ODEint_2, label='step-by-step')
plt.plot(t, y_exact, label='exact')
plt.legend()
plt.ylim([-20, 20])
plt.grid()

Important to notice that both methods are unstable, but the step-by-step explodes slightly before than the single odeint call.

With, for example dy_dt(y,t): -(1/2)*y + 2*np.sin(3*t) the integration becomes more stable, for instance, there is no noticeable error after integrating from zero to 200.