Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

4th ODE runge kutta (hiemenz equation)

  1. Dec 17, 2015 #1
    Hi all
    i want to write a matlab code by runge kutta solution for hiemenz equation.
    F''' + FF'' + 1 - F'^2 = 0
    BCs
    F(0)=F'(0)=0 and F'(inf)=1

    I have programmed for RK Fehlberg, RK4 and RK5 method but the results of these three methods are not matching with actual values.
    In the cod I defined like that
    F=f3
    F'=f2
    F''=f1
    and we dont know the initial value for F'' i assumed as 1 and started to make iteration.
    In attachment you can find all these 3 codes.

    Could you please help me?
    Thanks in advance

    Seyfi
    Code (Text):
    %Runge Kutta 4th Order

    clc;
    clear all;
    close all;
    n=0;
    nokta=1;
    h=0.01;
    f1=1;
    f2=0;
    f3=0;
    error=1;
    while error>1e-5
    close all;
    clc;
    fprintf('_________________________________________\n');
    fprintf('    t    |    f3    |    f2    |    f1\n ');
    fprintf('_________________________________________\n');
    fprintf('%6.6f | %6.6f | %6.6f | %6.6f\n ',n,f3,f2,f1);
    plot(n,f1,n,f2,n,f3);
    hold on
    f1_old=f1;
        while n<=(3-h)
        k1=((f2*f2)-(f3*f1)-1);
        l1=f1;
        m1=f2;
        k2=((f2+h*l1/2)*(f2+h*l1/2)-(f3+h*m1/2)*(f1+h*k1/2)-1);
        l2=(f1+h*k1/2);
        m2=(f2+h*l1/2);
        k3=((f2+h*l2/2)*(f2+h*l2/2)-(f3+h*m2/2)*(f1+h*k2/2)-1);
        l3=(f1+h*m2/2);
        m3=(f2+h*l2/2);
     
        k4=((f2+h*l3)*(f2+h*l3)-(f3+h*m3)*(f1+h*k3)-1);
        l4=(f1+h*k3);
        m4=(f2+h*l3);
        f1=f1+(k1+2*k2+2*k3+k4)*h/6;
        f2=f2+(l1+2*l2+2*l3+l4)*h/6;
        f3=f3+(m1+2*m2+2*m3+m4)*h/6;
     
     
        if mod(nokta,10)==0
            fprintf('%6.6f | %6.6f | %6.6f | %6.6f\n ',n+h,f3,f2,f1);
        end
        plot(n,f1,n,f2,n,f3);
        hold on
        n=n+h;
        nokta=nokta+1;
        end
     
        axis([0 10 -3 3]);
        error_f1=1-f2;
        error=max(abs(error_f1));
        if error_f1>0
        f1=f1_old*(error_f1+1);
        else
        f1=f1_old/(abs(error_f1)+1);
        end
        f2=0;
        f3=0;
        n=0;
        nokta=1;
    end
     
    Code (Text):
    %Runge Kutta 5th Order

    clc;
    clear all;
    close all;
    n=0;
    nokta=1;
    h=0.01;
    f1=1;
    f2=0;
    f3=0;
    error=1;
    while error>1e-6
    close all;
    clc;
    fprintf('_________________________________________\n');
    fprintf('    t    |    f3    |    f2    |    f1\n ');
    fprintf('_________________________________________\n');
    fprintf('%6.6f | %6.6f | %6.6f | %6.6f\n ',n,f3,f2,f1);
    plot(n,f1,n,f2,n,f3);
    hold on
    f1_old=f1;
        while n<=(3-h)
        k1=((f2*f2)-(f3*f1)-1);
        l1=f1;
        m1=f2;
        k2=h*((f2+l1/2)*(f2+l1/2)-(f3+m1/2)*(f1+k1/2)-1);
        l2=h*(f1+k1/2);
        m2=h*(f2+l1/2);
        k3=h*((f2+(3*l1+l2)/16)*(f2+(3*l1+l2)/16)-(f3+(3*m1+m2)/16)*(f1+(3*k1+k2)/16)-1);
        l3=h*(f1+(3*m1+m2)/16);
        m3=h*(f2+(3*l1+l2/16));
     
        k4=h*((f2+l3/2)*(f2+l3/2)-(f3+m3/2)*(f1+k3/2)-1);
        l4=h*(f1+k3/2);
        m4=h*(f2+l3/2);
        k5=h*((f2+(-3*l2+6*l3+9*l4)/16)*(f2+(-3*l2+6*l3+9*l4)/16)-(f3+(-3*m2+6*m3+9*m4)/16)*(f1+(-3*k2+6*k3+9*k4)/16)-1);
        l5=h*(f1+(-3*k2+6*k3+9*k4)/16);
        m5=h*(f2+(-3*l2+6*l3+9*l4)/16);
     
        k6=h*((f2+(l1+4*l2+6*l3-12*l4+8*l5)/7)*(f2+(l1+4*l2+6*l3-12*l4+8*l5)/7)-(f3+(m1+4*m2+6*m3-12*m4+8*m5)/7)*(f1+(k1+4*k2+6*k3-12*k4+8*k5)/7)-1);
        l6=h*(f1+(k1+4*k2+6*k3-12*k4+8*k5)/7);
        m6=h*(f2+(l1+4*l2+6*l3-12*l4+8*l5)/7);
        f1=f1+(7*k1+32*k3+12*k4+32*k5+7*k6)/90;
        f2=f2+(7*l1+32*l3+12*l4+32*l5+7*l6)/90;
        f3=f3+(7*m1+32*m3+12*m4+32*m5+7*m6)/90;
     
     
        if mod(nokta,10)==0
            fprintf('%6.6f | %6.6f | %6.6f | %6.6f\n ',n+h,f3,f2,f1);
        end
        plot(n,f1,n,f2,n,f3);
        hold on
        n=n+h;
        nokta=nokta+1;
        end
     
        axis([0 2 -3 3]);
        error_f1=1-f2;
        error=max(abs(error_f1));
        if error_f1>0
        f1=f1_old*(error_f1+1);
        else
        f1=f1_old/(abs(error_f1)+1);
        end
        f2=0;
        f3=0;
        n=0;
        nokta=1;
    end
    Code (Text):
    %Runge Kutta Fehlberg

    clc;
    clear all;
    close all;
    n=0;
    nokta=1;
    h=0.01;
    f1=1;
    f2=0;
    f3=0;
    error=1;
    while error>1e-6
    close all;
    clc;
    fprintf('_________________________________________\n');
    fprintf('    t    |    f3    |    f2    |    f1\n ');
    fprintf('_________________________________________\n');
    fprintf('%6.6f | %6.6f | %6.6f | %6.6f\n ',n,f3,f2,f1);
    plot(n,f1,n,f2,n,f3);
    hold on
    f1_old=f1;
        while n<=(10-h)
        k1=h*((f2*f2)-(f3*f1)-1);
        l1=h*f1;
        m1=h*f2;
        k2=h*((f2+l1/5)*(f2+l1/5)-(f3+m1/5)*(f1+k1/5)-1);
        l2=h*(f1+k1/5);
        m2=h*(f2+l1/5);
        k3=h*((f2+3*l1/40+9*l2/40)*(f2+3*l1/40+9*l2/40)-(f3+3*m1/40+9*m2/40)*(f1+3*k1/40+9*k2/40)-1);
        l3=h*(f1+3*m1/40+9*m2/40);
        m3=h*(f2+3*l1/40+9*l2/40);
     
        k4=h*((f2+3*l1/10-9*l2/10+6*l3/5)*(f2+3*l1/10-9*l2/10+6*l3/5)-(f3+3*m1/10-9*m2/10+6*m3/5)*(f1+3*k1/10-9*k2/10+6*k3/5)-1);
        l4=h*(f1+3*k1/10-9*k2/10+6*k3/5);
        m4=h*(f2+3*l1/10-9*l2/10+6*l3/5);
        k5=h*((f2+11*l1/54+5*l2/2-70*l3/27+35*l4/27)*(f2+11*l1/54+5*l2/2-70*l3/27+35*l4/27)-(f1+11*k1/54+5*k2/2-70*k3/27+35*k4/27)*(f3+11*m1/54+5*m2/2-70*m3/27+35*m4/27)-1);
        l5=h*(f1+11*k1/54+5*k2/2-70*k3/27+35*k4/27);
        m5=h*(f2+11*l1/54+5*l2/2-70*l3/27+35*l4/27);
     
        k6=h*((f2+1631*l1/55296+175*l2/512+575*l3/13824+44275*l4/110592+253*l5/4096)*(f2+1631*l1/55296+175*l2/512+575*l3/13824+44275*l4/110592+253*l5/4096)-(f1+1631*k1/55296+175*k2/512+575*k3/13824+44275*k4/110592+253*k5/4096)*(f3+1631*m1/55296+175*m2/512+575*m3/13824+44275*m4/110592+253*m5/4096)-1);
        l6=h*(f1+1631*k1/55296+175*k2/512+575*k3/13824+44275*k4/110592+253*k5/4096);
        m6=h*(f2+1631*l1/55296+175*l2/512+575*l3/13824+44275*l4/110592+253*l5/4096);
     
        f1=f1+(2825*k1/27548+18575*k3/48384+13525*k4/55296+277*k5/14336+k6/4);
        f2=f2+(2825*l1/27548+18575*l3/48384+13525*l4/55296+277*l5/14336+l6/4);
        f3=f3+(2825*m1/27548+18575*m3/48384+13525*m4/55296+277*m5/14336+m6/4);
     
     
        if mod(nokta,50)==0
            fprintf('%6.6f | %6.6f | %6.6f | %6.6f\n ',n+h,f3,f2,f1);
        end
        plot(n,f1,n,f2,n,f3);
        hold on
        n=n+h;
        nokta=nokta+1;
        end
     
        axis([0 10 -3 3]);
        error_f1=1-f2;
        error=max(abs(error_f1));
        if error_f1>0
        f1=f1_old*(error_f1+1);
        else
        f1=f1_old/(abs(error_f1)+1);
        end
        f2=0;
        f3=0;
        n=0;
        nokta=1;
    end
     
     
  2. jcsd
  3. Dec 19, 2015 #2

    BvU

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    Hi,
    I see you also posted on Matlab Central and they advise against shooting and recommend BVP4C.
    Are you re-doing the master thesis work of http://www.ewp.rpi.edu/hartford/users/papers/engr/ernesto/kaufme3/ES/update/project_draft_041508.pdf [Broken]? (summary) She claims to succeed with shooting and has all the source code in the thesis, so perhaps you can check for differences with your program.

    [edit] wolframalpha kind of likes the y''(0)=1.2326 too...
     
    Last edited by a moderator: May 7, 2017
  4. Dec 20, 2015 #3
    I run her code and its working perfectly, if I can able to reach her via mail I would thank her.
    But here the above especially Runge Kutta Fehlberg, I used it for Blasius equation and it gave me the right results.Now I modified for the F''' + FF'' + 1 - F'^2 = 0 equation but it did not work.I guess i am doing something wrong
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: 4th ODE runge kutta (hiemenz equation)
  1. ODE matlab (Replies: 2)

  2. Mathematica: ODE's (Replies: 2)

  3. Convergence of an ODE (Replies: 1)

Loading...