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.