如何将 RK4 ODE Solver 从一阶调整为二阶答案

【问题标题】：How to adapt RK4 ODE Solver from first order to 2nd order如何将 RK4 ODE Solver 从一阶调整为二阶
【发布时间】：2021-05-18 23:28:15
【问题描述】：

我想在 matlab 程序中模拟这个方程组。

我需要帮助来设置代码。我有一个 Matlab 函数，用于对一阶 ODE 进行 Runge-Kutta4k 逼近，我想对其进行调整以适用于二阶 ODE。

这里有一些代码：

RK4 功能：

    function yt = RK4_2ndOrder(dxdt,y0,h,tfinal)
    yt = zeros(2,length(h:tfinal)); %Memory Allocation
    yt(:,1) = y0; %Initial condition
    i = 2;
    for t = h : h : tfinal %RK4 loop
        k1 = dxdt(t-h,yt(:,i-1));
        k2 = dxdt(t - (0.5*h), yt(:,i-1) + 0.5*k1*h);
        k3 = dxdt(t - (0.5*h), yt(:,i-1) + 0.5*k2*h);
        k4 = dxdt(t, yt(:,i-1) + (k3 * h));
        yt(:,i) = yt(:,i-1) + (1/6 * (k1 + 2*k2 + 2*k3 + k4)* h);
        i = i + 1;
    end
end

主要：

<pre><code>clc;
clear;
h = 0.01; %Time step
y0 = [0,0,0,1,0,0]; %Initial condition dx1/dt = 1 m.s
tfinal = 20;
tarray = 0:h:tfinal;
ytRK4 = RK4(y0,h,tfinal);
ytRK4_2ndOrder = RK4_2ndOrder(@dxdt,y0,h,tfinal);
plot(tarray,ytRK4_2ndOrder,'g');
</code></pre>

dxdt() 函数：

<pre><code>function dx = dxdt(t,x)
m1 = 12000;
m2 = 10000;
m3 = 8000;
k1 = 3000;
k2 = 2400;
k3 = 1800;
dx = [x(4);x(5);x(6);-k1/m1*x(1)+k2/m1*(x(2)-x(1));k2/m2*(x(1)-x(2))+k3/m2(x(3)-x(2));k3/m3*(x(2)-x(3))];

end</pre></code>

请帮我为我的问题设置代码:) 我知道我必须将我的 3 个二阶 ODE 拆分为 6 个一阶 ODE，但我不知道如何。在此先感谢您的帮助！问候。

【问题讨论】：

请先纠正明显的错误，例如将函数调用匹配到函数接口，并删除不应存在的硬编码常量，例如状态空间维度不总是 2。
为什么你在代码中已经实现了一阶系统的问题？请指出哪些代码片段来自您或至少已被您完全消化。
除了我在 Internet 上看到的 dxdt() 之外，它们都来自我。我宣布将我的一阶 ODE 求解器（有效）调整为二阶 ODE 求解器存在问题。
请参阅stackoverflow.com/questions/60405185/is-there-a-better-way，了解以分区方式将 RK4 应用于二阶 DE 系统。但是你现有的框架，经过适当的修正，也应该可以工作。
我查看了链接，但对我来说（我是初学者 ^^）它离 matlab 编程语言太远了，无法将它应用到我的程序中......你能帮我设置一下吗符合我的方法吗？

标签： matlab

【解决方案1】：

清理您拥有的代码，对于您正确给出的一阶系统，给出一个脚本（如果您愿意或重复使用，可以分解为单个模块）

function so67567216_RK_spring_train()
  clc; clear;
  h = 0.01; %Time step
  y0 = [0,0,0,1,0,0]; %Initial condition dx1/dt = 1 m.s
  tfinal = 20;
  tarray = 0:h:tfinal;
  ytRK2 = RK2(@(t,x)dxdt(t,x),tarray,y0);
  ytRK4 = RK4(@(t,x)dxdt(t,x),tarray,y0);
  plot(tarray,ytRK4(1:3,:),'g');

  function dx = dxdt(t,x)
    m1 = 12000; m2 = 10000; m3 = 8000;
    k1 = 3000; k2 = 2400; k3 = 1800;
    dx = [x(4); x(5); x(6); -k1/m1*x(1)+k2/m1*(x(2)-x(1)); k2/m2*(x(1)-x(2))+k3/m2*(x(3)-x(2)); k3/m3*(x(2)-x(3))];
  end

  function yt = RK2(dxdt,tarray,y0)
    N = length(tarray);
    yt = zeros(length(y0),N); %Memory Allocation
    yt(:,1) = y0; %Initial condition
    for i = 2:N %RK4 loop
      h = tarray(i)-tarray(i-1);
      k1 = dxdt(tarray(i-1), yt(:,i-1));
      k2 = dxdt(tarray(i) - (0.5*h), yt(:,i-1) + 0.5*k1*h);
      yt(:,i) = yt(:,i-1) +  k2 * h;
    end%for
  end%function

  function yt = RK4(dxdt,tarray,y0)
    N = length(tarray);
    yt = zeros(length(y0),N); %Memory Allocation
    yt(:,1) = y0; %Initial condition
    for i = 2:N %RK4 loop
      h = tarray(i)-tarray(i-1);
      k1 = dxdt(tarray(i-1), yt(:,i-1));
      k2 = dxdt(tarray(i) - (0.5*h), yt(:,i-1) + 0.5*k1*h);
      k3 = dxdt(tarray(i) - (0.5*h), yt(:,i-1) + 0.5*k2*h);
      k4 = dxdt(tarray(i), yt(:,i-1) + (k3 * h));
      yt(:,i) = yt(:,i-1) + (1/6 * (k1 + 2*k2 + 2*k3 + k4)* h);
    end%for
  end%function
end%script

【讨论】：

非常感谢！我喜欢你对我的代码进行更改的方式（例如传递 tarray，我们可以从中找到 h 和 tfinal 而不是同时传递它们）