## Implicit and Semi Implicit Methods for [1 and for Similar Coupled Systems

Let us first show how the idea of the previous paragraph can be used to obtain a numerical method for a single equation  where the coefficients D, g, and f depend only on x and t (i.e., do not depend on the unknown u). Equations with this property are referred to as linear. Generalizations to nonlinear equations and systems of linear or nonlinear coupled equations will be considered later. The second-order accurate approximations to the terms in  at the virtual node (xm, tn + (t/2)) are produced by the following discretization schemes:

g (x, t )ux ! - (¿( um+1—un—o + ¿+1 ( jn+1—um+1))

to hold just to ensure good accuracy of the method.) Let us rewrite  as u — I + E

where the evolutions associated with I and E need to be computed implicitly and explicitly, respectively. (In the case of , I — Dd2xu and E — \$ xu + f) Then a method that has the accuracy 0(r2) + O(h2) and the stability restriction of the form  is

Putting these terms together produces an unconditionally stable (i.e., with no restriction on r for the stability) method with the overall accuracy 0(r ) + <Xh ). Naturally, it is assumed that the coefficients D, g and f change little over distances of order h and times of order r. Also, the above discretizations are not unique in the sense that replacing both Dnm and 1 with (1/2) = D(xm, tn + (r/2)) (and similarly for g and f) yields another scheme with the same accuracy and the stability property as .

When generalizing the previous technique to the case where D, g, and f depend on the u, one faces the problem that the RHSs in  become nonlinear functions of the unknown variables U n+1 (see ) on the (n + 1)th time level. (If one uses D"m+1/2 instead of D"m+1 as suggested above, one still needs to find Um+ (1/2)n.) Solving systems of nonlinear algebraic equations requires considerably more effort than solving their linear counterparts. A common approach is to do so by the Newton-Raphson method, described in textbooks on numerical analysis. This approach needs to be used when one wants to obtain an unconditionally stable method with the overall accuracy O(r2) + O(h2) and, in addition, either of the following holds: (1) D depends on the unknown u; or (2) the reaction term f s stiff in the sense described in Numerical Methods for Local Models. (The latter case may take place not for the single equation  but for two or more coupled equations of that form.)

However, when both the diffusion coefficient D does not depend on the unknown variable (although it may still depend on x and t) and the reaction term is not stiff, a simpler approach can be used. Its idea is to use the CN method (or any method with analogous accuracy and stability; e.g., ) to discretize the diffusion term, while using an explicit method to discretize the other two terms on the RHS of . Indeed, since the diffusion term is computed by the unconditionally stable CN method, then the restriction on the time step is imposed by the explicit calculation of the other two terms, which requires only that un+1 _

and is much less restrictive than  when h is small. (In fact, as we noted earlier, the condition r < const • h needs

 [3 1 _2 — 27_ Im + 7— 2

where 1/2 < 7 < 1 and, in our example, un —iun+ un nn um+1 2 um + um— 1

un um+1

2h etc. Method  is an example of a semi-implicit, or IMEX, method. It is a two-step method and, as  above, it needs to be started by computing Um1 by a single-step method (e.g., by  with 7 = 1/2 and the second parentheses on the RHS being replaced with Em). When 7 = 9/16, method  is most (among all 7's) efficient for smoothing our nonsmooth initial conditions, while when 7 = 3/4, its stability boundary is most extended along the imaginary (rA)-axis. As shown in Numerical Methods for Local Models, the latter feature of the method is essential for stability when the explicitly treated terms in  describe no or little dissipation. Mathematically, this is the case when, for example, the advection coefficient g is on the order of or greater than the reaction term f.

In connection with the last remark, let us also note that in that case, one should not use IMEX methods where the E-term is computed by any method in the Runge-Kutta family. This may sound surprising given that the latter methods were advocated in Numerical Methods for Local Models precisely because they have part of their stability region boundary extended along the imaginary (rA)-axis, which is required for successful numerical solution of ODEs with no or little dissipation. However, when a Runge-Kutta method is used to compute the E-term in an IMEX method, the stability region of such an IMEX method turns out to be close to that of the explicit Euler method, and the latter region's boundary does not extend along the imaginary (rA)-axis.

In many cases one needs to deal not with a single reaction-advection-diffusion equation but with a system of coupled such equations. However, in most applications, un m n n the diffusion terms are decoupled; that is, a term proportional to 8xUj enters only into the equation for hj but not into the equations for other unknown variables u with i—j. In such a case, one can straightforwardly generalize method  (or any other IMEX method). Indeed, moving the Im+ J-term to the LHS of the equation for the jth variable yields a tridiagonal matrix multiplying only that variable and not involving other unknown variables at the (n + 1)th time level. The RHS of the equation does involve all of the variables (because they are coupled via the E-term), but they are evaluated at the earlier time levels and therefore have been already calculated at previous steps. Thus, each of the equations can be solved time-efficiently by the Thomas algorithm. The above approach can be used for both linear and nonlinear equations of form . For linear equations only, the block Thomas algorithm (described in more advanced textbooks or in the Internet resources) can be used as an alternative. Finally, let us note that if terms DipxUj enter into the equation for uj with i—j where all of the Dj do not depend on x and u (but are allowed to depend on t), the I-part of such a system can be made uncoupled by a change ofthe unknown variables that diagonalizes matrix Dij. The corresponding methods are considered in all textbooks on 'linear algebra'.

Let us note that Matlab has a built-in solver for coupled equations of form  in one spatial dimension and with arbitrary D, g, and f. Matlab's setup allows one to also handle radially symmetric solutions in regions with circular symmetry. The user must provide the set of discretization points in the spatial interval and also code in the coefficients D, g, and f and the boundary and initial conditions. See Matlab's help for pdepe. 