4th order Runge Kutta Matlab with 2 2nd order ode

[noelll]

Homework Statement

Hi There!
MX''=Fn(sin θ - uCos θ )
MZ''=Fn(cos θ + uSin θ ) - Mg
Fn,M,θ,u is constant
fn/M = 0.866

M = 6000
θ = 30
u = 0.5774

i split my motion equation into 2 individual 1st ode,
X' = Vx
Z' = Vz

Vx'=[fn*(sin θ - uCos θ )]/M
Vz'={[fn(cos θ + uSin θ )]/M} - g

After trying many tutorials i still don't understand how to implement my motion equation into matlab,
i tried to do the following for Vz' motion but it seems wrong, my right hand side of the equation does not have t or y. please help !

clc; % Clears the screen
clear all;
h=0.1; % step size
thete=30;
g=9.81;
x = 0:h:3; % Calculates upto y(3)
y = zeros(1,length(x));

y(1) = 5; % initial condition

F_z = @(t,y) 0.866*(cos(thete)+0.5774*(sin(thete)))-g; % function

for i=1:(length(x)-1) % calculation loop
k_1 = F_z(x(i),y(i));
k_2 = F_z(x(i)+0.5*h,y(i)+0.5*h*k_1);
k_3 = F_z((x(i)+0.5*h),(y(i)+0.5*h*k_2));
k_4 = F_z((x(i)+h),(y(i)+k_3*h));
y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h; % main equation

end

hold on %Plot graph
plot(x,y,'+-', 'Linewidth', 1.5, 'color', 'blue')
xlabel('x')
ylabel('y')
legend('RK4')

Thank you

 

[DrClaude]

It shouldn't be too hard to find the analytical solution to the ODEs and compare to what you are calculating. Your RK algorithm appears correct.

Next time, please use code tags around your code to make it easier to read, eg,

F_z = @(t,y) 0.866*(cos(thete)+0.5774*(sin(thete)))-g; % function

gives

F_z = @(t,y) 0.866*(cos(thete)+0.5774*(sin(thete)))-g; % function

 

[noelll]

Hi thanks so much for your reply,

However my knowledge on MATLAB functions is still limited , base on my four 1st order ode

X' = Vx
Z' = Vz

Vx'=[fn*(sin θ - uCos θ )]/M
Vz'={[fn(cos θ + uSin θ )]/M} - g

Do I need to write a function for each ode? How do I write it?

Currently I only have one function which is

F_z = @(t,y) 0.866*(cos(thete)+0.5774*(sin(thete)))-g; % function

Please help

Thank you

 

[DrClaude]

The Matlab documentation available online is very good. For functions, see http://www.mathworks.com/help/matlab/ref/function.html

 

[noelll]

The Matlab documentation available online is very good. For functions, see http://www.mathworks.com/help/matlab/ref/function.html

Hi thanks for your reply,

clc;                                               % Clears the screen
clear all;

thete=30;
g= 9.81;                                             %Constant

h=0.001;                                              % step size
tfinal = 3;
N=ceil(tfinal/h);                                      %ceil is round up                              

t(1) = 0;% initial condition
Vx(1)=0;%initial accleration
X(1)=0;
Vz(1)=0;
Z(1)=0;%initial velocity

F_X  = @(t,X,Vx) Vx;
F_Z  = @(t,Z,Vz) Vz;
F_Vx = @(t,X,Vx)(0.866*(sin(thete)-0.5774*(cos(thete))));
F_Vz = @(t,Z,Vz)(0.866*(cos(thete)+0.5774*(sin(thete)))-g);for i=1:N-1;                                         % calculation loop
   
    %update time
    t(i+1)=t(i)+h;
    %update motions main equation
    k_1 = F_X (t(i)      ,X(i)          ,Vx(i));
    L_1 = F_Vx(t(i)      ,X(i)          ,Vx(i));
    k_2 = F_X (t(i)+0.5*h,X(i)+0.5*h*k_1,Vx(i)+0.5*h*L_1);
    L_2 = F_Vx(t(i)+0.5*h,X(i)+0.5*h*k_1,Vx(i)+0.5*h*L_1);
    k_3 = F_X (t(i)+0.5*h,X(i)+0.5*h*k_2,Vx(i)+0.5*h*L_2);
    L_3 = F_Vx(t(i)+0.5*h,X(i)+0.5*h*k_2,Vx(i)+0.5*h*L_2);
    k_4 = F_X (t(i)+h    ,X(i)+h    *k_3,Vx(i)+    h*L_3);
    L_4 = F_Vx(t(i)+h    ,X(i)+h    *k_3,Vx(i)+    h*L_3);

   
    k_11 = F_Z(t(i)      ,Z(i)          ,Vz(i));
    L_11 = F_Vz(t(i)     ,Z(i)          ,Vz(i));
    k_22 = F_Z(t(i)+0.5*h,Z(i)+0.5*h*k_11,Vz(i)+0.5*h*L_11);
    L_22 = F_Vz(t(i)+0.5*h,Z(i)+0.5*h*k_11,Vz(i)+0.5*h*L_11);
    k_33 = F_Z(t(i)+0.5*h,Z(i)+0.5*h*k_22,Vz(i)+0.5*h*L_22);
    L_33 = F_Vz(t(i)+0.5*h,Z(i)+0.5*h*k_22,Vz(i)+0.5*h*L_22);
    k_44 = F_Z(t(i)+    h,Z(i)+ h*k_33,Vz(i)+   h*L_33);
    L_44 = F_Vz(t(i)+   h,Z(i)+ h*k_33,Vz(i)+   h*L_33);    X(i+1) = X(i) + (h/6)*(k_1+2*k_2+2*k_3+k_4);
    Vx(i+1) = X(i) + (h/6)*(L_1+2*L_2+2*L_3+L_4);
    Z(i+1) = Z(i) + (h/6)*(k_11+2*k_22+2*k_33+k_44);
    Vz(i+1) = Z(i) + (h/6)*(L_11+2*L_22+2*L_33+L_44);
end
plot(t,X,'--', 'Linewidth', 1.5, 'color', 'red')
hold on
plot(t,Vx,'--', 'Linewidth', 1.5, 'color', 'black')
hold on
plot(t,Z,'--', 'Linewidth', 1.5, 'color', 'blue')
hold on
plot(t,Vz,'--', 'Linewidth', 1.5, 'color', 'yellow')
xlabel('time')
ylabel('height')
legend('X','Vx','Z','Vz')

i did it , however now the problem is on my graph when i generate it, i am doing a trajectory simulation of a free-fall lifeboat falling from a height of 5m above water level, the current graph i generate out starts from 0

and i should get something like the picture i attached
can you help me see is there anything wrong with my code?

thank you