# 4th ODE runge kutta (hiemenz equation)

Tags:
1. Dec 17, 2015

### seyfi

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.

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. Dec 19, 2015

### BvU

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.

 wolframalpha kind of likes the y''(0)=1.2326 too...

Last edited by a moderator: May 7, 2017
3. Dec 20, 2015

### seyfi

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