How to adapt RK4 ODE Solver from first order to 2nd order

Question

I'd like to simulate this system of equations in a matlab program.

I'm need help to set up the code. I have a Matlab function for doing Runge-Kutta4k approximation for first-order ODE's, and I want to adapt it to work for second-order ODE's.

Here's some code :

RK4 function:

    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

Main :

<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() function :

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

Please help me set up the code for my problem :) I know I have to split my 3 2nd Order ODE into 6 first ODE but I can't figure out how. Thanks in advance for the help ! Regards.

Please start by correcting the obvious errors like matching function calls to function interfaces and removing hard-coded constants where there should be none, such as in the state space dimension not always being 2. — Lutz Lehmann, May 17 '21 at 11:13
And why do you proclaim problems with the first-order system when in the code you already have implemented it? Please indicate what of the code fragments are from you or at least are completely digested by you. — Lutz Lehmann, May 17 '21 at 13:43
They're all from me except dxdt() which I saw on Internet. I proclaim problems with adapting my first order ODE solver (which works) into a second order ODE solver. — Antoine Coty, May 17 '21 at 18:26
See https://stackoverflow.com/questions/60405185/is-there-a-better-way for applying RK4 in a partitioned way to second order DE systems. But your existing framework, suitably corrected, should also work. — Lutz Lehmann, May 17 '21 at 18:56
I've looked at the link but for me (I'm a beginner ^^) it's too far away from matlab programmation langage to apply it to my program... Could you help me set it up with my method ? — Antoine Coty, May 17 '21 at 22:49

score 0 · Accepted Answer · answered May 18 '21 at 08:11

Cleaning up the code you have, for the first-order system that you have correctly given, gives a script (decompose in single modules if you prefer that or for reuse)

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

Thanks a lot ! I like the way you made changes in my code (for example passing tarray from which we can find h and tfinal instead of passing them both) — Antoine Coty, May 18 '21 at 09:42

How to adapt RK4 ODE Solver from first order to 2nd order

1 Answers1