Adapt RK4 ODE Solver from first order to 2nd order

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. This is what I have for first order RK4 ( and it's not working x) ) Would anyone be able to help me find a solution ?

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

This is my main

clc;
clear;

h = 0.01; %Time step
y0 = 1; %Initial condition dx1/dt = 1 m.s
tfinal = 20;

tarray = 0:h:tfinal;
ytRK4 = RK4(@dxdt,y0,h,tfinal);
plot(tarray,ytRK4_2ndOrder,'g');

And my function

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

When I try to run my program, it says "Index exceeds the number of array elements" and "Error in dxdt (line 8)

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))];". As a beginner on Matlab, I don't really understand what that means...

Thanks in advance for the help !

[James Tursa]

This is hard to follow. Your main function calls RK4( ), which you don't show us. I don't see anywhere where you call the RK4_2ndOrder( ) function. I don't see anywhere where you calculate the ytRK4_2ndOrder variable you use for plotting. yt is dimensioned for two state variables (first dimension is 2), but you are clearly assuming 6-element state vectors in your dxdt( ) function, and the y0 initial value is only a scalar. Also your k3 line should be using k2 on the right hand side, not k1.

Can you clarify this?

[Nicolas Caro]

My bad about RK4 I forgot it :

function yt = RK4(y0,h,tfinal)
    yt = y0;
    i = 2;
    for t = h : h : tfinal
        k1 = f(t-h,yt(i-1));
        k2 = f(t - (0.5*h), yt(i-1) + 0.5*k1*h);
        k3 = f(t - (0.5*h), yt(i-1) + 0.5*k2*h);
        k4 = f(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

I'll call RK4_2ndOrder( ) in the main, like this : (my bad I forgot to copy/paste that line)

ytRK4_2ndOrder = RK4_2ndOrder(@dxdt,y0,h,tfinal);

About yt, I probably messed up because I was trying to understand the logic of other codes to apply it on mine but as you can see, it didn't work ^^.

About dxdt() I'm not sure if I'm doing things correctly, my 3 functions are :

I have to find x1, x2 and x3 and the only initial value I have id dx1/dt (and I don't know how to translate it in my code ^^).

I hope it's a bit clearer now :)

Don't hesitate to tell me if it isn't I'll try my best to explain what I think I'm doing !

[James Tursa]

I can help you set up the code for this, but first you must realize that you HAVE to have the following SIX initial conditions in order to solve this as an initial value problem: x1, x2, x3, x1dot, x2dot, x3dot. You need to recheck your assignment and find all of these initial conditions. Maybe it is just a simple statement that all of the others start at 0?

Answer 1

[James Tursa]

Try this for starters

yt = zeros(numel(y0),length(h:tfinal)); %Memory Allocation

And then maybe this is what is intended for initial values:

y0 = [0;0;0;1;0;0]; %Initial condition dx1/dt = 1 m.s

And remember to fix this line for the k2:

k3 = dxdt(t - (0.5*h), yt(:,i-1) + 0.5*k2*h);

[Nicolas Caro]

I made the changes. Thanks for explaining about y0 and how to translate my initial condition ! I get it now !

x1dot is the only initial condition I have so yeah, all others at 0 in y0.

Now I have that error : Array indices must be positive integers or logical values.

I think the error is on dxdt() function but I can't figure out a solution..

[James Tursa]

Please post your current code.

[Nicolas Caro]

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

I didn't see your message earlier sorry, I don't know if i'm doing right to simulate my system with dxdt. For now, what I have is a first order ODE solver using RK4

Main :

clc;
clear;

h = 0.01; %Time step
y0 = 1;
%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(y0,h,tfinal);

plot(tarray,ytRK4,'g');

RK4 :

function yt = RK4(y0,h,tfinal)
    yt = y0;
    i = 2;
    for t = h : h : tfinal
        k1 = f(t-h , yt(i-1));
        k2 = f(t - (0.5*h) , yt(i-1)+0.5*k1*h);
        k3 = f(t - (0.5*h) , yt(i-1)+0.5*k2*h);
        k4 = f(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

f

function grad = f(t,y)
    grad = (0.1).^t - 3*y;
end