www.gusucode.com > demos工具箱matlab源码程序 > demos/kneeode.m
function kneeode %KNEEODE The "knee problem" with non-negativity constraints. % % For 0 < epsilon << 1, the solution of the initial value problem % % epsilon*y' = (1-x)*y - y^2, y(0) = 1 % % approaches null isoclines y = 1 - x and y = 0, for x < 1 and % x > 1, respectively. The numerical solution, computed with % default tolerances, follows the y = 1 - x isocline for the % whole interval of integration. Imposing non-negativity % constraints results in the correct solution. % % G. Dahlquist, L. Edsberg, G. Skollermo, G. Soderlind, Are the % Numerical Methods and Software Satisfactory for Chemical % Kinetics?, in Numerical Integration of Differential Equations % and Large Linear Systems, J. Hinze ed., Springer, Berlin, 1982, % pp. 149-164. % % See also ODE15S, ODE23T, ODE23TB, ODESET, FUNCTION_HANDLE. % Jacek Kierzenka and Lawrence F. Shampine % Copyright 1984-2014 The MathWorks, Inc. % Problem parameter epsilon = 1e-6; y0 = 1; xspan = [0, 2]; % Solve without imposing constraints options = []; [x1,y1] = ode15s(@odefcn,xspan,y0,options); % Impose non-negativity constraint options = odeset('NonNegative',1); [x2,y2] = ode15s(@odefcn,xspan,y0,options); figure plot(x1,y1,x2,y2) axis([0,2,-1,1]); title('The "knee problem"'); legend('No constraints','Non-negativity') xlabel('x'); ylabel('solution y') % ----------------------------------------------------------------------- % Nested function -- epsilon provided by the outer function. % function yp = odefcn(x,y) yp = ((1 - x)*y - y^2)/epsilon; end % ----------------------------------------------------------------------- end % kneeode