What modifications do I need to make in the following codes for solving the boundary value problem similar to the Blasius equation using Shooting method with R-K 4 numerical analysis

1 回表示 (過去 30 日間)
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 = 1
k1- 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-')

採用された回答

Torsten
Torsten 2017 年 11 月 8 日
Use MATLAB's "bvp4c" instead of a code that you don't understand.
Best wishes
Torsten.
  1 件のコメント
naygarp
naygarp 2017 年 11 月 8 日
How do I use 'bvp4c' for solving bvp using Shooting method with RK-4 numerical method?

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeFluid Dynamics についてさらに検索

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by