www.gusucode.com > matlab 案例源码 matlab代码程序 > matlab/NonnegativeODESolutionExample.m
%% Nonnegative ODE Solution % This topic shows how to constrain the solution of an ODE to be % nonnegative. Imposing nonnegativity is not always trivial, but sometimes % it is necessary due to the physical interpretation of the equations or % due to the nature of the solution. You should only impose this constraint % on the solution when necessary, such as in cases where the integration % fails without it, or where the solution would be inapplicable. % % If certain components of the solution must be nonnegative, then use % |odeset| to set the |NonNegative| option for the indices of these % components. This option is not available for |ode23s|, |ode15i|, or for % implicit solvers (|ode15s|, |ode23t|, |ode23tb|) applied to problems with % a mass matrix. %% Example: Absolute Value Function % Consider the initial value problem % % $$y' = -|y| ,$$ % % solved on the interval $[0,40]$ with the initial condition $y(0) = 1$. % The solution of this ODE decays to zero. If the solver produces a % negative solution value, then it begins to track the solution of the ODE % through this value, and the computation eventually fails as the % calculated solution diverges to $- \infty$. Using the |NonNegative| % option prevents this integration failure. %% % Compare the analytic solution of $y(t) = e^{-t}$ to a solution of the ODE % using |ode45| with no extra options, and to one with the |NonNegative| % option set. ode = @(t,y) -abs(y); % Standard solution with |ode45| options1 = odeset('Refine',1); [t0,y0] = ode45(ode,[0 40],1,options1); % Solution with nonnegative constraint options2 = odeset(options1,'NonNegative',1); [t1,y1] = ode45(ode,[0 40],1,options2); % Analytic solution t = linspace(0,40,1000); y = exp(-t); %% % Plot the three solutions for comparison. Imposing nonnegativity is % crucial to keep the solution from veering off toward $- \infty$. plot(t,y,'b-',t0,y0,'ro',t1,y1,'k*'); legend('Exact solution','No constraints','Nonnegativity', ... 'Location','SouthWest') %% Example: The Knee Problem % Another example of a problem that requires a nonnegative solution is the % _knee problem_ coded in the example file |kneeode|. The equation is % % $$\epsilon y' = \left(1-x\right)y - y^2,$$ % % solved on the interval $[0,2]$ with the initial condition $y(0) = 1$. The % parameter $\epsilon$ generally is taken to satisfy $0 < \epsilon \ll 1$, % and this problem uses $\epsilon = 1 \times 10^{-6}$. The solution to this % ODE approaches $y = 1-x$ for $x<1$ and $y=0$ for $x>1$. However, % computing the numerical solution with default tolerances shows that the % solution follows the $y=1-x$ isocline for the whole interval of % integration. Imposing nonnegativity constraints results in the correct % solution. %% % Solve the knee problem with and without nonnegativity constraints. epsilon = 1e-6; y0 = 1; xspan = [0 2]; odefcn = @(x,y,epsilon) ((1-x)*y - y^2)/epsilon; % Solve without imposing constraints [x1,y1] = ode15s(@(x,y) odefcn(x,y,epsilon), xspan, y0); % Impose a nonnegativity constraint options = odeset('NonNegative',1); [x2,y2] = ode15s(@(x,y) odefcn(x,y,epsilon), xspan, y0, options); %% % Plot the solutions for comparison. plot(x1,y1,'ro',x2,y2,'b*') axis([0,2,-1,1]) title('The "knee problem"') legend('No constraints','Non-negativity') xlabel('x') ylabel('y') %% References % [1] Shampine, L.F., S. Thompson, J.A. Kierzenka, and G.D. Byrne, % "Non-negative solutions of ODEs," _Applied Mathematics and Computation_ % Vol. 170, 2005, pp. 556-569.