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.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Fluid Dynamics에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by