用于耦合 ODE 的 Runge-kutta

Question

我正在 Octave 中构建一个函数，可以求解 N 类型的耦合常微分方程：

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

使用这三种方法中的任何一种（Euler、Heun 和 Runge-Kutta-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

注意：RK4 公式中的 h 太多了：

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 在那里。

但除了那个小错误，我不认为你实际上做错了什么。

如果我绘制由 ode45 中实现的更高级的自适应 4ᵗʰ/5ᵗʰ 阶 RK 生成的解决方案：

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 在您的案例中表现如此出色，好吧，纯属运气，如上面的符号反转示例所示。

欢迎来到数值数学世界 - 从来没有一种方法可以在所有情况下解决所有问题:)

Answer 2

除了旧答案中描述的错误外，实施中确实存在根本的方法论错误。首先，实现对于标量一阶微分方程是正确的。但是，当您尝试在耦合系统上使用它时，Runge-Kutta 方法中阶段的去耦合处理（请注意，Heun 只是 Euler 步骤的副本）将它们简化为一阶方法。

具体来说，从

开始

      k2=h*dfuns{j-1}(E, [sol(i-1,1)+(0.5*h), sol(i-1,2:end)+(0.5*h*k1)]);

在sol(i-1,2:end)中加入0.5*k1的意思是将第一阶段的斜率向量相加，而不是将相同的斜率值添加到位置向量的所有分量中。

考虑到这一点导致对实现的更改

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

可以看出，向量分量上的循环频繁出现。可以通过对耦合 ODE 系统的右侧使用矢量值函数的全矢量化来隐藏这一点。

具有这些更改的解的第二个分量的图给出了步长 1 更合理的图

并细分为 120 个间隔，步长为 0.2

其中 RK4 的图表变化不大，而其他两个从下方和上方向它移动。

用于耦合 ODE 的 Runge-kutta

Runge-kutta for coupled ODEs

matlab

octave

numerical-methods