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

Tags:
1. Mar 3, 2016

### noelll

1. The problem statement, all variables and given/known data
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

2. Mar 3, 2016

### Staff: Mentor

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.

[CODE="matlab"]
F_z = @(t,y) 0.866*(cos(thete)+0.5774*(sin(thete)))-g; % function
[/CODE]
gives
Code (Matlab M):

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

3. Mar 3, 2016

### noelll

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

Code (Matlab M):

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

Thank you

4. Mar 4, 2016

5. Mar 5, 2016

### noelll

Code (Matlab M):

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

File size:
6.2 KB
Views:
64