Python 中的逐步时间积分器

Question

我正在求解以下形式的一阶初始值问题：

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

我想从给定的数值方案中获得 y(n+1)，例如：

使用显式欧拉方案，我们有

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

示例代码：

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

我目前遇到的问题是 python 中的 ODEint 等函数提供了整个 y(t) 而不是 y(i)。就像在行“y_ODEint = odeint(dy_dt,y0,t)”

请参阅代码，我是如何编写显式方案的，它为每个时间步长提供 y(i)。我想对 ODEint 做同样的事情，我尝试了一些但没有奏效（所有注释行）

我想使用 ODEint 获取 y(i) 而不是所有 ys。这可能吗？

Answer 1

您的系统是时变的，因此您无法将时间步长从 (t[i-1], t[i]) 转换为 (0, delt)。

虽然你的微分方程的逐步积分是不稳定的

这是我得到的



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

重要的是要注意这两种方法都不稳定，但是 step-by-step 比单个 odeint 调用稍微早一点爆炸。

例如dy_dt(y,t): -(1/2)*y + 2*np.sin(3*t)积分变得更加稳定，例如从零积分到200后没有明显的错误。

Python 中的逐步时间积分器

Step by step time integrators in Python

python

scipy