Stiff Equations

Models where there are two (or more) disparate timescales give rise to so-called stiff ODEs. A simple special case exhibiting such behavior is a system of two interacting species which, when considered independently of each other, evolve at their respective equilibria at very different rates. An example is a model (with exaggerated parameters) of interacting phytoplankton and zooplankton populations (P and Z), proposed by J. E. Truscott and J. Brindley in 1994:

where v, 7, and ! are some positive constant parameters of the model. For ^ = 2000, v = 0.05, 7 = 0.4,

y0=[0.5; 0.5]; °/0 initial condition for u and du/dt options=odeset(,RelTol;,10"(-3),'AbsTol',10"(-6));

Figure 5 The MATLAB codes used to obtain the solution shown in Figure 4.

! = 0.7, and the initial condition P(0) = Z(0) = 0.5, the evolutions of P and Z are shown in Figure 4. Note that when the two populations do not interact, P and Z tend to their equilibria (1 and 0) with rates proportional to the greatly different parameters ft and 7!, respectively. Mathematically, these rates are proportional to the eigenvalues of the corresponding equations linearized about any 'point' (P, Z) of the solution. In the presence of the interaction, the 'memory' about these disparate timescales is reflected by the existence of both smooth and abrupt changes of P and Z. Let us emphasize that even though those abrupt changes (with the rate proportional to ft) occur relatively rarely and 'most of the time', the evolution is smooth, the two eigenvalues A1 and A2 of the linearized system [26] have magnitudes on the orders of ft >> 1 and 7! = O(1) for all times. Then any explicit numerical method will require, as its stability condition, that t max{|A1|, |A2|} = tO (ft) be less than a constant of order one. Thus, one needs to take steps t = O(ft— 1) <<1

to resolve even the smooth motion occurring on the time-scale O ((7!) ) = O (1). A common approach to avoid such an inefficient use of computational resources is to use implicit methods (of a form more sophisticated than [12] and [14]) which, as stated earlier, have much greater stability regions than explicit methods and hence proceed in time with much larger steps. Since the implementation of such methods requires the solution of (a system of) nonlinear algebraic equations, their programming takes much more effort than that of explicit methods. Fortunately, modern computational software have built-in commands for solving (systems of) stiff ODEs. For example, the solution shown in Figure 4 was obtained by a MATLAB script shown in Figure 5. (The nonstiff solver ode4 5 for this problem would run about 50 times slower, and for larger values of ft, it would simply fail.) Even though the example above uses unrealistic values of parameters, it still illustrates the issue which occurs in many known systems of a large number of coupled equations (e.g., equations describing chemical reactions among trace gases in the atmosphere).

The periodic solutions in the plankton model [26] and Lotka-Volterra equations are different in the following regard. Any (meaningful) initial condition in the plankton model will be attracted to the unique periodic solution shown in Figure 4, while for the Lotka-Volterra model, any initial condition will be repeated periodically, that is, a continuum of periodic solutions exists. Correspondingly, the respective models are often referred to as excitable and conservative. Conservative models with four and more interacting species may have imaginary eigenvalues of disparate magnitude. (According to our definition, such systems would be also called stiff, but the convention is to use that name only for systems with Re(A) < 0.) Therefore, a suitable numerical method must have a section of its stability boundary following the imaginary axis. The corresponding MATLAB's stiff solver is ode2 3t. Its performance can be tested on the following simple model:

which has no ecological interpretation but serves to illustrate the above point. Indeed, the 'free' solution of [27] oscillates 1000 times faster than the 'forced' one, hence two disparate rates of conservative evolution are present. Other stiff and nonstiff MATLAB's solvers do not perform satisfactorily for this problem. Thus, the conclusion drawn from this example is similar to that stated earlier, namely: methods that can be applied to nonconservative systems (such as [26]) may not work for conservative ones, and vice versa.

Was this article helpful?

0 0
Project Earth Conservation

Project Earth Conservation

Get All The Support And Guidance You Need To Be A Success At Helping Save The Earth. This Book Is One Of The Most Valuable Resources In The World When It Comes To How To Recycle to Create a Better Future for Our Children.

Get My Free Ebook

Post a comment