4th ODE runge kutta (hiemenz equation)

AI Thread Summary
The discussion revolves around implementing a MATLAB code using the Runge-Kutta method to solve the Hiemenz equation, expressed as F''' + FF'' + 1 - F'^2 = 0, with specified boundary conditions. The user has attempted various Runge-Kutta methods (RK4, RK5, and RK Fehlberg), but the results do not align with expected values. They assumed an initial value for F'' as 1 and are seeking assistance to correct their approach. Another participant suggests exploring the BVP4C function for boundary value problems and references a thesis that successfully used shooting methods for similar equations. The user expresses interest in the thesis and its source code, indicating a desire to compare and troubleshoot their implementation.
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
Views
1K
Wrong solution order using Runge Kutta 4
Replies
15
Views
3K
Problems with Runge-Kutta method
Replies
5
Views
2K
MATLAB Resolving MATLAB Plotting Issue with Velocity & Position
Replies
4
Views
2K
MATLAB Does anyone had this plotting problem with MATLAB?
Replies
4
Views
2K
Runge-Kutta 4 order method. Severe 161: Array bounds exceeded.
Replies
3
Views
3K
Mathematica Mathematica: Solving system of three equations
Replies
2
Views
3K
I Two-Mass Oscillator: Plotting Amplitudes Over Frequency in Hertz
Replies
10
Views
2K
MATLAB Solving Orbital Motion Equations with RK4 Using MATLAB
Replies
4
Views
3K
MATLAB How can I solve a 2nd order ODE in matrix form using Matlab?
Replies
2
Views
4K

Hot Threads

Recent Insights

Back
Top