matlab - 使用 MATLAB 计算 ODE

Question

微分方程：

α'(t)=s(β-βα+α-qα^2)

β'(t)=(s^-1)(-β-αβ+γ)

γ'(t)=w(α-γ)

初始值

α (0) = 30.00

β (0) = 1.000

γ (0) = 30.00

计算

我想解决从 t_0=0 到 t=10 的问题，同时使用值 s = 1、q = 1 和 w = 0.1610

我不知道如何为 ODE 编写函数，非常感谢您的帮助！

Answer 1

我通常不习惯解决别人的作业，但我猜今天是你的幸运日。

所以，你有一个耦合常微分方程系统：

α'(t) = s(β-α(β+1)-qα²)

β'(t) = (-β-αβ+γ)/s

γ'(t) = w(α-γ)

你想解决

y = [α(t) β(t) γ(t)]

0 < t < 10，s = 1，q = 1，w = 0.1610。在 matlab 中执行此操作的方法是定义一个计算导数 ([α'(t) β'(t) γ'(t)]) 的函数，并将其放入其中一个 ODE 求解器ode45中赌注）：

s = 1;
q = 1; 
w = 0.1610;

% Define y(t) = [α(t) β(t) γ(t)] = [y(1) y(2) y(3)]:

deriv = @(t,y) [...
    s * (y(2) - y(1)*(y(2)+1) - q*y(1)^2)   % α'(t)
    (-y(2) - y(1)*y(2) + y(3))/s            % β'(t)
    w * (y(1)-y(3))                         % γ'(t)
];

% initial value
y0 = [30 1 10];

% time span to integrate over
tspan = [0 10];

% solve ODE numerically
[t, y] = ode45(deriv, tspan, y0)

这将输出

y =
   30.0000    1.0000   10.0000
   28.5635    0.9689   10.0049   % numerical solutions at times t
   27.2558    0.9413   10.0094
   26.0603    0.9166   10.0136
   ...        ...      ...
   = α(t)     = β(t)   = γ(t)

t =
         0
    0.0016
    0.0031   % corresponding times
    0.0047
    ...

您可以像这样绘制所有这些：

figure, clf, hold on
plot(t, y(:,1), 'r')
plot(t, y(:,2), 'g')
plot(t, y(:,3), 'b')
legend('\alpha(t)', '\beta(t)', '\gamma(t)')

这导致了这个数字：