• Support PF! Buy your school textbooks, materials and every day products Here!

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

  • #1
13
0

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
 

Answers and Replies

  • #2
DrClaude
Mentor
7,050
3,201
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,
[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
 
  • #3
13
0
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

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

Thank you
 
  • #5
13
0
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,

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

Related Threads for: 4th order Runge Kutta Matlab with 2 2nd order ode

Replies
0
Views
1K
Replies
0
Views
5K
  • Last Post
Replies
0
Views
2K
Replies
3
Views
5K
  • Last Post
Replies
2
Views
741
Replies
2
Views
3K
  • Last Post
Replies
2
Views
743
Replies
1
Views
2K
Top