gusucode.com > demos工具箱matlab源码程序 > demos/threebvp.m
function threebvp(solver) %THREEBVP Three-point boundary value problem % This example of multipoint BVPs comes from the study of a physiological % flow problem described in C. Lin and L. Segel, Mathematics Applied to % Deterministic Problems in the Natural Sciences, SIAM, Philadelphia, PA, % 1988. For x in [0, lambda], the equations for the problem are % % v' = (C - 1)/n % C' = (vC - min(x,1))/eta % % for known parameters n, kappa, and lambda > 1. eta = lambda^2/(n*kappa^2). % The term min(x,1) in the equation for C'(x) is not smooth at x = 1. % Lin and Segel describe this BVP as two problems, one set on [0, 1] % and the other on [1, lambda], connected by the requirement that v(x) % and C(x) be continuous at x = 1. The solution is to satisfy the boundary % conditions v(0) = 0 and C(lambda) = 1. BVP4C solves this problem as % a three-point BVP, imposing internal conditions at the interface point % x = 1. % % This example solves the problem for n = 5e-2, lambda = 2, and a range % kappa = 2,3,4,5. The solution for one value of kappa is used as guess % for the next. % % By default, this example uses the BVP4C solver. Use syntax % THREEBVP('bvp5c') to solve this problem with the BVP5C solver. % % See also BVP4C, BVP5C, BVPSET, BVPGET, BVPINIT, DEVAL, 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 n = 5e-2; % Initial mesh - duplicate the interface point x = 1. lambda = 2; xinit = [0, 0.25, 0.5, 0.75, 1, 1, 1.25, 1.5, 1.75, 2]; % Constant initial guess for the solution yinit = [1; 1]; % The initial profile sol = bvpinit(xinit,yinit); % For each kappa, the quantity of interest is the emergent osmolarity. % This example compares the value computed by BVP4C, Os = 1/v(lambda), % with Lin and Segel's asymptotic approximation. fprintf(' kappa computed Os approximate Os \n') for kappa = 2:5 eta = lambda^2/(n*kappa^2); % use in nested functions sol = bvpsolver(@f,@bc,sol); K2 = lambda*sinh(kappa/lambda)/(kappa*cosh(kappa)); approx = 1/(1 - K2); computed = 1/sol.y(1,end); fprintf(' %2i %10.3f %10.3f \n',kappa,computed,approx); end figure plot(sol.x,sol.y(1,:),sol.x,sol.y(2,:),'--') legend('v(x)','C(x)') title(['A three-point BVP solved with ',upper(func2str(bvpsolver))]) xlabel(['\lambda = ',num2str(lambda),', \kappa = ',num2str(kappa),'.']) ylabel('v and C') % ----------------------------------------------------------------------- % Nested functions -- n and eta are provided by the outer function. % function dydx = f(x,y,region) % Derivative function -- share n and eta with the outer function. dydx = zeros(2,1); dydx(1) = (y(2) - 1)/n; % The definition of C'(x) depends on the region. switch region case 1 % x in [0 1] dydx(2) = (y(1)*y(2) - x)/eta; case 2 % x in [1 lambda] dydx(2) = (y(1)*y(2) - 1)/eta; otherwise error('MATLAB:threebvp:BadRegionIndex','Incorrect region index: %d',region); end end % ----------------------------------------------------------------------- function res = bc(YL,YR) % Boundary (and internal) conditions res = [ YL(1,1) % v(0) = 0 YR(1,1) - YL(1,2) % continuity of v(x) at x = 1 YR(2,1) - YL(2,2) % continuity of C(x) at x = 1 YR(2,end) - 1 ]; % C(lambda) = 1 end % ----------------------------------------------------------------------- end % threebvp