1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Homework Help: 4th order Runge Kutta Matlab with 2 2nd order ode

  1. Mar 3, 2016 #1
    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
    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


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

    Thank you
  2. jcsd
  3. Mar 3, 2016 #2


    User Avatar

    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.

    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
    Code (Matlab M):

    F_z = @(t,y) 0.866*(cos(thete)+0.5774*(sin(thete)))-g; % function
  4. Mar 3, 2016 #3
    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

    Code (Matlab M):

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

    Thank you
  5. Mar 4, 2016 #4


    User Avatar

    Staff: Mentor

  6. Mar 5, 2016 #5

    Hi thanks for your reply,

    Code (Matlab M):

    clc;                                               % Clears the screen
    clear all;

    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
    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
        %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);
    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')

    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

    Attached Files:

Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted