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

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

Related Engineering and Comp Sci Homework Help News on Phys.org
DrClaude
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
Matlab:
F_z = @(t,y) 0.866*(cos(thete)+0.5774*(sin(thete)))-g; % function

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

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

Thank you

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

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

#### Attachments

• 6.2 KB Views: 467