計算ODE與MATLAB

Question

微分方程：

α '（T）= S（β-βα+α-Qα^ 2）

β'（T）=（S^-1 ）（ - β-αβ+γ）

γ'（T）= W（α-γ）

初使值

α（0）= 30.00

β（0）= 1.000

γ（0）= 30.00

計算

我要解決從T_0 = 0的問題至t = 10，而使用數值• = 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)') 

導致這個數字：