用于耦合ODEs的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 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 formul中有一些太多的RK4：

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方法中的各个阶段的去耦合处理(注意Heun只是Euler步骤的一个副本)将它们减少到一个顺序一的方法。

具体来说，从

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

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

考虑到这一点，改变了执行的结果

  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提供了更合理的图。

​

​

并将步骤尺寸0.2细分为120个间隔。

​

​

在这里，RK4的图形并没有发生太大的变化，而另外两个则从下面和上面向它移动。