龙格 - 库塔用于耦合常微分方程

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对您的案例表现如此出色的事实恰恰是，幸运的是，正如上面的符号倒置示例所展示的那样。

欢迎来到数学数学的世界 - 从来没有一种方法最适合所有情况下的所有问题:)