耦合 ODE 的 Runge-kutta [英] Runge-kutta for coupled ODEs

问题描述

我正在用 Octave 构建一个函数，它可以解决 N 类型的耦合常微分方程:

I’m building a function in Octave that can solve N coupled ordinary differential equation of the type:

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).

With any of these three methods (Euler, Heun and 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最差，我一直在检查和检查我不知道我做错了什么.

When I use the function for a single ordinary differential equation, the RK4 method is the best as expected, but when I ran the code for a couple system of differential equation, RK4 is the worst, I've been checking and checking and I don't know what I am doing wrong.

以下代码是如何调用该函数的示例

The following code is an example of how to call the function

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");

推荐答案

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

Careful: there's a few too many h's in the 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 已删除).

(last h's removed).

但是，这与您提供的示例没有区别，因为 h=1 在那里.

However, this makes no difference for the example that you provided, since h=1 there.

但是除了那个小错误之外，我认为您实际上并没有做错任何事情.

But other than that little bug, I don't think you're actually doing anything wrong.

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

If I plot the solution generated by the more advanced, adaptive 4ᵗʰ/5ᵗʰ order RK implemented in 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');    

然后我们有更可信的东西可以比较.

Then we have something more believable to compare to.

现在，简单明了的 RK4 确实在这种孤立的情况下表现得非常糟糕:

Now, plain-and-simple RK4 indeed performs terribly for this isolated case:

但是，如果我简单地翻转最后两个函数中最后一项的符号:

However, if I simply flip the signs of the last term in the last two functions:

%                       ± 
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 是您可以采取的最大步骤，而不会引入不可接受的错误.

The main reason RK4 performs badly for your case is because of the step size. The adaptive RK4/5 (with a tolerance set to 1 instead of 1e-12 as above) produces an average δt = 0.15. This means that basic error analysis has indicated that for this particular problem, h = 0.15 is the largest step you can take without introducing unacceptable error.

但是您采用了 h = 1，这确实会产生很大的累积误差.

But you were taking h = 1, which then indeed gives a large accumulated error.

Heun 和 Euler 在您的案例中表现如此出色的事实只是运气好，如上面的符号反转示例所示.

The fact that Heun and Euler perform so well for your case is, well, just plain luck, as demonstrated by the sign inversion example above.

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

Welcome to the world of numerical mathematics - there never is 1 method that's best for all problems under all circumstances :)

这篇关于耦合 ODE 的 Runge-kutta的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！