High order semi-implicit multistep methods for time dependent partial differential equations

We consider the construction of semi-implicit linear multistep methods which can be applied to time dependent PDEs where the separation of scales in additive form, typically used in implicit-explicit (IMEX) methods, is not possible. As shown in Boscarino, Filbet and Russo (2016) for Runge-Kutta methods, these semi-implicit techniques give a great flexibility, and allows, in many cases, the construction of simple linearly implicit schemes with no need of iterative solvers. In this work we develop a general setting for the construction of high order semi-implicit linear multistep methods and analyze their stability properties for a prototype linear advection-diffusion equation and in the setting of strong stability preserving (SSP) methods. Our findings are demonstrated on several examples, including nonlinear reaction-diffusion and convection-diffusion problems.


Introduction
In many situations we have to deal with dynamical system arising from the spatial discretization (or more in general discretization in the phase space) of time dependent partial differential equations. For problems where the various terms have different time scales which can be easily separated the resulting system of ordinary differential equations has the form du(t) dt = f (t, u(t)) + 1 ε g(t, u(t)), with ε > 0 a small parameter emphasizing the stiffness in the system. In (1) the solution u(t) is a vector in R m which initially satisfies the condition u(0) = u 0 . Typically, the term f contains some nonlinearity or complexity that we do not want to integrate implicitly, whereas the term g/ε is stiff and requires an implicit integration. For systems of the type (1) implicit-explicit (IMEX) schemes are nowadays a very popular choice [5,6,13,27]. However, in some cases, this separation is not possible and more in general we are lead to a dynamical system in the form where the right hand side has a stiff dependence only on the last argument.
Similarly to (1) for problem (2) it is highly desirable to construct a numerical method based on evaluating implicitly only the stiff component u(t)/ε by keeping explicit the non stiff one in order to reduce the computational complexity of a fully implicit solver. Following [9] we shall call them semi-implicit methods, since they can be used in a more general context than implicit-explicit methods.
A simple first order method which realizes the above idea reads as follows u n+1 = u n + ∆tH t n+1 , u n , u n+1 ε . ( There are many circumstances in which this semi-implicit approach leads to considerable advantages, for example if the function H is linear with respect to the stiff component the numerical solution can be computed solving only linear systems of equations [8]. The extension of this simple idea to higher order, however, is not straightforward. Here, following the approach recently introduced in [9] for Runge-Kutta methods we construct high order semi-implicit discretization based on the use of linear multistep methods. To this aim, by setting v(t) = u(t)/ε in (2) we obtain the equivalent formulation as a partitioned system where v(0) = u 0 /ε. Thus, the formal equivalence among (2) and (4) allows us to adopt IMEX techniques for partitioned systems to more general cases [9,10]. We refer to [9] for a detailed discussion on the equivalence between the various forms of system usually treated with IMEX methods and to [5,13] for general references on IMEX methods based on Runge-Kutta schemes. A large literature in this direction has been devoted to the construction of IMEX Runge-Kutta schemes satisfying the asymptotic-preserving (AP) property in the case of hyperbolic problems [10,12,24,27] and for kinetic equations [11,14,15,17]. For the case of IMEX linear multistep methods, we refer to [1,2,4,6,18,20,21,30,31] for results on the construction and properties of the schemes for various types of PDEs and to [3,16] for the construction of schemes satisfying the AP property. The rest of the manuscript is organized as follows. In the next Section we introduce the IMEX multistep methods for partitioned systems and derive the corresponding semiimplicit formulation for problem (2). Next in Section 3 we detail the derivation of the schemes and analyze their stability properties for a prototype linear advection-diffusion equation and in the setting of strong stability preserving (SSP) methods. Section 4 is then devoted to present several numerical applications that confirm the validity of the present approach. The manuscript ends with some conclusions in Section 5.

Semi-implicit multistep methods
In this Section, we first introduce the general class of IMEX linear multistep schemes for partitioned systems together with some preliminary definitions. Next we recall some general results on the order conditions and, subsequently, we apply the schemes to system in the general form (2) to derive the corresponding semi-implicit formulation.

IMEX linear multistep methods for partitioned systems
Let us consider a general partitioned system in the form where y(t) ∈ R p and z(t) ∈ R q , p, q ≥ 1 and y(0) = y 0 , z(0) = z 0 . For partitioned system in the form (5) we consider schemes based on solving the first component with an explicit linear multistep method and the second with an implicit one where b −1 = 0. Implicit methods for which b j = 0, j = 0, . . . , s − 1 are referred to as backward differentiation formula (BDF). Another important class is represented by Adams methods, for whichã 0 = −1, a 0 = −1,ã j = 0, a j = 0, j = 1, . . . , s − 1.

Semi-implicit methods in predictor-corrector form
Let us now apply the previous general formulation to the case of system (2) in the partitioned form (4), i.e.
where in order to simplify notations, we removed the dependence of the parameter ε in the second argument, keeping in mind that this dependence is stiff. We obtain the semi-implicit multistep solver where initially we assume v n−j = u n−j , j = 0, . . . , s−1. Let us point out that even if the scheme doubles the number of unknown the number of evaluations of the function H(t, u(t), v(t)) is not doubled since both schemes use the same time levels. In particular, if for notation simplicity we assume the system to be autonomous and the function H to depend linearly from the second argument where the function K : R m → R m and A(u(t)) is an invertible m × m matrix, the resulting scheme can be solved without any need of an iterative solver. In fact, the second equation in (10) can be rewritten as since u n+1 is computed from the first equation in (10). For semi-implicit Runge-Kutta methods (see [9]) in the case of autonomous systems it is possible to construct the scheme in such a way that the two solutions furnished by the system at time n+1 coincide. At variance, for scheme (10) at time level n+1 we will have two distinct numerical solutions u n+1 and v n+1 approximations of the true solution u(t n+1 ) of problem (2). Note, however, that both these solutions are used by the scheme to advance in time.
In order to define a unique solution, it is natural to consider the scheme (10) as a predictorcorrector multistep method for the non stiff component, where the explicit scheme is used to predict u n+1 which is then used by the implicit solver as a corrector for v n+1 .
As an example, reverting back to the notations used initially, we can write the semi-implicit scheme for (2) asû The above scheme uniquely identifies the numerical solution at time u n+1 . In the following we will discuss the general order conditions for (6) and present different types of semi-implicit multistep methods of various order.
Remark 2.2. As a consequence of the above predictor-corrector formulation for the non stiff component, if the implicit solver has order p it is typically enough to consider an explicit solver of order p − 1. We emphasize that this predictor-corrector interpretation holds true only for the non stiff component, since the stiff one is treated implicitly by the resulting scheme.

Order conditions, stability and derivation of the schemes
In this Section we provide the details of the schemes that will be used in our numerical results. First we recall the order condition and then some general definitions concerning the stability properties. The stability properties for IMEX multistep methods are usually discussed in the case of additive systems for simple one-dimensional convection-diffusion problems [6,18,20,21]. We also refer to the recent stability analysis in [3] for the case of two dimensional partitioned systems.

Order conditions
For a partitioned system in the form (6) an order p scheme is obtained if both schemes are of order p, namely the following conditions are satisfied We recall that a s-step implicit multistep scheme can achieve order s + 1, while an s-step explicit methods has only order at most s. We refer to [6,18,20,23] for more details on the order conditions for IMEX multistep schemes in the case of additive systems. Here, we remark that the order conditions for the partitioned scheme are simpler than in the case of additive schemes where there is a coupling between the explicit and the implicit solver. In addition, in the case of system (12) only p − 1 order of accuracy is required by the explicit predictor solver to guarantee that an order p implicit corrector step yields the desired order p accuracy.

Stability properties
The stability properties are usually analyzed for simple one-dimensional linear problems of the form where λ, µ ∈ R and i is the imaginary unit. These kind of problems typically are originated by the space discretization of convection diffusion equations, hyperbolic balance laws or linear kinetic models [6,16,18,20,23]. For example, in the case of one-dimensional linear convection diffusion problems ∂u ∂t = a ∂u ∂x + D ∂ 2 u ∂x 2 , discretized by standard central differences of mesh ∆x we have (14) with where k is the frequency of the corresponding Fourier mode. In partitioned form system (14) will correspond to system Note that, since the direct application of an IMEX multistep method to a system in the additive form (14) is not equivalent to the application of the combination of an explicit and an implicit multistep scheme to the partitioned form (15) even the resulting stability analysis will differ. Applying the semi-implicit scheme (10) to (15) yields where u n−j = v n−j , j = 0, . . . , s − 1. By direct substitution of the first equation into the second we obtain the explicit form This leads to the characteristic equation Stability corresponds to the requirement that all roots have modulus less or equal one and that all multiple roots have modulus less than one.

Derivation of the schemes
In the sequel to simplify notations, we restrict to autonomous systems (this is always possible simply by augmenting the dimension of the system by one), i.e. the function H does not depend explicitly on time. Let us first point out that the simplest first order method (3) is common to multistep methods and Runge-Kutta methods and in the case of system (9) reads with v n = u n , which corresponds to a simple identity as explicit predictor for the non stiff component and backward Euler as implicit corrector for the stiff one.
Second order methods. The general form of second order schemes for (9) reads as follows with u n = v n , u n−1 = v n−1 and the solver for the non stiff component is represented by forward Euler. Popular choices for the implicit solver are obtained for α = 1/2 and β = 0 which corresponds to Crank-Nicholson, the resulting scheme will be referred to as FE-CN2, and α = 1 and β = 0 corresponding to the second order BDF scheme, we will refer to this scheme as FE-BDF2. In the case of α = 1/2 the value of β = 1/8 yields the best damping properties [6], the resulting scheme requires the additional storage of level n − 1 and is referred to as FE-MCN2.
Third order methods and higher. The most natural way to obtain third order methods is to combine as explicit solver a two-step Adams-Bashforth method with a two-step Adams-Moulton method. This same strategy actually can be used to obtain methods of higher order. We will refer to this general class of schemes as AB-AMp, where p is the order of the resulting scheme. Except for AB-AM2, which is the same as FE-CN2, these schemes in general suffer of poor stability properties when µ 0 (see Figure 2). Replacing the Adams-Moulton methods with BDF schemes with the same order yields a class of schemes with better stability properties referred to as AB-BDFp, where p is the order of the resulting scheme. In this way, AB-BDF2 is the same as FE-BDF2.
We report in Figure 1 the stability contours for various semi-implicit multistep methods up to order four.
Strong stability preserving methods. Often the time integration of PDEs requires some monotonicity properties to be satisfied. An important class of methods in this direction is represented by the so-called strong stability preserving (SSP) methods. These methods were designed specifically for solving the ODEs coming from a semi-discrete, spatial discretization of time dependent partial differential equations, especially hyperbolic PDEs and convectiondiffusion problems [19]. In summary, these schemes are stable for a certain (semi) norm under a suitable time step restriction ∆t ≤ ∆t 0 . Typically this schemes are applied to the convective part and ∆t 0 refers to the stability constraint, usually referred to as CFL condition,  Figure 1: Contours of the stability region in the case of problem (14) for some semi-implicit multistep methods of order one to four. Increasing the order of accuracy the stability constraints on ∆t become more severe for µ 0.
which links ∆t 0 = C∆t F E , where C is the CFL coefficient and ∆t F E the stability constraint of the SSP property in the Forward Euler scheme. In [19] it is shown that there are no implicit multistep SSP schemes of order higher than 1. Therefore, we can combine optimal explicit multistep methods satisfying the SSP property (17) as predictor for the non stiff component with implicit methods to improve the overall stability region in our semi-implicit schemes.
The explicit multistep SSP schemes are of the general form The optimal second order two steps explicit SSP method (C = 1/2) corresponds to the choices α = (4/5, 1/5) T and β = (8/5, −2/5) T , whereas the optimal third order explicit four steps scheme (C = 1/3) is obtained for α = (16/27, 0, 0, 11/27) T and β = (16/9, 0, 0, 4/9) T (see [19]). The corresponding schemes are denoted as SSP-AM3, SSP-BDF3 and SSP-BDF4. If one increases the number of steps, then SSP methods can be found to have larger SSP regions. Because there is no significant increase in the computational cost when the number of steps is increased, if storage is not a consideration, it may be advantageous to use a SSP multi-step methods with more steps and larger stability domain. For example, the optimal second order four steps SSP method (C = 2/3) is obtained with α = (8/9, 0, 0, 1/9) T and β = (4/3, 0, 0, 0) T . The corresponding third order semi-implicit schemes are denoted as SSP2-AM3 and SSP2-BDF3. We report in Figure 2, the stability contours of AB-AM3 and SSP-AM3. The importance of the optimal SSP property in the explicit scheme is evident and improves dramatically the stability properties of the resulting method. For large values of |µ| the semi-implicit method based on AM3 becomes stable under a reasonable CFL condition of ∆tz I ≤ 1/2. Similarly the use of SSP2 predictor increase this stability constraint to approximatively ∆tz I ≤ 3/5. We also gave in the same Figure the optimal third and fourth order methods constructed using the BDF formulae. Again the use of SSP2 predictor slightly improves the stability properties for large |µ|.

Numerical results
In what follows we investigate semi-implicit schemes for non-linear problems with stiff terms, where the stiffness can be solve efficiently by linear solvers.
We will discuss non-linear diffusion-reaction and convective-diffusion problems, following different examples extracted from [9] and [22].

Test 1: Order of convergence for reaction-diffusion system
Following the validation example of [9] we consider the non-autonomous diffusion-reaction system, where ω = (ω 1 , ω 2 ) : R + × [0, 2π) 2 → R 2 is solution of the following system the time dependent factors are α(t) = 2e t/2 and f (t) = −2e −t/2 . Accounting periodic boundary conditions, the initial data is extracted from the exact solution In order to apply the semi-implicit multi-step scheme (10) we reformulate system (19) introducing u = (u 1 , u 2 ) and v = (v 1 , v 2 ), and the operator For the spatial discretization of the diffusion operator we apply a 6th order central finite difference, on an uniform grid for the periodic domain [0, 2π) 2 with ∆x = ∆y. For u(t, x, y) evaluated on a point (t n , x i , y j ) we consider the operators  Figure 2: Contours of the stability region in the case of problem (14) for third order and fourth order SSP methods. The use of the SSP property permits to improve the stability behavior for µ 0.
(a) ω = (ω1, ω2) and similarly D yy for the y − direction, taking into account periodic boundary conditions. In order to estimate the order of accuracy in time we proceed refining the space step and time step with stability conditions ∆t = λ∆x, where we choose λ = 0.5. Thus we monitor the 1 error decay of the numerical solution ω ij (t) at final time T = 2 is the numerical solution evaluated on N x = 2 k points, with k = 6, 7, 8 and ω(x i , y j , t n ) is the exact solution evaluated at (x i , y j , T ).
In Figure 3 we report the 1 error decay for several schemes. We compare the accuracy of the solution (a) ω = (ω 1 , ω 2 ), on the left side, with respect to (b) the accuracy of ω 2 , observing better convergence properties for the second component which non-autonomous and independent of ω 1 . Nonetheless in both cases 5th order of convergence is observed.
Initial data and parameters are selected in order to match the test proposed in [22] and inspired from [26]. In order to apply the semi-implicit multi-step scheme (10) we reformulate system (21) introducing u = (u 1 , u 2 ) and v = (v 1 , v 2 ), and the operator where now we treat explicitly the diffusion terms, since D 1 , D 2 are negligible with respect to the reaction coefficients. Thus the linear reaction terms are taken implicitly, whereas the non-linear term is taken implicit only in the ω 1 component. We employ LM scheme SSP3-BDF4 for time integration and forth order central difference for the space discretization, introducing the operators and similarly D yy for the y − direction, taking into account periodic boundary conditions. Figure 4 reports the evolution of the component ω 2 (t, x, y) for different times. Starting from a symmetric concentration Gray-Scott model produces spot multiplication, resembling cell division process. The computational domain is discretized with N x = N y = 200 points in space, final time is T = 1500 with uniform step ∆t = ∆x/2.

Test 3: Nonlinear convection-diffusion
We finally consider the nonlinear convection diffusion model defined on the full plan x ∈ R 2 as follows ∂ t ω + (E + µ∇ log(ω)) · ∇ω = µ∆ω ω 0 (0, x) = e − x 2 /2 (24) where E = (1, 1) and µ = 0.5. The initial data is extracted from the exact solution given by In order to solve numerically (24) we introduce the operator H(t, u, v) = − (E + µ∇ log(u)) · ∇v + µ∆v where we treat the convection and diffusion terms implicitly, on the computational domain [−10, 10] up to final time T = 1. We choose uniform space step ∆x = ∆y = 0.1 and ∆t = ∆x/2 For time integration we employ SSP3-BDF4, and 4th order central difference for diffusion operator, as in (23), and for the convective term we account the operator and equivalently D y for the y−direction.
In Figure 5 we report the numerical solution at different times from t = 0 to t = 1.

Conclusions
We derived high order semi-implicit schemes based on linear multistep methods. The schemes have been constructed following the approach recently introduced in [9] for Runge-Kutta methods. The resulting time discretizations have a predictor-corrector structure and, compared with Runge-Kutta methods, do not require additional order conditions so that they can easily reach high order accuracy. Numerical tests for schemes up to fifth order accuracy have been presented in the case of nonlinear reaction-diffusion and convection-diffusion problems.