龍格 - 庫塔用於耦合常微分方程

Question

我建立在八度的功能，可以解決的類型的N耦合常微分方程：

dx/dt = F(x,y,…,z,t) 
dy/dt = G(x,y,…,z,t) 
dz/dt = H(x,y,…,z,t) 

與任何這三種方法（歐拉，威享和龍格 - 庫塔的-4）。

下面的代碼對應於該函數：

function sol = coupled_ode(E, dfuns, steps, a, b, ini, method) 
    range = b-a; 
    h=range/steps; 
    rows = (range/h)+1; 
    columns = size(dfuns)(2)+1; 
    sol= zeros(abs(rows),columns); 
    heun=zeros(1,columns-1); 
    for i=1:abs(rows) 
    if i==1 
     sol(i,1)=a; 
    else 
     sol(i,1)=sol(i-1,1)+h;  
    end 
    for j=2:columns 
     if i==1 
     sol(i,j)=ini(j-1); 
     else 
     if strcmp("euler",method) 
      sol(i,j)=sol(i-1,j)+h*dfuns{j-1}(E, sol(i-1,1:end));  
     elseif strcmp("heun",method) 
      heun(j-1)=sol(i-1,j)+h*dfuns{j-1}(E, sol(i-1,1:end));   
     elseif strcmp("rk4",method) 
      k1=h*dfuns{j-1}(E, [sol(i-1,1), sol(i-1,2:end)]); 
      k2=h*dfuns{j-1}(E, [sol(i-1,1)+(0.5*h), sol(i-1,2:end)+(0.5*h*k1)]); 
      k3=h*dfuns{j-1}(E, [sol(i-1,1)+(0.5*h), sol(i-1,2:end)+(0.5*h*k2)]); 
      k4=h*dfuns{j-1}(E, [sol(i-1,1)+h, sol(i-1,2:end)+(h*k3)]); 
      sol(i,j)=sol(i-1,j)+((1/6)*(k1+(2*k2)+(2*k3)+k4));  
     end 
     end 
    end 
    if strcmp("heun",method) 
     if i~=1 
     for k=2:columns 
      sol(i,k)=sol(i-1,k)+(h/2)*((dfuns{k-1}(E, sol(i-1,1:end)))+(dfuns{k-1}(E, [sol(i,1),heun]))); 
     end 
     end 
    end  
    end 
end 

當我使用功能對單個普通微分方程，所述RK4方法是如預期的最好的，但是當我跑的代碼爲一對夫婦系統的微分方程，RK4是最差的，我一直在檢查和檢查，我不知道我做錯了什麼。

下面的代碼是如何調用該函數

F{1} = @(e, y) 0.6*y(3); 
F{2} = @(e, y) -0.6*y(3)+0.001407*y(4)*y(3); 
F{3} = @(e, y) -0.001407*y(4)*y(3); 

steps = 24; 

sol1 = coupled_ode(0,F,steps,0,24,[0 5 995],"euler"); 
sol2 = coupled_ode(0,F,steps,0,24,[0 5 995],"heun"); 
sol3 = coupled_ode(0,F,steps,0,24,[0 5 995],"rk4"); 

plot(sol1(:,1),sol1(:,4),sol2(:,1),sol2(:,4),sol3(:,1),sol3(:,4)); 
legend("Euler", "Heun", "RK4"); 

Answer 1

細心的例子：有是一個數太多h在RK4的formulæ：

k2 = h*dfuns{ [...] +(0.5*h*k1)]); 
k3 = h*dfuns{ [...] +(0.5*h*k2]); 

應該

k2 = h*dfuns{ [...] +(0.5*k1)]); 
k3 = h*dfuns{ [...] +(0.5*k2]); 

（最後h被刪除）。

但是，這對您提供的示例沒有影響，因爲h=1那裏。

但除了那個小錯誤之外，我認爲你實際上並沒有做錯什麼。

如果我繪製由更先進，適應性4ᵗʰ/5ᵗʰ爲了RK所產生的解決方案ode45實施：

F{1} = @(e,y) +0.6*y(3); 
F{2} = @(e,y) -0.6*y(3) + 0.001407*y(4)*y(3); 
F{3} = @(e,y)   -0.001407*y(4)*y(3); 

tend = 24; 
steps = 24; 
y0 = [0 5 995]; 
plotN = 2; 

sol1 = coupled_ode(0,F, steps, 0,tend, y0, 'euler'); 
sol2 = coupled_ode(0,F, steps, 0,tend, y0, 'heun'); 
sol3 = coupled_ode(0,F, steps, 0,tend, y0, 'rk4'); 

figure(1), clf, hold on 
plot(sol1(:,1), sol1(:,plotN+1),... 
    sol2(:,1), sol2(:,plotN+1),... 
    sol3(:,1), sol3(:,plotN+1)); 

% New solution, generated by ODE45 
opts = odeset('AbsTol', 1e-12, 'RelTol', 1e-12); 

fcn = @(t,y) [F{1}(0,[0; y]) 
       F{2}(0,[0; y]) 
       F{3}(0,[0; y])]; 
[t,solN] = ode45(fcn, [0 tend], y0, opts);  
plot(t, solN(:,plotN)) 

legend('Euler', 'Heun', 'RK4', 'ODE45'); 
xlabel('t');  

然後我們有更多的東西可信來比較。

現在，普通而簡單的RK4確實可怕這種孤立的情況下執行：

但是，如果我只是翻轉的最後一項的跡象在最近兩個功能：

%      ± 
F{2} = @(e,y) +0.6*y(3) - 0.001407*y(4)*y(3); 
F{3} = @(e,y)   +0.001407*y(4)*y(3); 

然後我們得到這個：

RK4對你的情況表現不好的主要原因是步長。自適應RK4/5（公差設置爲1而不是1e-12）產生的平均值δt= 0.15。這意味着基本的錯誤分析表明，對於這個特定的問題，h = 0.15是您可以採取的最大步驟，而不會引入不可接受的錯誤。

但是你採取的是h = 1，這確實給出了一個很大的累積誤差。

Heun和Euler對您的案例表現如此出色的事實恰恰是，幸運的是，正如上面的符號倒置示例所展示的那樣。

歡迎來到數學數學的世界 - 從來沒有一種方法最適合所有情況下的所有問題:)