4th ODE runge kutta (hiemenz equation)

In summary, the conversation was about writing a MATLAB code using the Runge Kutta method for solving the Hiemenz equation. The equation and boundary conditions were provided, and the person had already tried using RK Fehlberg, RK4, and RK5 methods but their results did not match the actual values. They shared their code and asked for help. The expert advised against using shooting and recommended using BVP4C instead. They also mentioned a previous similar work and suggested checking for differences in the code. WolframAlpha was also mentioned as a tool to help with finding a potential initial condition.
  • #1
seyfi
3
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
  • #2
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:
  • #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
 

FAQ: 4th ODE runge kutta (hiemenz equation)

1. What is the Runge-Kutta method for solving 4th order differential equations?

The Runge-Kutta method is a numerical method for solving ordinary differential equations (ODEs). It is a 4th order method, meaning that it uses 4 function evaluations per step to approximate the solution of the ODE. The Hiemenz equation is a specific type of 4th order ODE that can be solved using the Runge-Kutta method.

2. How does the Runge-Kutta method work?

The Runge-Kutta method works by breaking down the ODE into smaller steps and using a weighted combination of these steps to approximate the solution. It involves calculating intermediate values of the function at different points and using these values to update the solution at each step. This allows for a more accurate approximation compared to simpler methods, such as Euler's method.

3. What is the Hiemenz equation?

The Hiemenz equation is a specific type of 4th order ODE that describes the flow of an incompressible viscous fluid over a semi-infinite flat plate. It is commonly used in fluid mechanics and aerodynamics to model boundary layer flows. It takes the form of a nonlinear differential equation and can be solved using numerical methods such as the Runge-Kutta method.

4. What are the advantages of using the Runge-Kutta method for solving 4th order ODEs?

The Runge-Kutta method is a more accurate and efficient method compared to simpler methods, such as Euler's method. It can handle stiff equations, meaning that the solution changes rapidly over a small range of values. It also allows for adaptive step sizes, meaning that the step size can be adjusted to maintain a desired level of accuracy.

5. Are there any limitations to using the Runge-Kutta method for solving 4th order ODEs?

The main limitation of the Runge-Kutta method is that it can be computationally intensive for complex ODEs with a large number of variables. It also requires a known initial value, which may not always be available. In addition, it may not be suitable for stiff equations with rapidly changing solutions, as it may require a very small step size to maintain accuracy.

Similar threads

Wrong solution order using Runge Kutta 4
Replies
15
Views
2K
Problems with Runge-Kutta method
Replies
5
Views
2K
MATLAB Resolving MATLAB Plotting Issue with Velocity & Position
Replies
4
Views
2K
Runge-Kutta 4 order method. Severe 161: Array bounds exceeded.
Replies
4
Views
3K
MATLAB Does anyone had this plotting problem with MATLAB?
Replies
4
Views
2K
Mathematica Mathematica: Solving system of three equations
Replies
2
Views
2K
I Two-Mass Oscillator: Plotting Amplitudes Over Frequency in Hertz
Replies
10
Views
1K
MATLAB Coupled systems of ODE's solved by matlabs ODE45 (Runge Kutta Alg.)
Replies
0
Views
1K
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