1 Introduction
Stiff ODEs are present in many applications, notably when the method of lines is used for the semidiscretization of PDEs. Both the parabolic stiffness due to diffusion and the hyperbolic stiffness due to the CFL condition are common sources of stiffness in PDEs. Stiffness imposes small timestep sizes in order to avoid numerical instabilities. In general, implicit timestepping allows a larger timestep size but it requires the solution of an implicit equation at every step/stage of the method. IMEX methods decompose the righthandside of the ODE as the sum of two functions
(1) 
where is the explicit part of the decomposition and is the implicit part. An efficient decomposition should place nonstiff terms in and stiff terms in while keeping the implicit equation simple to solve. For IMEX, there is the restriction that the implicit equation needs to be solved up to a given precision in order to avoid introducing errors that could overwhelm the local truncation error; i.e., there is a smallenoughresidual restriction to be fulfilled. The objective of the present article is to remove this restriction.
A new method called shortcutIMEX (SIMEX) is introduced. It allows an iterative solver applied to the implicit equations to be interrupted at a given iteration. The reason SIMEX does not introduce errors is the use of a residual balanced decomposition: after the iterative solver is interrupted, any remaining residual of the implicit equation is accounted for in the explicit timestep by a suitable redefinition of the implicitexplicit decomposition. This withinstep adjustment of the decomposition requires a modification of traditional algorithms used with IMEX. Many IMEX timeintegrators can be adapted to SIMEX; namely, singlydiagonallyimplicit RungeKutta schemes (SDIRK), explicit first stage SDIRK schemes (ESDIRK)[2, 26, 16, 22, 5, 29, 17, 27, 14, 1], and multistep schemes [3, 23, 13, 10, 6, 28, 30, 4]. Here, we refer to IMEXRK and IMEXMS to distinguish RungeKutta (RK) methods from multistep (MS) methods.
Recently, IMEX methods have been intensely researched in many fields of science. The choice of stiff terms being placed in may vary. For example, the stiff term can be a parabolic term [20, 6, 27], or it can be a term connected to acoustic waves in the atmosphere [23, 10, 29, 14, 4, 31], or it can be a term associated to the CFL condition on a refined grid [16, 22]. In fluid mechanics, when there is parabolic stiffness, IMEXRK allows highorder timestepping to be used [16]. Along with discontinuous Galerkin, IMEX has been used in numerous applications [16, 27, 25, 11]. The stiffness of the CFL condition has been addressed by applying IMEX methods to the relaxation of hyperbolic equations [5]. IMEX has been compared to the deferred correction method [24] in atmospheric flows. In simulations of thermal convection in the Earth’s liquid outer core, IMEX [13, 21] has been compared to exponential integrators [12] where IMEX is found to have a computational advantage.
Originally, the IMEX method was proposed in [9] as a loworder multistep method; short after, it appeared in [8, 7] associated to additive RungeKutta (ARK) methods. An IMEXRK method can also be called an ARK method; more precisely, it can be called an method, where an method partitions the ODE in parts. A comprehensive review of ARK methods can be found in [20, 18]. Both IMEXRK and IMEXMS methods can be blended into the generallinearmethod (IMEXGLM) by combining multiple stages and steps [31].
In some applications, the implicit part is linearized in order to avoid solving a nonlinear system at each step [22, 14]. For example, one may rewrite the ODE as
(2) 
and place only the linear term in the implicit part of the IMEX decomposition. can be an approximation to the Jacobian derivative of with respect to . In IMEXRK, the choice of the decomposition can be performed at the beginning of each step [14] or remain fixed throughout timeintegration [22].
In this article we explore how SIMEX overcomes parabolic stiffness using RK methods. In examples, we apply it to scalar advectiondiffusionreaction PDEs. The paper is organized as follows. A summary of IMEXRK methods, its notation and its tableau, is found in Section 2. The SIMEXRK method is introduced in Section 3 along with the algorithm; some theoretical considerations are also given in this section. Examples where SIMEXRK is applied to ODEs are given in Section 4, these ODE examples arise from a coarse spatial discretization of advectiondiffusionreaction PDEs. The stability regions of the new method is discussed in Section 5. An example of application of the new method for an advectiondiffusion PDE is given in Section 6 where its numerical convergence with respect to spacetime grid refinement is studied. Conclusions are drawn in Section 7.
2 IMEX RungeKutta methods
When considering IMEX, there are two large classes of timestep methods, RK methods and MS methods. In this section, the IMEXRK methods with ESDIRK schemes are summarized.
A number of IMEXRK schemes are available to choose from, or to be tailored to one’s need. The choice of scheme may be guided, for example, by its classical order, or by the stepsizecontrol with embedded RKpairs, or by its economy in computational storage [15], or if it has denseoutput [19]. There are schemes with all these properties within the class of ESDIRK schemes. They are commonly found in IMEXRK methods [19, 18, 20] connected to PDE applications [2, 26, 16, 22, 5, 29, 17, 27, 14, 1]. Because the implicit equation keeps the same structure at every implicit stage, ESDIRK schemes save some computational effort by reusing Jacobians and matrix factorizations when solving implicit equations at each successive RKstage.
An IMEXRK method with an ESDIRK scheme consists of two specially crafted RKmethods where stages are computed in alternation, each stage “acting” exclusively either on or . The stages are then combined at the end of each step to integrate the full ODE. Below we write the general form of a joint Butcher’s tableau for an ESDIRK method with stages:
(3) 
The empty entries are zero. The explicit method is on the righthand side. The coefficients and are different on each side of the tableau but the values of and are the same. For ESDIRK, the same element repeats along the diagonal, for , except for which is zero. Details about the order condition for an IMEXRK tableau can be found in [18] and [20].
With the above tableau, the IMEXRK algorithm is given in Algorithm 1. It describes the algorithm for a single step of size where the input is the current approximation of denoted by . The algorithm calls a solver procedure denoted by . The solver must return an approximate solution of the implicit equation,
(4) 
which must be solved for .
Observe that only the righthand side of equation Eq. 4 changes from stage to stage, this is a property of both ESDIRK and SDIRK schemes. A comment about line 5 of this algorithm: the evaluation of at this line is unnecessary whenever is evaluated at line 4 with the same arguments; this occurs when the residual of Eq. 4 is evaluated inside the function .
An example of a simple tableau that can be used with the above algorithm is the following:
(5) 
This is a second order stable scheme which is a combination of CrankNicolson and Heun’s methods henceforth called the CNH method. For CNH, Algorithm 1 is equivalent to the following formulas: given , solve for the implicit equation
(6) 
and compute
(7) 
The three IMEXRK schemes being used in this article are: (i) CNH given above, (ii) ARK548 (reference [18], page 49, labeled as ARK5(4)8L[2]SA), and (iii) ARK436 (reference [18], pages 47 and 48, labeled as ARK4(3)6L[2]SA). These two ARK schemes have embedded pairs of order and ; the number of RKstages is . Both ARK schemes are stiffly accurate with stage order 2 where the implicit tableau is stable.
Motivated by the structure of Eq. (4) where both and are for , this equation can be rewritten as
(8) 
where , , and the remaining terms are accounted for in the righthand side . The important remark is that, in an IMEX method, the solution of Eq. Eq. 8 needs to be at least as accurate as the timestepping method. This restriction is removed in the next section with the introduction of the SIMEX method.
3 SIMEX for RungeKutta methods
In this section we define the shortcutIMEX method for RK methods where ESDIRK is used in the implicit timestepping tableau.
A SIMEXRK method is defined with three elements: a tableau like (3), an implicitexplicit decomposition as in Eq. (1), and an implicitstepfilter which we refer simply as a filter. A filter can be defined, for example, from an iterative solver by stopping the iterative solver after a fixed number of iterations. Here, the distinction between a solver and a filter is necessary because neither precision nor convergence are required from a filter. Indeed, a filter may yield a low precision approximation of Eq. Eq. 8. This is not a problem for SIMEX. Also, the filter can be changed from step to step but not within a step, i.e., the same filter must be used in all RKstages belonging to a single step. We use the notation when it is necessary to emphasize that the filter changes between steps.
The input arguments of a filter are the same elements that define Eq. Eq. 8
: the vectors
and , the real number , and the function . The output of is a vector with the same dimensions as . As a mathematical function, we use the notation and we simplify the notation to when some arguments are constant. Notice that depends on time according to the dependence on only.In SIMEX, timestep precision is obtained by the residual balanced decomposition which we explain next. Henceforth we call the decomposition , , of Eq. (1) as the prototype decomposition, or protodecomposition. We introduce a new decomposition
(9) 
which we call the residual balanced decomposition (RBD). RBD is formally defined by choosing
(10) 
where is the inverse function of the filter with respect to the argument . Also, as in Eq. (8), . Consequently, the definition of must be
(11) 
Fortunately, there is an easy way to evaluate and without the need of computing . Yet, the existence of is important for the theoretical justification of the method. A discussion about the existence of is given after the SIMEXRK algorithm.
In order to obtain a simple way to evaluate and , we observe that
is an exact solution of the implicit equation for the RDB decomposition:
(12) 
Thus, the actual value of can be computed from Eq. (12) by the formula
(13) 
where . Once is computed using Eq. (13), is computed using Eq. (11) with . Thus, RBD provides both an exactly solvable implicitstep equation, Eq. (12), and a simple way to evaluate by using Eq. (13).
RBD has an interpretation relating it to the original protodecomposition. The implicit equation for the protodecomposition, Eq. Eq. 8, can be rewritten as
(14) 
and is an approximate solution where the residual is
Thus, the interpretation of the RBD is that it balances this residual of Eq. Eq. 8 by transferring it to the explicit part of the RDB decomposition, namely to Eq. (11).
Algorithm 2 shows how RBD can be implemented, this is called the SIMEXRK algorithm; it describes a single step of size where the current value is an input. The decomposition and is computed in lines 7 and 8 of the algorithm; they follow from Eqs. (13) and (11) respectively.
A few remarks about the filter : Observe that should be true for all filters, it is unreasonable to assume otherwise. Also, observe that if no information about Eq. Eq. 8 is used, then the choice of filter should be the identity ; henceforth this is called the defaultfilter. Of course, the defaultfilter does not help to improve numerical stability, it guarantees the “best/easiest” inversion . In contrast, a filter that exactly solves Eq. Eq. 8 is called an exactfilter. When there is a square matrix such that , the filter is called linear; otherwise it is called nonlinear. A linear filter can be used when Eq. Eq. 8 is nonlinear also. For example, the filter can be defined as a single Newtonlike iteration applied to Eq. (8). When is linear, RBD recasts the protodecomposition of Eq. Eq. 1 into the decomposition of Eq. Eq. 2 where and are given by Eq. Eq. 10.
Even though is not directly evaluated in Algorithm 2, a proof of its existence implies that RBD is a valid IMEX decomposition for Eq. (1). Such proof places Algorithm 2 in the same theoretical grounds as Algorithm 1
. Some heuristic aspects of the existence and smoothness of
are discussed in the remainder of this section. We argue that these properties are easily obtainable when is derived from the usual range of iterative solvers using a small number of iterations. In this discussion, assume does not depend on ; thus can be removed from the argument of all functions. Denote the function on the lefthand side of Eq. (8) by , , and denote its inverse by . Also, denote the Jacobian by , this leads to and .First, consider the case where is a linear filter. If is an approximate solution of Eq. (8) for any with small norm, then one should expect the matrix norm to be small. This expresses that is close to the nonsingular matrix , thus should be nonsingular too. Even if is merely a rough approximation to the solution of Eq. (8), it is still reasonable to expect that should have an inverse. This is specially clear when is small because then . This heuristic argument can be quantified in terms of matrix norms; for example, assume that is close enough to so that
(15) 
Then is nonsingular because, given an arbitrary unitary vector , it follows that
In the limit as approaches zero, inequality Eq. 15 reduces to which is particularly easy to satisfy.
Second, when is a nonlinear filter, the discussion about the existence and smoothness of becomes more complex. For instance, observe that even if is very close to this does not guarantee that is close to ; i.e., one cannot directly conclude that is close to the nonsingular matrix . Thus, unlike linear filters, the existence of may depend on bounds of the second derivative . Nevertheless, if one has a reason to believe that is sufficiently close to a linear filter, then the existence of
should be expected without apprehension. Another aspect of nonlinear filters is the smoothness of
. It needs to have at least as many continuous derivatives as the derivatives necessary to bound the local truncation error of the RKmethod being used. From the inverse function theorem, it is known that has as many continuous derivatives as . When is computed by subsequent iterations of Newton’s method, the Jacobian of one iteration uses the answer of a previous iteration as argument; thus, the order of differentiation with respect to increases at every iteration. Consequently, the regularity of and could have less continuous derivatives than . Here, we explore in Section 4 an example of nonlinear filter where things work nicely for SIMEXRK even though we do not verify if is either invertible or smooth.The last remark about Algorithm 2 is that the filter can be chosen from a sequence of filters, , , , , according to a stabilization criterion. This choice is made at the first execution of line number 6 in Algorithm 2, i.e., when . Once chosen, the same filter must be used for . This is summarized in Algorithm 3 which can replace line 6 in Algorithm 2:
The sequence of filters can be defined by the successive iterations of an iterative solver. Here, we do not explore Algorithm 3 in numerical experiments; instead, we experiment with each filter individually in order to understand its properties. Nevertheless, it is clear that Algorithm 2 produces almost the same calculations as Algorithm 1 whenever the stabilization criterion is chosen equal to the stopping criterion used within the iterative solver .
4 SIMEX for ODEs
Timestep convergence is studied in this section using examples of ODEs that originate from the semidiscretization of PDEs. The purpose of these examples is to show that SIMEXRK can maintain accuracy even if the computational effort placed in the filter is small compared to the effort placed in the solver of IMEXRK. In this section, the numerical convergence rate is studied as is decreased while keeping the spatial discretization parameter fixed.
For the first example, we consider a 1D forced heatequation
where has zero Dirichlet boundary conditions in and is the forcing term such that the exact solution is
(16) 
The equation is solved by the method of lines where the spatial derivative is discretized by secondorder finite differences at the points which are apart. The resulting ODE is
(17) 
where is the differentiation matrix for the second derivative and the forcing term is a timedependent vector . The protodecomposition is chosen with and .
Let denote the matrix of the linear system found in Eq. Eq. 8, this matrix is tridiagonal. Denote the diagonal of by and denote the remaining part by . The filters used in SIMEXRK are computed from Jacobi iterations. The filter is defined as the th element of the sequence where .
In the numerical experiment shown in Fig. 1, the solutions of Eq. Eq. 17 are computed using Algorithm 2, each solution with a different filter , . In this example, Eq. Eq. 17
has 9 degrees of freedom and
; this is a lowdimensional ODE which is not particularly stiff. The tableau used is the fifth order ARK548 scheme. The reference exact solution for Eq. Eq. 17 is obtained from the “ode45” routine in Matlab with absolute and relative tolerances set to .The numerical convergence rate for the SIMEXRK method is fifth order for all filters; this is shown in Fig. 1(a). Regardless the number of iterations, the SIMEXRK method maintains the correct convergence rate. The numerical convergence of IMEXRK in Algorithm 1 was carried out by replacing the solver (line 4 of the algorithm) with each one of the filters ; the result is shown in Fig. 1(b). IMEXRK needs at least three Jacobi iterations in order to maintain fifth order convergence.
The second example uses nonlinear filters. The advectionreactiondiffusion equation,
(18) 
is defined with compatible with the exact solution (16). As in the first example, this equation also has zero Dirichlet boundary conditions in the same domain and it is discretized using second order finite difference with (centered difference is used for the first derivative). The discretization leads to the ODE system
(19) 
where represents the nonlinear terms and is the timedependent forcing term. The reference exact solution of this ODE system is obtained with high precision using the “ode45” routine in Matlab.
In this example, the protodecomposition is chosen with
and with . We remark that this choice of may not be the most efficient; nevertheless, this numerical example addresses the possible scenario where is nonlinear because the stiff term is embedded within a nonlinear term and, therefore, not easily identifiable. The tableau used in this example is the ARK548 also. The implicit equation, Eq. (8), is nonlinear and Newton’s method is used. The filter is defined as the element of the sequence generated by Newton’s method where followed by
and where, as before, . This is a standard Newton iteration where Jacobians are computed exactly and the linear system is solved with LU decomposition.
The numerical convergence rate is shown in Fig. 2. The convergence of SIMEXRK is unaffected by the choice of filter ; in contrast, IMEXRK needs either or to maintain timestep precision.
For the third example, we perform a numerical experiment to show that it is necessary to have a well defined smooth filter in order for SIMEXRK to work properly (this third example could be called a “counter example” instead). In this experiment, the setting is the same used for Eq. (17) (shown in Fig. 1) the only exception being the definition of the filter. A “purposely bad filter” is defined as being either equal to when is even in Algorithm 2 or equal to when
is odd. This means that either 2 or 3 Jacobi iterations are being alternately used to approximate the solution of the linear system at each RKstage. The result of this numerical experiment (not shown graphically) is the following: the convergence rate of SIMEXRK drops to fourth order instead of maintaining the fifth order shown in
Fig. 1(a). It is the RKstage alternation that causes the loss of precision. For completion, consider the following experiment: in Algorithm 2, use an “ok” filter where is either equal to when is even or when is odd; i.e., the number of Jacobi iterations alternates between steps rather than between RKstages. The result of this experiment (not shown graphically) is that SIMEXRK converges with fifth order just as shown in Fig. 1(a). In summary, if a filter is comprised of an iterative solver then the number of iterations of the iterative solver must remain fixed at every RKstage of each timestep.5 Timestep stability of SIMEX for PDEs
In this section we explore timestep stability of the SIMEXRK algorithm when applied to a stiff model PDE. The numerical stability is modified by the choice of tableau and filter.
For methods that do not decompose the ODE, i.e. monolithic methods, the method’s stability can be described by a region in the complex plane. The scalar model equation, with and , is used in order to depict the region of where the method yields an asymptotically decreasing sequence of values as timesteps progress. The analysis of this scalar model equation can be justified using the canonical form of a general linear ODE system with constant coefficients.
For SIMEX methods, the stability depends on many things: beside the choice of the IMEX tableau and the choice of the protodecomposition, the stability depends on the filter also. The scalar model equation is not useful for studying SIMEX stability because the influence of the filter is lost. Here we propose a generalization of the scalar model equation by considering the stability of the SIMEX method for an ODE system with the form
(20) 
where and
is a matrix with one eigenvalue is equal to
and (spectral radius). The stability region is defined as the subset of such that if, and only if, the numerical sequence generated by the timestepping method applied to Eq. Eq. 20 whith converges asymptotically to zero as for every unitsize initial condition . With this definition, the region defined by the scalar model equation is equal to .The matrix should represent a typical application of interest. Here we consider the equation
(21) 
where , and where is defined in the domain with periodic boundaries. This is a suitable model for addressing both diffusive and dispersive second order PDEs. The solutions are asymptoticaly stable when and oscillatory when is imaginary. The discretized Laplacian operator is obtained with the standard secondorder, five points, finitedifference stencil on a uniform grid. Let denote the number of discretization points along each edge of the domain and let be the spacing between grid points. We define the matrix from the discretized Laplacian according to . This multiplicative constant ensures the eigenvalues of are within the interval ; moreover, the largest eigenvalue of approaches the value 1 with the rate . For example, the largest eigenvalue of is . For the purpose of depicting the stability region, we can simply set in Eq. Eq. 20.
Henceforth, we study the stability of SIMEXRK methods by applying it to Eq. Eq. 20 with the matrix . We choose the protodecomposition of Eq. Eq. 20 with everything implicit: , and . When comparing Eq. Eq. 20 with the discretized version of Eq. Eq. 21, is related to according to the formula
(22) 
By changing the filter and the tableau, the stability region attains different sizes and shapes; this is shown in the figures that follow. These regions are computed by placing a grid on the complexplane and then applying the algorithm with for up to 30 steps for each value of on the grid. A random initial condition is chosen and, after 30 steps, the amplification factor of the last step is stored. The resulting data is then smoothed with a lowpass filter and the level curve of value 1 is plotted as the boundary of the region. For values of outside each depicted region, a random initial condition becomes larger as the steps increase. The regions shown here are obtained with . When using or larger, changes in these regions are very small.
Fig. 3 shows stability regions obtained using SIMEXRK with the CNH tableau and Fig. 4 using the ARK436 tableau. The filters used in these figures are the following: (i) is the defaultfilter. (ii) computes one iteration of an alternatedirection tridiagonal solver (ATS). (iii) applies five GaussSeidel iterations. (iv) computes 3 iterations of ATS. (v) approximately solves the linear system with an incompleteLU factorization (ILU). (vi) computes successivelyrelaxed GS iterations. The parameters defining these filters are detailed next.
Each ATS iteration used in solves two tridiagonal systems. The matrix of the linear system is decomposed as where is tridiagonal. Then, an ATS iteration is
where is the permutation of produced by interchanging the roles of the horizontal and the vertical directions in the natural ordering of nodes.
The ILU filter, , approximately solves the linear system with an incompleteLU factorization using the “ilu” routine in Matlab; the droptolerance is set to . This means any element less than 2% the size of the pivot is dropped.
Each iteration in “relaxes” one GS iteration where the relaxation parameter is set to ; this is actually an “underrelaxation”. This value of the relaxation parameter was tested against other values and it produced the largest stable interval along the real axis in the case of the CNH tableau. Different from other filters, is not based on a convergent linear solver. This filter has a bigger stability region than in the case of the CNH tableau but both filters have similar stability regions in the case of ARK436.
Several interesting features can be observed in Figs. 3 and 4. Primarily, it is visible that ARK436 obtains a larger stability region than CNH for the same filter. Of course, the computational effort placed in each step is different; ARK436 has five implicit RKstages while CNH has only one. The filter, region (5), produces a remarkably large stability region in ARK436; it reaches up to in the real axis. We should remark that is not a constanteffort filter; the fillin rate changes as changes because the droptolerance is being kept fixed in the ILU decomposition. Nevertheless, the fillin at is modest when compared with the matrix of the linear system : the number of nonzero entries in the lowertriangular matrix is times the number of nonzero entries in .
Another interesting feature is that a segment of the imaginary axis is inside the stability region for some filters in Fig. 4. For instance, the segment is contained in the stability region of the filter . Moreover, the regions associated with filters and contain the segment . This may be a consequence of the stable implicit scheme used in ARK436.
In SIMEXRK the size of the stability region is largely controlled by the choice of the filter. Combinations of tableau and filter pairs provide an immense range of possibilities to explore. The most appropriate combination depends on the particular PDE and on the precision to be attained; Eq. (21) should be viewed as a standard example where filters can be compared.
6 Spacetime grid convergence for a linear advectiondiffusion PDE
This section shows how SIMEXRK can be used for the solution of an advectiondiffusion equation in 2D where parabolic stiffness is predominant as the grid gets refined. Parabolic PDEs require ever larger stability regions when is being reduced while keeping the ratio fixed.
Here, we obtain numerical solutions of
(23) 
where is defined in the spatial domain with periodic boundaries and where is a constant unitary vector. The r.h.s. is consistent with the exact solution . All partial derivatives are discretized with fourthorder finitedifference formulas which use stencils with five points in each Cartesian direction. In this way, the discretized Laplacian has a 9point stencil with a cross shape. As before, denotes the number of discretization points along each Cartesian direction; the spatial grid is uniform with spacing between grid points. Using the method of lines, this leads to an ODE system
(24) 
where and and denote the matrices resulting from the discretization of the Laplacian and of the advection term respectively; is the forcing term where is evaluated at grid points. The protodecomposition is and . At every implicit RKstage, the matrix of the linear system has nonzero elements. For timeintegration, we use the fourth order ARK436 tableau so that the resulting scheme is both fourth order in time and space.
Numerical results are shown in Fig. 5 for seven different spacetime grids: ; for . The ratio of timestep size to gridspace is kept fixed, . Ten different filters were tested with SIMEXRK, these filters can be grouped in two categories: (i) GS filters, denoted by , which perform GaussSeidel iterations where ( being the defaultfilter); and (ii) preconditioned Conjugate Gradient Squared (CGS) filters, denoted by , which apply iterations of the CGS method using an ILU decomposition as preconditioner. The droptolerance for the ILU decomposition is . The number of CGS iterations is either or . The fillin is small when using this droptolerance: at the finest grid the number of nonzero elements in and in divided by is equal to ; this means an increase in memory use.
Fig. 5 shows the numerical error of the best filter for SIMEXRK at each grid. The best filter is the one with the smallest computational time among all filters with stable timeintegration up to . In this experiment, numerical instability is detected when any component of the numerical solution exceeds at any timestep along the algorithm. SIMEXRK data points are marked with a plus sign and joined by a dashed line. A line with slope 4 is placed in the diagonal of this figure. SIMEXRK’s errors follow closely the fourth order convergence rate.
For comparison, Fig. 5 brings the error of the “best filter” for IMEXRK also. The data points are marked with circles (). This data is obtained replacing the solver in IMEXRK by each one of the filters. As before, the best filter is the one with the least computational time among the filters with stable timeintegration. Of course, IMEXRK may become imprecise because each filter has a fixed number of iterations. Nevertheless, using this criterion we can distinguish whether a filter used in SIMEXRK is solving Eq. Eq. 8 accurately or not.
SIMEXRK reaches the correct numerical convergence rate even though Eq. Eq. 8 is not being accurately solved at some of these grids. At grids and of Fig. 5, SIMEXRK is more precise than IMEXRK while using the same number of iterations. This advantage can be attributed to the extra precision provided by RBD. Furthermore, at grids and , IMEXRK needs more iterations than SIMEXRK in order to reach a stable, albeit imprecise, solution. It may be a fortunate coincidence that the improved numerical precision provided by RBD results in a stable outcome for SIMEXRK with less iterations than for IMEXRK.
In order to show how different filters behave, part of the data of Fig. 5 is expanded in Fig. 6. Each marker in Fig. 6 shows the numerical error vs. CPU time for a given filter/method pair at a given grid. This figure shows the filters , , , and , when used with the grids for . Only the combinations that yield stable results are shown. SIMEXRK allows every filter to be used in more grids than IMEXRK. In particular, at grid , SIMEXRK can reach an accurate solution using either or ; in contrast, IMEXRK needs to use . At this grid, is slightly less precise than ; perhaps this is a sign of borderline stability which indicates that does not bear more grid refinements. At the finest grid, , IMEXRK needs to use more than two CGS iterations in order to reach an accurate solution.
We remark that, in view of Algorithm 3, any stopping criterion used in IMEXRK can be used as a stabilization criterion in SIMEXRK also. According to the numerical experiments shown here, the stabilization criterion used with SIMEXRK is likely to be less strict than the stopping criterion used with IMEXRK.
7 Conclusion
In this article the residual balanced decomposition (RBD) is introduced. Given an implicitexplicit decomposition and a filter, RBD provides a new decomposition where the implicit equation is solved exactly without additional computational work. This remarkable property allows great freedom for exploring the numerical properties of the resulting algorithm. The decomposition itself is rather abstract, Eqs. Eqs. 11 and 10, but the computational steps for its evaluation are straightforward, Eq. Eq. 13. Here, RBD is used with IMEXRK methods that have ESDIRK implicit schemes, the resulting method is the SIMEXRK method, Algorithm 2.
One way to construct a filter is to arbitrarily fix the number of iterations of an iterative solver applied to the implicit equation of the original IMEX method, Eq. Eq. 8. Computational experiments confirm that SIMEXRK maintains the correct numerical convergence rate even though its filter yields an inaccurate solution for Eq. Eq. 8. Thus, the goal of the filter is to ensure timestep stability and not an accurate solution of Eq. Eq. 8. Numerical experiments in Section 6 show a spatialtemporal grid refinement of an advectiondiffusion PDE. In these experiments SIMEXRK can produce an accurate solution as soon as its filter is strong enough to yield a stable outcome.
In order to study the stability of the SIMEXRK method, the stabilityregion in the complexplane is generalized to ODE systems where the largest eigenvalue of the matrix is equal to one. As shown in Section 5, the stability region can vary significantly depending on the combination of RK scheme and filter being used. The size of the stability region clearly depends on the computational effort placed in the filter. This region can reach far into the negative real axis as observed with the ILU filter. Also, when the implicit RK scheme is Lstable, SIMEXRK can produce stability regions that contain a significant interval of the imaginary axis.
SIMEX can be extended to multistep (MS) methods. This occurs because equations Eq. 12 through Eq. 14 are valid as long as the implicit equation remains the same at each step of the MS method. A careful initialization procedure is necessary when using SIMEXMS method (an example is available as supplementary material). RBD can be used with some types of generallinearmethods (GLM) also. One such GLM is found in a recent study [31]. The RK part of the GLM method needs to have the singlydiagonallyimplicit property; such property ensures that the same implicit equation appears at every stage.
The search for performance improvements with SIMEXRK can start with Algorithms 3 and 2 where the sequence of filters and the stabilization criterion are provided by the incumbent IMEXRK method. At this point, it is likely that SIMEXRK will closely replicate the incumbent IMEXRK method. Then, SIMEXRK may be improved by relaxing the stabilization criterion. Further improvements may result either from using less expensive preconditioners or from constructing preconditioners that target specific subspaces where the numerical instability is present. SIMEX also opens the possibility of using filters that are not based on iterative solvers.
References
 [1] X. Antoine, C. Besse, and V. Rispoli, Highorder imexspectral schemes for computing the dynamics of systems of nonlinear Schrodinger/GrossPitaevskii equations, J. Comp. Phys., 327 (2016), pp. 252–269.

[2]
U. M. Ascher, S. J. Ruuth, and R. J. Spiteri,
Implicitexplicit RungeKutta methods for timedependent partial differential equations
, Applied Numerical Mathematics, 25 (1997), pp. 151–167.  [3] U. M. Ascher, S. J. Ruuth, and B. T. R. Wetton, Implicitexplicit methods for timedependent partial differential equations, SIAM J. Numer. Anal., 32 (1995), pp. 797–823.
 [4] S. Blaise, J. Lambrechts, and E. Deleersnijder, A stabilization for threedimensional discontinuous Galerkin discretizations applied to nonhydrostatic atmospheric simulations, International Journal for Numerical Methods in Fluids, 81 (2016), pp. 558–585.
 [5] S. Boscarino, L. Pareschi, and G. Russo, Implicitexplicit RungeKutta schemes for hyperbolic systems and kinetic equations in the diffusion limit, SIAM J. Sci. Comput., 35 (2013), pp. A22–A51.
 [6] J. H. Chaudhry, D. Estep, V. Ginting, J. N. Shadid, and S. Tavener, A posteriori error analysis of IMEX multistep time integration methods for advectiondiffusionreaction equations, Computer Methods in Applied Mechanics and Engineering, 285 (2015), pp. 730–751.

[7]
G. J. Cooper and A. Sayfy,
Additive methods for the numerical solution of ordinary differential equations
, Math. Comp., 35 (1980), pp. 1159–1172.  [8] G. J. Cooper and A. Sayfy, Addtive RungeKutta methods for stiff ordinary differential equations, Math. Comp., 40 (1983), pp. 207–218.
 [9] M. Crouzeix, Une méthode multipas impliciteexplicite pour l’approximation des équations d’evolution paraboliques, Numer. Math., 35 (1980), pp. 257–276.
 [10] D. R. Durran and P. N. Blossey, Implicitexplicit multistep methods for fastwaveslowwave problems, Monthly Weather Review, 140 (2012), pp. 1307–1325.
 [11] B. Froehle and P.O. Persson, A highorder discontinuous Galerkin method for fluidstructure interaction with efficient implicitexplicit time stepping, J. Comp. Phys., 272 (2014), pp. 455–470.
 [12] F. Garcia, L. Bonaventura, M. Net, and J. Sánchez, Exponential versus imex highorder time integrators for thermal convection in rotating spherical shells, J. Comp. Phys., 264 (2014), pp. 41––54.
 [13] F. Garcia, M. Net, B. GarcíaArchilla, and J. Sánchez, A comparison of highorder time integrators for thermal convection in rotating spherical shells, J. Comp. Phys., 229 (2010), pp. 7997––8010.
 [14] D. Ghosh and E. Constantinescu, Semiimplicit time integration of atmospheric flows with characteristicbased flux partitioning, SIAM Journal on Scientific Computing (SISC), 38 (2016), pp. A1848–A1875, doi:10.1137/15M1044369.
 [15] I. Higueras and T. Roldán, Construction of additive semiimplicit Runge–Kutta methods with lowstorage requirements, J. Sci. Comput., 67 (2016), pp. 1019–1042.
 [16] A. Kanevsky, M. H. Carpenter, D. Gottlieb, and J. S. Hesthaven, Application of implicitexplicit high order RungeKutta methods to discontinuousGalekin schemes, J. Comp. Physics, 225 (2007), pp. 1753–1781.
 [17] V. KazemiKamyab, A. van Zuijlen, and H. Bijl, Analysis and application of high order implicit rungekutta schemes for unsteady conjugate heat transfer: A stronglycoupled approach, J. Comput. Physics, 272 (2014), pp. 471–486.
 [18] C. A. Kennedy and M. H. Carpenter, Additive RungeKutta schemes for convectiondiffusionreaction equations, Technical report NASA/TM2001211038, (2001).
 [19] C. A. Kennedy and M. H. Carpenter, Diagonally implicit RungeKutta methods for ordinary differential equations. a review, Technical report NASA/TM2016219173, (2016).
 [20] C. A. Kennedy, M. H. Carpenter, and R. M. Lewis, Additive RungeKutta schemes for convectiondiffusionreaction equations, Appl. Numer. Math., 44 (2003), pp. 139–181.
 [21] P. Mati, M. A. Calkins, and K. Julien, A computationally efficient spectral method for modeling core dynamics, Geochemeistry, Geophysics, Geosystems, 17 (2016), pp. 3031–3053.
 [22] P.O. Persson, Highorder les simulations using implicitexplicit rungekutta schemes, in 49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, AIAA 2011684, (2011), doi:10.2514/6.2011684.
 [23] M. Restelli and F. X. Giraldo, A conservative discontinuous Galerkin semiimplicit formulation for the NavierStokes equations in nonhydrostatic mesoscale modeling, SIAM J. Sci. Comput., 31 (2009), pp. 2231–2257.
 [24] D. Ruprecht and R. Speck, Spectral deferred corrections with fastwave slowwave splitting, SIAM J. Sci. Comput., 38 (2016), pp. A2535–A2557.
 [25] M. P. Ueckermann and P. F. J. Lermusiaux, Hybridizable discontinuous Galerkin projection methods for NavierStokes and Boussinesq equations, J. Comp. Physics, 306 (2016), pp. 390–421.
 [26] A. H. van Zuijlen and H. Bijl, Implicit and explicit higher order time integration schemes for structural dynamics and fluidstructure interaction computations, SIAM J. Sci. Comput., 38 (2016), pp. A1430–A1453.

[27]
H. Wang, C.W. Shu, and Q. Zhang,
Stability and error estimates of local discontinuous galerkin methods with implicitexplicit timemarching for advectiondiffusion problems
, SIAM J. Numer. Anal., 53 (2015), pp. 206–227.  [28] H. Wang, C.W. Shu, and Q. Zhang, Stability analysis and error estimates of local discontinuous galerkin methods with implicit–explicit timemarching for nonlinear convection–diffusion problems, Applied Mathematics and Computation, 272 (2016), pp. 237–258.
 [29] H. Weller, S.J. Lock, and N. Wood, Runge–Kutta IMEX schemes for the horizontally explicit/vertically implicit (HEVI) solution of wave equations, J. Comp. Physics, 252 (2013), pp. 365–381.
 [30] M. Zayernouri and A. Matzavinos, Fractional AdamsBashforth/Moulton methods: An application to the fractional KellerSegel chemotaxis system, J. Comp. Phys., 317 (2016), pp. 1–14.
 [31] H. Zhang, A. Sandu, and S. Blaise, High order implicitexplicit general linear methods with optimized stability regions, SIAM J. Sci. Comput., 38 (2016), pp. A1430–A1453.