matlab求解ode系统

Question

我有 9 个具有时间相关系数的方程 g

% MY M file
function dy =tarak(t,y)
G= 3.16;
g =  0.1*exp(-((t-200)/90).^2);
dy=zeros(9,1);
dy(1)=-2*2*y(1)+2*G*y(5)+2*g*y(7);
dy(2)=2*y(1)-2*G*y(5);
dy(3)=2*y(1)-2*g*y(7);
dy(4)=-2*y(4)+g*y(9);
dy(5)=-2*y(5)+G*(y(2)-y(1))+g*y(8);
dy(6)=-2*y(6)-G*y(9);
dy(7)=-2*y(7)+g*(y(3)-y(1))+G*y(8);
dy(8)=-G*y(7)-g*y(5);
dy(9)=G*y(6)-g*y(4);

然后在命令窗口中：

[T,Y] = ode45(@tarak,[0 ,500],[0 0 1 0 0 0 0 0 0])

其中系数G = 3.16和g = 0.1*exp(-((t-200)/90).^2)是一个时间相关系数和时间t = 0:500；初始条件[0 0 1 0 0 0 0 0 0]。

我的输出 y(1)、y(2) 得到了错误的负值。有人可以尝试用ode45 解决上述方程，以便我可以比较结果。

Answer 1

通过 RK4 的简单应用，我得到了这张图片

非常积极，y(1) 组件有一个奇怪的初始跳跃。但请注意规模，总的来说y(1) 相当小。此时系统似乎很僵硬，所以 rk45 可能会出现问题，隐式 Runge-Kutta 方法会更好。

以及初始振荡的缩放

---

Python 代码

import numpy as np
import matplotlib.pyplot as plt
from math import exp

def dydt(t,y):
    dy = np.array(y);

    G = 3.16;
    g = 0.1*exp(-((t-200)/90)**2);

    dy[0]=-2*2*y[0]+2*G*y[4]+2*g*y[6];
    dy[1]=   2*y[0]-2*G*y[4];
    dy[2]=   2*y[0]-2*g*y[6];
    dy[3]=  -2*y[3]+  g*y[8];
    dy[4]=  -2*y[4]+  G*(y[1]-y[0])+g*y[7];
    dy[5]=  -2*y[5]-  G*y[8];
    dy[6]=  -2*y[6]+  g*(y[2]-y[0])+G*y[7];
    dy[7]=  -G*y[6]-  g*y[4];
    dy[8]=   G*y[5]-  g*y[3];
    return dy;

def RK4Step(f,x,y,h):
    k1=f(x      , y         )
    k2=f(x+0.5*h, y+0.5*h*k1)
    k3=f(x+0.5*h, y+0.5*h*k2)
    k4=f(x+    h, y+    h*k3)
    return (k1+2*(k2+k3)+k4)/6.0


t= np.linspace(0,500,200+1);
dt = t[1]-t[0];
y0=np.array([0, 0, 1, 0, 0, 0, 0, 0, 0]);

y = [y0]

for t0 in t[0:-1]:
    N=200;
    h = dt/N;
    for i in range(N):
        y0 = y0 + h*RK4Step(dydt,t0+i*h,y0,h);
    y.append(y0);

y = np.array(y);

plt.subplot(121);
plt.title("y(1)")
plt.plot(t,y[:,0],"b.--")
plt.subplot(122);
plt.title("y(2)")
plt.plot(t,y[:,1],"b-..")
plt.show()

Answer 2

在Matlab:

options = odeset('AbsTol', 1e-12);
[T,Y] = ode45(@tarak, [0, 500], [0 0 1 0 0 0 0 0 0], options);