MX''=Fn(cosΦ−usinΦ)(adsbygoogle = window.adsbygoogle || []).push({});

MZ''=Fn(sinΦ+ucosΦ)−Mg

MΦ''=Fn(Bxx+uBz)

I tried using Runge-Kutta methods to approximate motion equations in matlab but it turn out wrong.

I also tired finding and researching forums and web for solution but to no avail.

Fn,M,θ,u is constant fn/M = 0.866

it seems that i did my coupled runge kutta wrongly. my phi and omega does not change with time. please help

Code (Matlab M):

clc; % Clears the screen

clear all;

thete=30;

g= 9.81; %Constant

h=0.01; % 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)=20;%initial velocity

fn = 0.866;

Bz = 0.35;

Phi(1)= 30;

W(1)= 0;

%Bxx = ((Z-Z(1))+5.8*cos(thete)+Bz*(sin(Phi)-sin(thete)));

Bxx = 0.5;

%sliding phase

F_X = @(t,X,Vx,Phi) Vx;

F_Z = @(t,Z,Vz,Phi) Vz;

F_Phi= @(t,Phi,W) W;

F_Vx = @(t,X,Vx,Phi)(fn*(cos(Phi)-0.05*(sin(Phi))));

F_Vz = @(t,Z,Vz,Phi)(fn*(sin(Phi)+0.05*(cos(Phi)))-g);

F_W = @(t,Phi,W)((fn*(Bxx+Bz*0.05))/20000);

for i=1:N; % calculation loop

%update time

t(i+1)=t(i)+h;

%update motions main equation

%rotation phase

k_1 = F_X (t(i) ,X(i) ,Vx(i) ,Phi(i));

L_1 = F_Vx(t(i) ,X(i) ,Vx(i) ,Phi(i));

k_2 = F_X (t(i)+0.5*h,X(i)+0.5*h*k_1,Vx(i)+0.5*h*L_1,Phi(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,Phi(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,Phi(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,Phi(i)+0.5*h*L_2);

k_4 = F_X (t(i)+h ,X(i)+h *k_3,Vx(i)+ h*L_3,Phi(i)+h*L_3);

L_4 = F_Vx(t(i)+h ,X(i)+h *k_3,Vx(i)+ h*L_3,Phi(i)+h*L_3);

k_11 = F_Z (t(i) ,Z(i) ,Vz(i) ,Phi(i));

L_11 = F_Vz(t(i) ,Z(i) ,Vz(i) ,Phi(i));

k_22 = F_Z (t(i)+0.5*h,Z(i)+0.5*h*k_11,Vz(i)+0.5*h*L_11,Phi(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,Phi(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,Phi(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,Phi(i)+0.5*h*L_22);

k_44 = F_Z (t(i)+ h,Z(i)+ h*k_33,Vz(i)+ h*L_33,Phi(i)+h*L_33);

L_44 = F_Vz(t(i)+ h,Z(i)+ h*k_33,Vz(i)+ h*L_33,Phi(i)+h*L_33);

k_5 = F_Phi (t(i) ,Phi(i) ,W(i));

L_5 = F_W (t(i) ,Phi(i) ,W(i));

k_6 = F_Phi (t(i)+0.5*h,Phi(i)+0.5*h*k_5,W(i)+0.5*h*L_5);

L_6 = F_W (t(i)+0.5*h,Phi(i)+0.5*h*k_5,W(i)+0.5*h*L_5);

k_7 = F_Phi (t(i)+0.5*h,Phi(i)+0.5*h*k_6,W(i)+0.5*h*L_6);

L_7 = F_W (t(i)+0.5*h,Phi(i)+0.5*h*k_6,W(i)+0.5*h*L_6);

k_8 = F_Phi (t(i)+ h,Phi(i)+ h*k_7,W(i)+ h*L_7);

L_8 = F_W (t(i)+ h,Phi(i)+ h*k_7,W(i)+ h*L_7);

X(i+1) = X(i) + (h/6)*(k_1+2*k_2+2*k_3+k_4);

Vx(i+1) = Vx(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) = Vz(i) + (h/6)*(L_11+2*L_22+2*L_33+L_44);

Phi(i+1)= Phi(i) + (h/6)*(k_5+2*k_6+2*k_7+k_8);

W(i+1) = W(i) + (h/6)*(L_5+2*L_6+2*L_7+L_8);

end

figure

plot(X,Z,'--', 'Linewidth', 1.5, 'color', 'red')

xlabel('time')

ylabel('height')

legend('position')

figure

plot(Phi,W,'--', 'Linewidth', 1.5, 'color', 'blue')

xlabel('time')

ylabel('height')

legend('rotation')

figure

plot(t,Vx,'--', 'Linewidth', 1.5, 'color', 'black')

hold on

plot(t,Vz,'--', 'Linewidth', 1.5, 'color', 'yellow')

hold on

plot(t,W,'--','Linewidth',1.5,'color','green')

xlabel('time')

ylabel('height')

legend('Vx','Vz','W')

fprintf('time %.3f\n x %.3f\n z %.3f\n ',t,X,Z);

fprintf('W %.3f\n',W);

fprintf('Phi %.3f\n',Phi);

**Physics Forums | Science Articles, Homework Help, Discussion**

Join Physics Forums Today!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

# I 4th order Runge-Kutta with system of coupled 2nd order ode motion equations

Have something to add?

Draft saved
Draft deleted

Loading...

**Physics Forums | Science Articles, Homework Help, Discussion**