 # 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.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 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)), (1)

with a small parameter emphasizing the stiffness in the system. In (1) the solution

is a vector in

which initially satisfies the condition . Typically, the term contains some nonlinearity or complexity that we do not want to integrate implicitly, whereas the term 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

 du(t)dt=H(t,u(t),u(t)ε) (2)

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 by keeping explicit the non stiff one in order to reduce the computational complexity of a fully implicit solver. Following  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

 un+1=un+ΔtH(tn+1,un,un+1ε). (3)

There are many circumstances in which this semi-implicit approach leads to considerable advantages, for example if the function is linear with respect to the stiff component the numerical solution can be computed solving only linear systems of equations.

The extension of this simple idea to higher order, however, is not straightforward. Here, following the approach recently introduced in  for Runge-Kutta methods we construct high order semi-implicit discretization based on the use of linear multistep methods. To this aim, by setting in (2) we obtain the equivalent formulation as a partitioned system

 du(t)dt = H(t,u(t),v(t)) (4) εdv(t)dt = H(t,u(t),v(t)),

where . 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  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, 17, 15, 14]. 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 semi-implicit 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.

## 2 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.

### 2.1 IMEX linear multistep methods for partitioned systems

Let us consider a general partitioned system in the form

 dy(t)dt = F(t,y(t),z(t)) (5) dz(t)dt = G(t,y(t),z(t)).

where and , and , .

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

 yn+1=−s−1∑j=0~ajyn−j+Δts−1∑j=0~bjF(tn−j,yn−j,zn−j) (6) zn+1=−s−1∑j=0ajzn−j+Δts−1∑j=−1bjG(tn−j,yn−j,zn−j)

where . Implicit methods for which , are referred to as backward differentiation formula (BDF). Another important class is represented by Adams methods, for which , , , , .

###### Remark 2.1.

Classical systems in additive form (1), i.e.

 du(t)dt=f(t,u(t))+1εg(t,u(t)) (7)

with initial data can be also written in partitioned form by defining , , and rewriting

 dy(t)dt = F(t,y(t),z(t)) (8) dz(t)dt = G(t,y(t),z(t)),

for any initial data such that .

### 2.2 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.

 du(t)dt = H(t,u(t),v(t)) (9) dv(t)dt = H(t,u(t),v(t)),

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

 un+1=−s−1∑j=0~ajun−j+Δts−1∑j=0~bjH(tn−j,un−j,vn−j) (10) vn+1=−s−1∑j=0ajvn−j+Δts−1∑j=−1bjH(tn−j,un−j,vn−j)

where initially we assume , . Let us point out that even if the scheme doubles the number of unknown the number of evaluations of the function 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 to depend linearly from the second argument

 H(u(t),v(t))=K(u(t))+A(u(t))v(t), (11)

where the function and is an invertible 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

 vn+1= − s−1∑j=0ajvn−j+Δts−1∑j=0bjH(un−j,vn−j) + Δtb−1(K(un+1)+A(un+1)vn+1).

or equivalently in explicit form

 vn+1= (I−b−1ΔtA(un+1))−1(−s−1∑j=0ajvn−j+Δts−1∑j=0bjH(un−j,vn−j) +Δtb−1(K(un+1))),

since is computed from the first equation in (10).

For semi-implicit Runge-Kutta methods (see ) 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

coincide. At variance, for scheme (

10) at time level we will have two distinct numerical solutions and approximations of the true solution 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 predictor–corrector multistep method for the non stiff component, where the explicit scheme is used to predict which is then used by the implicit solver as a corrector for .

As an example, reverting back to the notations used initially, we can write the semi-implicit scheme for (2) as

 ^un+1 = −s−1∑j=0~ajun−j+Δts−1∑j=0~bjH(tn−j,un−j,un−jε) (12) un+1 = −s−1∑j=0ajun−j+Δts−1∑j=0bjH(tn−j,un−j,un−jε) +Δtb−1H(tn+1,^un+1,un+1ε).

The above scheme uniquely identifies the numerical solution at time . 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 it is typically enough to consider an explicit solver of order . 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.

## 3 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  for the case of two dimensional partitioned systems.

### 3.1 Order conditions

For a partitioned system in the form (6) an order scheme is obtained if both schemes are of order , namely the following conditions are satisfied

 1+s−1∑j=0~aj=0, 1+s−1∑j=0aj=0, (13) 1−s−1∑j=1j~aj=s−1∑j=0~bj, 1−s−1∑j=1jaj=s−1∑j=−1bj, 12+s−1∑j=1j22~aj=−s−1∑j=1j~bj, 12+s−1∑j=1j22aj=b−1−s−1∑j=1jbj ⋮ ⋮ 1p!+s−1∑j=1(−j)pp!~aj 1p!+s−1∑j=1(−j)pp!aj =s−1∑j=1(−j)p−1(p−1)!~bj, =b−1(p−1)!+s−1∑j=1(−j)p−1(p−1)!bj.

We recall that a -step implicit multistep scheme can achieve order , while an -step explicit methods has only order at most . We refer to [6, 20, 18, 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 order of accuracy is required by the explicit predictor solver to guarantee that an order implicit corrector step yields the desired order accuracy.

### 3.2 Stability properties

The stability properties are usually analyzed for simple one-dimensional linear problems of the form

 du(t)dt=iλu(t)+μu(t), (14)

where and 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, 20, 18, 23]. For example, in the case of one-dimensional linear convection diffusion problems

 ∂u∂t=a∂u∂x+D∂2u∂x2,

discretized by standard central differences of mesh we have (14) with

 λ=aΔxsin(kΔx),μ=2D(Δx)2(cos(kΔx)−1),

where is the frequency of the corresponding Fourier mode.

In partitioned form system (14) will correspond to system

 du(t)dt = iλu(t)+μv(t) (15) dv(t)dt = iλu(t)+μv(t).

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

 un+1= −s−1∑j=0~ajun−j+Δts−1∑j=0~bj(iλun−j+μvn−j) (16) vn+1= (11−μb−1Δt)(−s−1∑j=0ajvn−j+Δts−1∑j=0bj(iλun−j+μvn−j) + Δtb−1iλun+1),

where , . By direct substitution of the first equation into the second we obtain the explicit form

 vn+1= (11−μb−1Δt)(−s−1∑j=0ajvn−j+Δts−1∑j=0bj(iλvn−j+μvn−j) + Δtb−1iλ(−s−1∑j=0~ajvn−j+Δts−1∑j=0~bj(iλvn−j+μvn−j))).

This leads to the characteristic equation

 ζs(1−b−1zR)+ρ(ζ)−(zR+zI)σ(ζ)+b−1zI(~ρ(ζ)−(zR+zI)~σ(ζ))=0,

where , and

 ρ(ζ) = s−1∑j=0ajζs−j−1,σ(ζ)=s−1∑j=0bjζs−j−1, ~ρ(ζ) = s−1∑j=0~ajζs−j−1,~σ(ζ)=s−1∑j=0bjζs−j−1.

Stability corresponds to the requirement that all roots have modulus less or equal one and that all multiple roots have modulus less than one.

### 3.3 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 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

 un+1 = un vn+1 = vn+ΔtH(un+1,vn+1),

with , 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

 un+1 = un+ΔtH(un,vn) vn+1 = 12α+1(4αun−(2α−1)un−1)+Δt2α+1((2α+β)H(un+1,vn+1) +2(1−α−β)H(un,vn)+βH(un−1,vn−1)),

with , and the solver for the non stiff component is represented by forward Euler. Popular choices for the implicit solver are obtained for and which corresponds to Crank-Nicholson, the resulting scheme will be referred to as FE-CN2, and and corresponding to the second order BDF scheme, we will refer to this scheme as FE-BDF2. In the case of the value of yields the best damping properties , the resulting scheme requires the additional storage of level 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-AM, where 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 (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-BDF, where 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. 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.

#### 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 convection-diffusion problems . In summary, these schemes are stable for a certain (semi) norm

 ∥un+1∥≤∥un∥ (17)

under a suitable time step restriction . Typically this schemes are applied to the convective part and refers to the stability constraint, usually referred to as CFL condition, which links , where is the CFL coefficient and the stability constraint of the SSP property in the Forward Euler scheme.

In  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

 un+1=s−1∑j=0(αjun−j+ΔtβjH(un−j,vn−j)),αj≥0. (18)

The optimal second order two steps explicit SSP method () corresponds to the choices and , whereas the optimal third order explicit four steps scheme () is obtained for and (see ). 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 () is obtained with and . 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-AM and SSP-AM. 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 . Similarly the use of SSP2 predictor increase this stability constraint to approximatively . 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 . 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.

## 4 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  and .

### 4.1 Test 1: Order of convergence for reaction-diffusion system

Following the validation example of  we consider the non-autonomous diffusion-reaction system, where is solution of the following system

 ∂tω1 =Δω1−α(t)ω21+92ω1+ω2+f(t) (19) ∂tω2 =Δω2+72ω2,

the time dependent factors are and . Accounting periodic boundary conditions, the initial data is extracted from the exact solution

 ω1(t,x,y)=e−t/2(1+cos(x))ω2(t,x,y)=e−t/2cos(2x)(x,y)∈[0,2π)2 (20)

In order to apply the semi-implicit multi-step scheme (10) we reformulate system (19) introducing and , and the operator

 H(t,u,v)=⎛⎜ ⎜⎝Δv1−α(t)u1v1+92v1+v2+f(t)Δv2+52v2⎞⎟ ⎟⎠.

For the spatial discretization of the diffusion operator we apply a 6th order central finite difference, on an uniform grid for the periodic domain with . For evaluated on a point we consider the operators

 (Dxxun)ij=2uni+3j−27uni+2j+270uni+1j−240unij+270uni−1j−27uni−2j+2uni−3j180Δx2.

and similarly for the , 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 . Thus we monitor the error decay of the numerical solution at final time

 ℓ1(ω(k))=ΔxΔy∑i,j∥ω(k)ij(T)−ω(xi,yj,T)∥,

where is the numerical solution evaluated on points, with and is the exact solution evaluated at .

In Figure 3 we report the error decay for several schemes. We compare the accuracy of the solution (a) , on the left side, with respect to (b) the accuracy of , observing better convergence properties for the second component which non-autonomous and independent of . Nonetheless in both cases 5th order of convergence is observed.

### 4.2 Test 2: Non-linear reaction diffusion system

We consider the Gray-Scott model, as studied in [26, 32, 22],

 ∂tω1=σ1Δω1−ω1ω22+γ(1−ω1) (21) ∂tω2=σ2Δω2+ω1ω22−(γ+κ)ω2

where with periodic boundary conditions, and initial data

 ω1(0,x,y)=1−2ω2(0,x,y)ω2(0,x,y)=14sin2(4πx)sin2(4πy)(x,y)∈[−1/4,1/4]2 (22)

and zero otherwise on the squared domain . The diffusion and reaction parameters are

 σ1=8×10−5,σ2=4×10−5,γ=0.024,κ=0.06.

Initial data and parameters are selected in order to match the test proposed in  and inspired from . In order to apply the semi-implicit multi-step scheme (10) we reformulate system (21) introducing and , and the operator

 H(t,u,v)=(D1Δu1−v1u22+γ(1−v1)D2Δu2+v1u22−(γ+κ)v2),

where now we treat explicitly the diffusion terms, since 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 component. We employ LM scheme SSP3-BDF4 for time integration and forth order central difference for the space discretization, introducing the operators

 (Dxxun)ij =−uni+2j+16uni+1j−30unij+16uni−1j−uni−2j12Δx2. (23)

and similarly for the , taking into account periodic boundary conditions.

Figure 4 reports the evolution of the component for different times. Starting from a symmetric concentration Gray-Scott model produces spot multiplication, resembling cell division process. The computational domain is discretized with points in space, final time is with uniform step . Figure 4: Test 2: (Gray-Scott model). Solution of the reaction-diffusion system (21) in the component ω2 at different time frames obtained with scheme SSP-BDF4 on the square Ω=[−1,1]2 with uniform space-time discretization with Δx=0.01 and Δt=Δx, and periodic boundary conditions. Parameters of the model are σ1=8×10−5, σ2=4×10−5, γ=0.024 and κ=0.06.

### 4.3 Test 3: Nonlinear convection-diffusion

We finally consider the nonlinear convection diffusion model defined on the full plan as follows

 ∂tω+(E+μ∇log(ω))⋅∇ω=μΔω (24) ω0(0,x)=e−∥x∥2/2

where and . The initial data is extracted from the exact solution given by

 ω(t,x)=1√4μt+1e−∥x−Et∥28μt+2.

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 up to final time . We choose uniform space step and 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

 (Dxun)ij =−ui+2j+8ui+1j−8ui−1j+ui−2j12Δx,

and equivalently for the direction.

In Figure 5 we report the numerical solution at different times from to . Figure 5: Test 3: (Non-linear convection diffusion). Numerical solution of model (24) on the domain [−10,10] up to final time T=1. Space step Δx=Δy=0.1 and Δt=Δx/2. Time integration is performed with LM scheme SSP3-BDF4 and 4th order central difference for convective and diffusion operator.

## 5 Conclusions

We derived high order semi-implicit schemes based on linear multistep methods. The schemes have been constructed following the approach recently introduced in  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.

## Acknowledgments

This research is partially supported by PRIN (2019-2021): "Innovative numerical methods for evolutionary partial differential equations and applications" and by INdAM-GNCS (2019) grant: "Approssimazione numerica di problemi di natura iperbolica ed applicazioni".

## References

•  G. Akrivis, Implicit–explicit multistep methods for nonlinear parabolic equations, Math. Comput., 82 (2012) 45–68.
•  G. Akrivis, M. Crouzeix and C. Makridakis, Implicit–explicit multistep methods for quasilinear parabolic equations, Numer. Math., 82 (1999) 521–541.
•  G. Albi, G. Dimarco, L. Pareschi, Implicit-Explicit multistep methods for hyperbolic systems with multiscale relaxation, arXiv:1904.03865, 2019.
•  G. Albi, M. Herty, L. Pareschi, Linear multistep methods for optimal control problems and applications to hyperbolic relaxation systems, Appl. Math. Comp. 34, 2019, 460–477.
•  U. M. Ascher, S. J. Ruuth, R. J. Spiteri, Implicit-Explicit Runge-Kutta Methods for Time-Dependent Partial Differential Equations, Appl. Numer. Math, 1997, V. 25, pages 151–167.
•  U. M. Ascher, S. J. Ruuth, B. T. R. Wetton, Implicit-Explicit methods for time-dependent partial differential equations, SIAM J. Num. Anal. 32, (1995), 797–823.
•  A.L. Bertozzi, The mathematics of moving contact lines in thin liquid films. Not. AMS 45, 689–697, (1998).
•  S. Boscarino, R. Bürger, P. Mulet, G. Russo, L. M. Villada, On linearly implicit IMEX Runge-Kutta methods for degenerate convection-diffusion problems modeling polydisperse sedimentation, Bull. Braz. Math. Soc., New Series 47(1), (2016), 171–185.
•  S. Boscarino, F. Filbet, G. Russo, High Order Semi-implicit Schemes for Time Dependent Partial Differential Equations, J. Sci. Comp., Vol. 68, (2016), 975–1001.
•  S. Boscarino, L. Pareschi, On the asymptotic properties of IMEX Runge–Kutta schemes for hyperbolic balance laws, J. Comput. Appl. Math., 316 (2017), pp. 60-73.
•  S. Boscarino, L. Pareschi, G. Russo, Implicit-Explicit Runge-Kutta schemes for hyperbolic systems and kinetic equations in the diffusion limit, SIAM J. Sci. Comput., 35 (2011), pp. A22-A51.
•  S. Boscarino, L. Pareschi, G. Russo, A unified IMEX Runge-Kutta approach for hyperbolic systems with multiscale relaxation, SIAM J. Num. Anal., 55, No. 4, (2017), pp. 2085-2109.
•  M. H. Carpenter, C. A. Kennedy, Additive Runge-Kutta schemes for convection-diffusion-reaction equations Appl. Numer. Math. 44 (2003), no. 1-2, 139–181.
•  G. Dimarco, L. Pareschi, Numerical methods for kinetic equations, Acta Numerica 23, 369–520, 2014.
•  G. Dimarco, L. Pareschi, Asymptotic-preserving IMEX Runge-Kutta methods for nonlinear kinetic equations, SIAM J. Num. Anal. (2013), 1064–1087.
•  G. Dimarco, L. Pareschi, Implicit-Explicit linear multistep methods for stiff kinetic equations. SIAM J. Num. Anal. 55(2), (2017), 664–690.
•  F. Filbet, S. Jin, A class of asymptotic preserving schemes for kinetic equations and related problems with stiff sources, J. Comp. Phys., 229 (2010), pp. 7625-7648.
•  J. Frank, W. Hundsdorfer, J. G. Verwer, On the stability of implicit-explicit linear multistep methods. Appl. Num. Math. 25, (1997), 193–205.
•  S. Gottlieb, C-W. Shu, E. Tadmor, Strong Stability-Preserving High-Order Time Discretization Methods, SIAM Review 43, (2001), 89–112.
•  W. Hundsdorfer, S. J. Ruuth, IMEX extensions of linear multistep methods with general monotonicity and boundedness properties, J. Comp. Phys. 225, (2007), 2016–2042.
•  W. Hundsdorfer, S. J. Ruuth, R. J. Spiteri, Monotonicity-preserving linear multistep methods, SIAM J. Numer. Anal. 41, (2003) 605–623.
•  W. Hundsdorfer, J. Verwer, Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations, Springer Series in Computational Mathematics, (2003).
•  E. Hairer, G. Wanner, Solving Ordinary Differential Equation II: stiff and Differential Algebraic Problems. Springer Series in Comput. Mathematics, Vol. 14, Springer-Verlag 1991, Second revised edition 1996.
•  S. Jin, Efficient Asymptotic-Preserving (AP) Schemes For Some Multiscale Kinetic Equations, SIAM J. Sci. Comput., 21(2), (1999), 441–454.
•  G. Naldi, L. Pareschi, G. Toscani , Relaxation schemes for PDEs and applications to degenerate diffusion problems, Surveys on Mathematics for Industry (now European Journal for Applied Mathematics), 10, 315–343, (2002).
•  J. E. Pearson, Complex patterns in a simple system, Science., 261(5118), (1993) 189-192.
•  L. Pareschi, G. Russo, Implicit-Explicit Runge-Kutta schemes and applications to hyperbolic systems with relaxations, J. Sci. Comput., 25 (2005), 129–155
•  L. Pareschi, G. Russo, G. Toscani, A kinetic approximation to Hele–Shaw flow, C. R. Math., Acad. Sci. Paris, Ser. I 338(2), 178–182 (2004).
•  L. Pareschi, M. Zanella, Structure preserving schemes for nonlinear Fokker-Planck equations and applications, J. Sci. Comput. 74 (2018), no. 3, 1575–1600.
•  R. R. Rosales, B. Seibold, D. Shirokoff, D. Zhou, Unconditional stability for multistep IMEX schemes: Theory, SIAM J. Numer. Anal. 55(5), (2017), 2336–2360.
•  B. Seibold, D. Shirokoff, D. Zhou, Unconditional stability for multistep IMEX schemes: Practice, J. Comp. Phys. 376, (2019), 295–321.
•  K. Zhang, J.C.F. Wong, R. Zhang, Second-order implicit–explicit scheme for the Gray–Scott model, J. Comput. Appl. Math. 213, (2008), 559–581.