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-')
0 件のコメント
採用された回答
Torsten
2017 年 11 月 8 日
Use MATLAB's "bvp4c" instead of a code that you don't understand.
Best wishes
Torsten.
その他の回答 (0 件)
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!