Sunday, 13 March 2022

Step by step time integrators in Python

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 ?



from Step by step time integrators in Python

No comments:

Post a Comment