## Methods for Systems of ODEs

The extension of all of the above methods to higher-order equations or to K coupled equations, u»=fW (I, t), k = 1; ..., K; ! = («W , ... 

is straightforward. Note that any higher-order ODE can be put in the form . For example, «'=f(u, «, u, t) is written in this form as follows: u(1) =u, u(1) = u(2), «(2) = «(3), «l(3) = f(u(^1), u(2), tP\ t). The extension of, for example, the RK method  to  is, for k = 1,..., K

= Tf ®(Un + 2 U2; tn + 1 T^j R» = Tf (k)(Un + tn + t)

where Un, R1, etc., are defined analogously to « in . As for a single ODE, an important consideration when solving systems of ODEs is the stability of the numerical method. Here the small deviation Su from the solution « of  satisfies an equation analogous to  where the coefficient 8uf is replaced with the matrix whose (j, k)th entry is 8f(j)/0«(k). By diagonalizing this matrix (which is possible in the generic case), one finds that the component of Su 'aligned along' the Ith eigenvector of this matrix satisfies , with A = A/ being the corresponding eigenvalue. Thus, the stability analysis for systems of ODEs reduces to that for a single ODE, where A may be complex even though the original equations are real-valued.

In this regard, systems which have purely imaginary eigenvalues when linearized near a stable equilibrium require special attention. An example of such a system that is typical in the study of population dynamics is the Lotka-Volterra predator-prey model:

where u and v are the populations of the prey and predators and a, b, c, and d are positive constants. Linearization of  near one of its equilibrium points, (« = c/d, v = a/b), yields a pair of imaginary eigenvalues iiy/Oc. Then the stability results presented above show that the explicit Euler method  is unstable for the Lotka-Volterra model; see Figure 1. (The second-order Heun method ,  and even the fifth-order RK method used by MATLAB's ode45 are also unstable, but their instability is much weaker than that of .) Moreover, the implicit Euler method , which is stable for this model, also produces a blatantly incorrect solution. Indeed, the point t • i^fac is strictly inside its stability region (see Figure 1 ), and hence the numerical solution obtained by this method will spiral toward the equilibrium point. However, the exact analytical solution of  is known to orbit periodically about that point. As we show below, the following first-order-accurate modification of  yields a solution of  that stays near the exact periodic solution for all times:

The reason is that the stability region for this method can be shown to be the same as that of the leap-frog method , that is, a segment along the imaginary (TA)-axis (see Figure 1). Since, as stated above, eigenvalues of the linearization of  at any point sufficiently near the equilibrium are almost purely imaginary, then the corresponding value TA is on the boundary of the stability region of method . Then small deviations from the exact solution (as well as from the equilibrium) will

0.02

0.019 Time

Figure 3 Evolution of quantity Q = vuc exp[—v — cu], which is conserved for the exact solution of the Lotka-Volterra model with a = b = 1 and d = c. The time is measured in nondimensional units, as in . The initial condition is u0 = v0 = 2, and c = 2. The results are shown for methods  ('simple' Euler),  ('symplectic' Euler), and  (fourth-order RK). Values of the step size are indicated in the text. The smaller the deviation of Q(t) from Q(0), the better the method performs.

Time

Figure 3 Evolution of quantity Q = vuc exp[—v — cu], which is conserved for the exact solution of the Lotka-Volterra model with a = b = 1 and d = c. The time is measured in nondimensional units, as in . The initial condition is u0 = v0 = 2, and c = 2. The results are shown for methods  ('simple' Euler),  ('symplectic' Euler), and  (fourth-order RK). Values of the step size are indicated in the text. The smaller the deviation of Q(t) from Q(0), the better the method performs.

neither grow nor decay with time, but will oscillate around that solution. Thus, the important lesson to be drawn from this example is that the choice of the numerical method depends on the eigenvalues of the linearized system in question. That is, methods that work for systems with purely imaginary eigenvalues may not work for systems with real (and negative) eigenvalues, and vice versa. In passing, let us also mention that method  and the leap-frog method  are members of the family of symplectic methods which are designed specifically for stable numerical integration of models exhibiting oscillatory behavior without dissipation.

When dealing with an arbitrary system of ODEs, it is usually not possible to establish beforehand whether the eigenvalues of its linearization are purely imaginary or not. In such cases, one should use the RK methods (of order 4 or higher, if studying long-term evolution). Indeed, if any of the eigenvalues are not purely imaginary, one needs a method whose stability region extends into the left half of the complex (rA)-plane. On the other hand, if there are purely imaginary eigenvalues, the corresponding values of rA must be on the boundary of the stability region, in order for the numerical solution to be close to the exact solution for all times (see above). From Figure 2, one sees that the fourth- and fifth-order RK methods satisfy both of these requirements, since their right boundary is following a segment along the imaginary axis very closely (although, as one can show, not exactly). To illustrate the applicability of high-order RK methods to model oscillatory behavior, we follow the evolution of a numerically computed quantity that is conserved by the exact solution of . In Figure 3, we plot this quantity obtained with: the 'simple' Euler 40 60

Time

Figure 4 The solution of the stiff model  with the parameters indicated in the text.

40 60

Time

Figure 4 The solution of the stiff model  with the parameters indicated in the text.

method , the 'symplectic' Euler method , and the RK method . For a fair comparison, we used r = 10—3 for each of the first-order methods  and  and r = 2?10— for the fourth-order method . This choice of r makes it clear that the observed poor performance of the 'simple' Euler method is not due to its lower accuracy compared to the RK method, but due to the inferior stability of the former method. Indeed, the 'symplectic' Euler method, whose accuracy is the same as that of 'simple' Euler method but the stability region is 'just right' for this problem, stays close to the exact solution (as can be shown, even for longer times than the RK method). 