4th ODE runge kutta (hiemenz equation)

Click For Summary
SUMMARY

The forum discussion centers on implementing a MATLAB code using the Runge-Kutta methods to solve the Hiemenz equation, specifically the third-order differential equation F''' + FF'' + 1 - F'^2 = 0 with boundary conditions F(0)=F'(0)=0 and F'(∞)=1. The user, Seyfi, has attempted the RK4, RK5, and RK Fehlberg methods but reports discrepancies between the computed results and expected values. A suggestion was made to consider using the BVP4C function instead of the shooting method, referencing a successful implementation from a master's thesis.

PREREQUISITES
  • Familiarity with MATLAB programming and syntax.
  • Understanding of numerical methods, specifically Runge-Kutta techniques.
  • Knowledge of boundary value problems (BVP) and initial value problems (IVP).
  • Basic understanding of differential equations, particularly third-order equations.
NEXT STEPS
  • Research the MATLAB BVP4C function for solving boundary value problems.
  • Study the implementation of the Runge-Kutta Fehlberg method in MATLAB for accuracy improvements.
  • Examine the master's thesis referenced for insights on successful coding practices and potential pitfalls.
  • Explore the differences between shooting methods and BVP solvers in numerical analysis.
USEFUL FOR

Mathematics and engineering students, MATLAB programmers, and researchers working on fluid dynamics or numerical methods for differential equations.

seyfi
Messages
2
Reaction score
0
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 don't 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: 
%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: 
%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: 
%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
 
Physics news on Phys.org
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 ? (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:
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
 

Similar threads

MATLAB Problem with Runge-Kutta script file
  • · Replies 4 ·
Replies
4
Views
2K
Wrong solution order using Runge Kutta 4
  • · Replies 15 ·
Replies
15
Views
3K
Problems with Runge-Kutta method
  • · Replies 5 ·
Replies
5
Views
2K
MATLAB Resolving MATLAB Plotting Issue with Velocity & Position
  • · Replies 4 ·
Replies
4
Views
2K
MATLAB Does anyone had this plotting problem with MATLAB?
  • · Replies 4 ·
Replies
4
Views
2K
Runge-Kutta 4 order method. Severe 161: Array bounds exceeded.
  • · Replies 3 ·
Replies
3
Views
3K
Mathematica Mathematica: Solving system of three equations
  • · Replies 2 ·
Replies
2
Views
3K
Undergrad Two-Mass Oscillator: Plotting Amplitudes Over Frequency in Hertz
  • · Replies 10 ·
Replies
10
Views
2K
MATLAB Solving Orbital Motion Equations with RK4 Using MATLAB
  • · Replies 4 ·
Replies
4
Views
3K
MATLAB How can I solve a 2nd order ODE in matrix form using Matlab?
  • · Replies 2 ·
Replies
2
Views
4K