gusucode.com > demos工具箱matlab源码程序 > demos/fsbvp.m
function fsbvp(solver) %FSBVP Continuation by varying an end point. % Falkner-Skan BVPs arise from similarity solutions of viscous, % incompressible, laminar flow over a flat plate. An example is % f''' + f*f'' + beta*(1-(f')^2) = 0 % with f(0) = 0, f'(0) = 0, f'(infinity) = 1 and beta = 0.5. % % The BVP is solved by imposing the boundary condition at infinity % at a finite point 'infinity'. Continuation in this end point is % used to get convergence for large values of 'infinity' and to gain % confidence from consistent results that 'infinity' is big enough. % The solution for one value of 'infinity' is extended to a guess for % a bigger 'infinity' using BVPXTEND. % % By default, this example uses the BVP4C solver. Use syntax % FSBVP('bvp5c') to solve this problem with the BVP5C solver. % % See also BVP4C, BVP5C, BVPINIT, BVPXTEND, FUNCTION_HANDLE. % Jacek Kierzenka and Lawrence F. Shampine % Copyright 1984-2014 The MathWorks, Inc. if nargin < 1 solver = 'bvp4c'; end bvpsolver = fcnchk(solver); % Problem parameter, shared with nested functions. beta = 0.5; infinity = 3; maxinfinity = 6; % This constant guess satisfying the boundary conditions % is good enough to get convergence when 'infinity' = 3. solinit = bvpinit(linspace(0,infinity,5),[0 0 1]); sol = bvpsolver(@fsode,@fsbc,solinit); eta = sol.x; f = sol.y; % Reference solution from T. Cebeci and H.B. Keller, Shooting and parallel % shooting methods for solving the Falkner-Skan boundary-layer equation, J. % Comp. Phy., 7 (1971) p. 289-300. fprintf('\n'); fprintf('Cebeci & Keller report that f''''(0) = 0.92768.\n') fprintf('Value computed using infinity = %g is %7.5f.\n',infinity,f(3,1)) figure plot(eta,f(2,:),eta(end),f(2,end),'o'); axis([0 maxinfinity 0 1.4]); title(['Falkner-Skan equation, positive wall shear, \beta = ',... sprintf('%.1f.',beta)]) xlabel('\eta') ylabel('df/d\eta') hold on drawnow shg for Bnew = infinity+1:maxinfinity solinit = bvpxtend(sol,Bnew); % Extend the solution to Bnew. sol = bvpsolver(@fsode,@fsbc,solinit); eta = sol.x; f = sol.y; fprintf('Value computed using infinity = %g is %7.5f.\n',Bnew,f(3,1)) plot(eta,f(2,:),eta(end),f(2,end),'o'); drawnow end hold off % ----------------------------------------------------------------------- % Nested functions -- beta is provided by the outer function. % function dfdeta = fsode(eta,f) dfdeta = [ f(2) f(3) -f(1)*f(3) - beta*(1 - f(2)^2) ]; end % ----------------------------------------------------------------------- function res = fsbc(f0,finf) res = [f0(1) f0(2) finf(2) - 1]; end % ----------------------------------------------------------------------- end % fsbvp