The equation is (1+2M*eta)f'''+ 2Mf"+ f*f"- (f')^2- K1*f'= 0 ; f' is df/d(eta) 'eta' is a similarity variable
where M- is a curvature parameter, M = 1k1- is another parameter, k1= 0.1
boundary conditions are: eta=0, f(0)=0, f'(0)=1, f'(inf)=0, inf= 10
The following is the code for Blasius equation
% Fourth Runge Kutta for solving Blasius Equation
%———————————————–
%Blasius Equation : 2f"+ff'=0
%———————————————–
%Boundary Condition
%1. f'(0) = 0
%2. f'(inf) = 0
%3. f(0) = 0
%from shooting method
%4. f"(0) = 0.332
%———————————————–
%find: f,f" and f'''
%———————————————–
%Given
%eta = x
%f = y0,f'=y1,f"=y2,f'''= -(1/2)*y0*y2
%———————————————–
func1 = inline('y1','x','y0','y1','y2');
%———————————————–
func2 = inline('y2','x','y0','y1','y2');
%———————————————–
func3 = inline('-0.5*y0*y2','x','y0','y1','y2');
%———————————————–
%input: x = 0 , y0 = 0 , y1= 0
% y2 = 0.332 , total = 7 and h = 0.1
%———————————————–
x = input('\n Enter the value of x : ');
y0 = input( 'Enter the value of y0 : ');
y1 = input( 'Enter the value of y1 : ');
y2 = input( 'Enter the value of y2 : ');
total = input('Enter the value of total : ');
h = input( 'Enter the value of h : ');
fprintf('\n Solution with step size =%5.3f is:',h);
fprintf('\n x y0 y1 y2');
fprintf('\n%16.3f%16.3f%16.3f%16.3f',x,y0,y1,y2);
for i = 1:(total/h)
ak1y0 = func1(x,y0,y1,y2);
ak1y1 = func2(x,y0,y1,y2);
ak1y2 = func3(x,y0,y1,y2);
xx = x + h/2.;
yy0 = y0 + h*ak1y0/2.;
yy1 = y1 + h*ak1y1/2.;
yy2 = y2 + h*ak1y2/2.;
ak2y0 = func1(xx,yy0,yy1,yy2);
ak2y1 = func2(xx,yy0,yy1,yy2);
ak2y2 = func3(xx,yy0,yy1,yy2);
yy0 = y0 + h*ak2y0/2.;
yy1 = y1 + h*ak2y1/2.;
yy2 = y2 + h*ak2y2/2.;
ak3y0 = func1(xx,yy0,yy1,yy2);
ak3y1 = func2(xx,yy0,yy1,yy2);
ak3y2 = func3(xx,yy0,yy1,yy2);
all_x(i) = x;
all_y0(i) = y0;
all_y1(i) = y1;
all_y2(i) = y2;
xx = x + h;
yy0 = y0 + h*ak3y0;
yy1 = y1 + h*ak3y1;
yy2 = y2 + h*ak3y2;
ak4y0 = func1(xx,yy0,yy1,yy2);
ak4y1 = func2(xx,yy0,yy1,yy2);
ak4y2 = func3(xx,yy0,yy1,yy2);
y0 = y0 + (ak1y0 + 2.*ak2y0 + 2.*ak3y0 + ak4y0)*h/6.;
y1 = y1 + (ak1y1 + 2.*ak2y1 + 2.*ak3y1 + ak4y1)*h/6.;
y2 = y2 + (ak1y2 + 2.*ak2y2 + 2.*ak3y2 + ak4y2)*h/6.;
x = x + h;
fprintf('\n%16.3f%16.3f%16.3f%16.3f',x,y0,y1,y2);
end
plot(all_x, all_y0, 'k-', all_x, all_y1, 'b-', all_x, all_y2, 'g-')
Best Answer