# Time-Accurate and highly-Stable Explicit operators for stiff differential equations

Unconditionally stable implicit time-marching methods are powerful in solving stiff differential equations efficiently. In this work, a novel framework to handle stiff physical terms implicitly is proposed. Both physical and numerical stiffness originating from convection, diffusion and source terms (typically related to reaction) can be handled by a set of predefined Time-Accurate and highly-Stable Explicit (TASE) operators in a unified framework. The proposed TASE operators act as preconditioners on the stiff terms and can be deployed to any existing explicit time-marching methods straightforwardly. The resulting time integration methods remain the original explicit time-marching schemes, yet with nearly unconditional stability. The TASE operators can be designed to be arbitrarily high-order accurate with Richardson extrapolation such that the accuracy order of original explicit time-marching method is preserved. Theoretical analyses and stability diagrams show that the s-stages sth-order explicit Runge-Kutta (RK) methods are unconditionally stable when preconditioned by the TASE operators with order p ≤ s and p ≤ 2. On the other hand, the sth-order RK methods preconditioned by the TASE operators with order of p ≤ s and p > 2 are nearly unconditionally stable. The only free parameter in TASE operators can be determined a priori based on stability arguments. A set of benchmark problems with strong stiffness is simulated to assess the performance of the TASE method. Numerical results suggest that the proposed framework preserves the high-order accuracy of the explicit time-marching methods with very-large time steps for all the considered cases. As an alternative to established implicit strategies, TASE method is promising for the efficient computation of stiff physical problems.

## Authors

• 3 publications
• 4 publications
• 1 publication
• ### Parallel implicit-explicit general linear methods

High-order discretizations of partial differential equations (PDEs) nece...
02/03/2020 ∙ by Steven Roberts, et al. ∙ 0

• ### High-order partitioned spectral deferred correction solvers for multiphysics problems

We present an arbitrarily high-order, conditionally stable, partitioned ...
09/03/2019 ∙ by Daniel Z. Huang, et al. ∙ 0

• ### Fast IMEX Time Integration of Nonlinear Stiff Fractional Differential Equations

Efficient long-time integration of nonlinear fractional differential equ...
09/09/2019 ∙ by Yongtao Zhou, et al. ∙ 0

• ### Stable and accurate numerical methods for generalized Kirchhoff-Love plates

Efficient and accurate numerical algorithms are developed to solve a gen...
08/04/2020 ∙ by Duong T. A. Nguyen, et al. ∙ 0

• ### On the impact of explicit or semi-implicit integration methods over the stability of real-time numerical simulations

Physics-based animation of soft or rigid bodies for real-time applicatio...
11/20/2013 ∙ by Teodor Cioaca, et al. ∙ 0

• ### An efficient explicit approach for predicting the Covid-19 spreading with undetected infectious: The case of Cameroon

This paper considers an explicit numerical scheme for solving the mathem...
05/21/2020 ∙ by Eric Ngondiep, et al. ∙ 0

• ### On the explicit two-stage fourth-order accurate time discretizations

This paper continues to study the explicit two-stage fourth-order accura...
07/06/2020 ∙ by Yuhuan Yuan, et al. ∙ 0

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

Many scientific and engineering problems can be modeled as partial differential equations (PDEs) or ordinary differential equations (ODEs), and solved numerically with the assistance of modern computations

Dafermos and Mathématicien (2000); Anderson and Wendt (1995); Kim et al. (1987). The success of solving these problems strongly depends on the accurate, efficient, and stable numerical algorithms for both temporal and spatial discretizations. These are however not trivial due to the possible presence of stiffness originating from both numerics and physics Curtiss and Hirschfelder (1952); Lomax et al. (2013); Mortazavi et al. (2015). Numerical stiffness typically induced by regions of excessive spatial resolutions imposes that the time-step should be chosen even smaller than the smallest physical time scale of the system. A common situation is when structured meshes are used in cylindrical coordinates leading to excessive azimuthal resolution near the axis of symmetry. Physical stiffness occurs when there is a separation of scales between a fast (quasi-)steady physical phenomenon and slower unsteady physical processes with larger time scales. In both cases, the common feature is that the numerical solution exhibits fast decaying modes that the user is not interested in.

In wide range of multi-physics systems the problematic stiff terms are not known ahead of designing the numerical algorithm, e.g. due to nonlinearities in the governing differential equations, or because of the lack of prior insights about the phenomena to be simulated such as the identification of competing physics. In such cases it would be convenient to first develop a fully explicit numerical solver since explicit codes are often simple to program and debug. In other cases, one may have quickly modified an already available explicit code to simulate a new problem that is later identified to be stiff. The stiff processes frequently originate from only a few terms in the governing equations that compete in a quasi-steady manner. These terms can be identified by examining the impact of each term on the behavior of the solution and their characteristic time scales. In these scenarios users often desire to insert a quick fix to handle these specific stiff terms, while maintaining the structure of an existing verified code, rather than re-programing the entire code using traditional implicit schemes. TASE operators are particularly useful as quick programming remedies for such scenarios, as in many of these conditions the programming time may overwhelm the simulation time itself.

For temporal discretization for stiff PDEs and ODEs, unconditionally stable implicit time-marching methods are currently a common choice as they enable a flexible choice of the time-step size based on the desired accuracy Dahlquist (1963); Ascher et al. (1995); Boom and Zingg (2018). Taking aerodynamic applications for instance, the quantities of interests, e.g. drag and lift, are only relevant to the steady solution of Euler or Navier-Stokes (NS) equations Sclafani et al. (2008); Tinoco et al. (2017); Rumsey et al. (2018). For high computational efficiency, the nonlinear convection and diffusion terms of the governing laws can be handled implicitly, linearized in time with first-order approximation, and marched with the efficient Lower-upper symmetric-Gauss-Seidel (LU-SGS) Yoon and Jameson (1988); Chen and Wang (2000); Zhang and Wang (2004) method, which eliminates the need for block-diagonal inversions without using a diagonalization procedure. This first-order accurate implicit framework allows for stable computations of flows across a wide range of Mach numbers without time-step constraints. To recover the temporal accuracy for unsteady flows, a second-order implicit time-marching method, which is both A- and L-stable Dahlquist (1963), can be developed by further combing the dual time-stepping strategy Jameson (1991, 2015) with the second-order backward difference formula (BDF2). Another common numerical challenge is to accurately simulate the time dynamics of reactive flows with stiff chemistry, for which the time scale is determined by the physico-chemical process related to reaction and transport. In cases where the smallest chemical time scale is comparable to that of transport process, the coupled system can be marched simultaneously via explicit methods. In general however, the chemical time scale is orders of magnitude smaller than the transport-related time scale, therefore an explicit treatment of chemical terms leads to prohibitively small computational time steps and unaffordable computational costs. A popular strategy for this type of convection-diffusion-reaction equations is to treat the non-stiff convection-diffusion operators explicitly and the reaction-related stiff source terms implicitly by a class of operator-splitting methods Marchuk (1968); Strang (1968). The resulting ODEs for the reaction steps can be integrated efficiently by fully implicit multi-step variable-order BDF methods Gear (1971) or Rosenbrock-Krylov methods Tranquilli and Sandu (2014). Although these splitting methods are strongly stable Knio et al. (1999); Singer et al. (2006), the overall temporal accuracy is typically limited to second order Speth et al. (2013).

As it is well known, for scientific computations with temporal accuracy requirements, higher-order methods allow for reaching a target accuracy at a lower overall computational cost, despite a larger computational cost per time step. Consequently, many efforts have been devoted to the development of unconditionally stable time-marching methods of high order (higher than second order), e.g. the cyclic and composite linear multistep methods Donelson III and Hansen (1971); Sloate and Bickart (1973), the advanced-step-point methods Psihoyios and Cash (1998); Aiken (1985); Cash (1983), and hybrid methods England (1982). However, A-stable linear multi-step schemes are at most second-order accurate Dahlquist (1963). On the contrary, it is possible to derive both A- and L-stable implicit RK methods that yield higher-order accuracy with strong stability Butcher (1987, 2003, 1964)

. Implicit RK methods can be classified based on the form of their coefficient matrix (also called “Butcher tableaux”). Specifically, the implicit RK methods characterized by a lower triangular coefficient matrix with nonzero diagonal elements are called diagonally-implicit RK method (DIRK)

Alt (1971); Crouzeix (1975); Alexander (1977). Compared to those with full coefficient matrix, the solution at each RK substep can be solved sequentially with high efficiency as it is independent of subsequent stage values Alexander (1977). The efficiency can be further improved by the singly-diagonally implicit RK (singly-DIRK) methods Wanner and Hairer (1996), for which the diagonal coefficients of the lower triangular matrix are identical and the operations on the Jacobian can be reused for each RK substep throughout these iterations without any significant reduction in convergence when the Newton’s method is employed to solve the resulting nonlinear system Crouzeix (1975).

As another class of efficient time-marching methods for convection-diffusion equations, the implicit-explicit (IMEX) Ascher et al. (1997); Calvo et al. (2001); Ascher et al. (1995); Cooper and Sayfy (1983); Wang et al. (2015) methods consist of employing an implicit discretization for the diffusion terms and an explicit discretization for the convection terms. The rationale is that the diffusion terms are typically stiff and linear while the convection terms are not stiff but nonlinear. It is quite challenging to handle the nonlinear convection terms implicitly as the resulting nonlinear system needs to be solved by an iterative solver. When combined with the high-order linear multistep schemes Donelson III and Hansen (1971); Sloate and Bickart (1973), the IMEX methods still suffer from strong time-step constraint Ascher et al. (1997); Calvo et al. (2001). Instead, the IMEX methods are typically developed based on the high-order RK method with larger stability regimes Ascher et al. (1997); Wang et al. (2015). An -stages high-order IMEX-RK method Ascher et al. (1997); Calvo et al. (2001); Wang et al. (2015) can be defined by deploying an -stages DIRK Alexander (1977) method for diffusion terms and an -stages explicit RK method for convection terms. IMEX time-marching methods have been widely used to solve different problems Zhong (1996); Pareschi and Russo (2005); Gottlieb and Wang (2012); Kennedy and Carpenter (2003); Wang et al. (2016); Fu and Shu (2017).

In the present work, we propose an alternative framework to implicitly handle stiff ODEs or PDEs. A family of unique TASE operators with arbitrarily high-order accuracy is defined that can be deployed to any existing explicit time-marching methods. The resulting time integration remains explicit, yet it becomes nearly unconditional stable. The stiffness originating from different mechanisms, including convection, diffusion, and source terms, can be treated implicitly in the unified TASE framework. Moreover, the order of accuracy of the baseline explicit time-marching methods can be maintained. The stability properties of the TASE operators in combination with established RK methods is analyzed by theoretical proof and stability diagrams. The performance of TASE operators in terms of accuracy and stability are demonstrated by a set of canonical benchmark simulations.

The rest of the paper is organized as follows. The TASE framework is introduced in Section 2,which also includes the analytical description of a family of operators. The stability properties of the proposed TASE operators are analyzed in Section 3. Section 4 further discuss the practical deployment of TASE operators to challenging physical problems with source terms, multiple stiff mechanisms, non-linear equations, and grid-induced stiffness in polar coordinates. Concluding remarks and future work are discussed in the last section.

## 2 TASE operators

The problem of formulating TASE operators is described in Section 2.1. A specific example of a first-order solution is described in Section 2.2, and a family of generalized TASE operators is introduced in Section 2.3. The free parameter, the explicit definitions, and the stability diagrams of TASE operators are discussed in Sections 2.42.6.

### 2.1 Problem definition

Typical theoretical descriptions of physical and engineering problems involve conservation equations represented as

 ∂y∂t=f(x,t,y(x,t)), (1)

where denote space and time, respectively,

is a vector of unknown physical variables, and

is a function that describes the time-rate of change of . The numerical solution of Eq. (1) is generally sought on a discrete grid , whose resolution ought to be chosen carefully for the numerical solution to faithfully approximate the exact solution. This is a major requirement of any physics-based computational solver. The spatial mesh resolution should resolve all the physical length scales of the problem. These length scales are determined by the the problem geometry, the initial condition, and the physical features that may arise over time. After spatial discretization the continuous conceptual model is turned into a semi-discrete computer model

 dYdt=F(X,t,Y(t)), (2)

where the numerical approximation of the operator is denoted as , and represents a vector of spatially discretized unknowns. In order to obtain the numerical solution to Eq. (2) at different time, the time derivative is numerically approximated with a certain time-marching scheme and time step . The temporal resolution is constrained by both physical and numerical requirements.

Consider the linear model problem for the semi-discrete ordinary differential equation given in Eq. (2)

 dYdt=LY, (3)

where is a vector of unknowns, and is a matrix operator. The explicit time integration of Eq. (3) is subject to an upper-bound constraint for the maximum time step that can be used for the numerical solution to be stable, as described in Section 1. The objective of the present paper is to propose a family of operators such that the explicit time integration of

 dYdt=(T(p)LL)Y (4)

is unconditionally or nearly unconditionally stable irrespective of the choice of the time-marching scheme, and requiring that

 T(p)L=1+O(Δtp), (5)

which ensures that the exact solution to Eq. (4) converges to the exact solution of Eq. (3) with integer order in the limit of small time steps. In other words, the goal is to enforce that the numerical solution obtained with the TASE framework collapses to that that would be obtained with a unmodified explicit scheme in the limit of small time-steps, consistently with the fact that it is likely that the explicit time integration is stable on its own in this limit. Additionally, when the order is chosen to match the order of accuracy of the time-stepping scheme, the global order of accuracy of the numerical solver is not affected by replacing in the code with as in Eq. (4).

The TASE operator can be viewed as a preconditioner that transforms the exact operator into a modified operator

 T(p)L=T(p)LL, (6)

such that . For various reasons highlighted in the following sections, the analyses of the TASE framework can be more clearly clarified using , therefore for the rest of the paper we mostly retain in our notation.

### 2.2 Example of a first-order TASE operator

The time-marching of Eq. (3), , with implicit Euler is unconditionally stable, and yields the fully-discrete computer model

 Yn+1−YnΔt=LYn+1, (7)

where and . Using the matrix equality , and assuming that is invertible, Eq. (7) can be re-arranged as

 Yn+1−YnΔt=((1−ΔtL)−1L) Yn. (8)

Equation (8) resembles the fully-discrete equation that would be obtained by time-advancing Eq. (4), with explicit Euler, yet with unconditionally stable property by construction. Term by term identification of both equations suggest an example of first-order solution to the TASE problem definition formulated in Section 2.1, i.e.

 T(1)L=(1−ΔtL)−1. (9)

The time-marching of Eq. (4), , using the TASE operator defined in Eq. (9) is unconditionally stable if the time-marching scheme is explicit Euler, by construction. Furthermore, it is trivially first-order since its Laurent series is in the limit of small time step Markushevich (1977). This first-order operator satisfies only partly the design constraints described in Section 2.1, because its stability properties for other time-marching schemes are still unknown, and it is only first-order. The extension of this operator to general explicit time-marching schemes of arbitrary accuracy order is described in the next section.

### 2.3 Generalization and extension to arbitrary order

The TASE operator of order is recursively defined as

 T(p)L(α)=⎧⎪ ⎪⎨⎪ ⎪⎩(1−αΔtL)−1,if p=12p−1T(p−1)L(α/2)−T(p−1)L(α)2p−1−1,if p≥2 (10)

where is a free parameter that depends on the time-marching scheme and must satisfy

 α≥αmin=2p−1C. (11)

The denominator is the maximum value of that guarantees a stable numerical solution for a given explicit time-marching scheme when applied to the model problem , with and real unknown Eberly (2008); Lambert (1991). It can also be interpreted as the absolute value of the intersection of the stability diagram of a given explicit time-marching scheme with the negative real axis, which is a often known value for common schemes. Depending on the time-marching scheme and desired TASE accuracy order, the order-unity minimum value of can be easily obtained from Eq. (11). For instance, the values of for common explicit RK methods preconditioned by the TASE operators of varying order are given in Table 1. The motivation for the constraint given by Eq. (11) on the choice of becomes clearer when analyzing the stability properties of the defined TASE operators in Section 3.

The specific first-order operator derived in Section 2.2 corresponds to , where is greater than the threshold value given in Table 1 for the combination of the explicit Euler scheme with the first-order TASE operator. However, when deployed with a high-order time-marching scheme, it would reduce the global order of accuracy of the solver to one. For particular scenarios, where the steady solution is pursued, this degeneration of accuracy order is fine. However, for most common cases with temporal accuracy requirement, it is yet desirable that the accuracy order of explicit time-marching methods is preserved after preconditioned with the TASE operator. The family of TASE operators defined in Eq. (10) generalizes the example given in Section 2.2 to arbitrarily high order via Richardson extrapolation, and a broader class of time-marching schemes via the parameter .

The proof of the statement of accuracy for the TASE operator given in Eq. (5) is straightforward due to the recursive application of Richardson extrapolation in Eq. (10) Richardson (1911); Richardson and Gaunt (1927). Note that Richardson extrapolation has been previously employed to improve the accuracy order of the strongly stable -method to solve stiff ODEs implicitly Havasi (2010). One secondary yet useful outcome of examining the Taylor series of Eq. (5) and Eq. (10) is the conclusion that the leading error term increases with , which indicates that is a good choice in practice.

### 2.4 Rationale for the minimum value of α

A necessary condition for the time integration of

to be unconditionally stable is that the modified eigenvalue

falls within the stability diagram of the given explicit time-marching scheme as the time step goes to infinity. We show below that this necessary condition collapses to the constraint prescribed on in Eq. (11),

###### Theorem 2.1 (Asymptotic stability)

The arbitrary order explicit time integration of remains stable in the limit of very large time steps when (necessary condition for unconditional stability).

It follows from the Laurent series of in the limit of large time steps and recursive application of Eq. (10), that

 (T(p)λ(α)λ)Δt=−2p−1α+O(Δt−1). (12)

Using the definition of given in Eq. (11), the right-hand side term becomes . In conjunction with Eq. (12), it implies that in the limit of large ,

 −C≤(T(p)λ(α)λ)Δt<0. (13)

Equation (13) proves that in the limit of large time steps , the modified dimensionless eigenvalue falls on the real axis and within the stability diagram of all time-marching schemes whose stability diagram cross the real axis at . This ensures stability in the asymptotic limit of large , which concludes the proof of theorem (2.1

). The above asymptotic analysis shows that the prescribed real value of

given in Eq. (11) is not compatible with time-marching schemes whose stability diagram is restricted to purely imaginary . These schemes are excluded from the discussion in the present paper anyway, since they are only useful for purely oscillatory physics which require to be time resolved and consequently rule out the need for their implicit treatment. In Section 3, we will show that the necessary condition given in the present section for works well in practice.

### 2.5 Explicit definition of TASE operators

The formulation of TASE operator only includes the parameter , the time step and the operator . The explicit expressions can be obtained by applying the recurrence formula prescribed in Eq. (10). For the convenience of readers, the TASE preconditioners and operators are explicitly defined up to order as belows.

The TASE preconditioners are expressed as

 T(p)L(α)=p−1∑k=0βp,k(2k−αΔtL)−1 (14)

where the coefficients are given in Table 2.

The TASE operators are expressed as

 T(p)L(α)=1αΔt(−(2p−1)+p−1∑k=0γp,k(2k−αΔtL)−1) (15)

where the coefficients are given in Table 3.

In terms of the implementation of the TASE operator, , or equivalently , does not need to be computed explicitly, rather it should be avoided since it is often a full matrix. The only required routine is the one that computes , or equivalently . For many practical operators , this leads to a banded matrix to solve which is much less costly than directly computing the inverse of a matrix Banerjee and Roy (2014). In general, the direct evaluation of requires computing the inverse of matrices, which makes the cost of the present method at least comparable to that of other established implicit schemes. Nevertheless, the computational savings allowed by the use of a larger time step can outweigh the additional cost induced by the TASE operators at each time step.

### 2.6 Stability diagrams

The stability diagrams for a set of -stages th-order RK methods preconditioned by TASE operators up to fourth-order are given in Fig. 1 and Fig. 2. As shown in Fig. 1, the stability regions of the explicit RK methods are substantially enlarged by the corresponding TASE operators when the parameter is set as the minimum value defined in Table. 1. Specifically, all the RK methods preconditioned by the first- and second-order TASE operators are unconditionally stable, i.e. the entire left-half complex plane is mapped into the interior of the stability regions of the original explicit RK methods with TASE operators. On the other hand, the preconditioned RK methods become nearly unconditionally stable when combined with the TASE operators of order higher than second, as revealed in the zoomed-in Fig. 2. Nevertheless, these unstable regions are tiny and correspond to operators strongly dominated by convection, as illustrated by the axes range in Fig. 2 that shows two orders of magnitude larger imaginary eigenvalue compared to their real counterparts. This class of problems is not of primary concern in the present analysis since stiffness originating from convection alone cannot be alleviated by the use of large time steps. Therefore, the nearly unconditionally stable RK methods may still allow for very-large time step in practice and will be demonstrated in the following numerical validations in Section 4.

In order to investigate the dependence of the stability properties of the preconditioned RK methods on the choice of , the zoomed-in stability diagrams of the fourth-order RK method preconditioned by the TASE operators of first- to fourth-order with varying are plotted in Fig. 3. It is shown that all the preconditioned RK methods have large unstable regions in the left-half complex plane when . Meanwhile, when , the stability of the preconditioned RK methods is improved with larger , i.e. the tiny unstable region if it exists in the left-half complex plane becomes even smaller.

## 3 Analyses of stability properties

In this section, the stability properties of the proposed TASE operators applied to several established explicit time-marching methods are analyzed in detail. The notion of stability chosen in this paper is the traditional notion of A-stability, which enforces that the numerical method be stable when the exact numerical solution is stable too, as defined below Havasi (2010); Wanner and Hairer (1996).

###### Definition 3.1 (A-stability)

Consider the set containing all values of for which . A method with stability function is said to be A-stable if .

This section proves that the family of TASE operators defined in Section 2.3 satisfies the constraint that the time-integration of Eq. (4), , is either unconditionally or nearly unconditionally stable. The analysis considers a single eigenvalue of , and its TASE substitute . It is sufficient to do the analysis for individual eigenvalues of since the result extends to general matrices

by eigenvalue decomposition, as the eigenvectors basis of modified operator

is the same as that of the original operator .

In particular, we restrict the analysis to general -stages th-order RK methods, written in Shu-Osher form Shu (2002):

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩y(0)=yn,y(i)=i−1∑k=0(ai,ky(k)+Δtbi,kf(y(k),tn+ckΔt)),   ai,k≥0,   i=1,…,syn+1=y(s), (16)

where by consistency . We additionally assume that all ’s are non-negative, i.e. , which is the case for most RK methods including -stages th-order with ones of interest in the present study. For these specific methods, the stability diagram does not depend on the coefficients themselves, only on the order, or equivalently the number of stages, therefore there is no need to specify the exact coefficients here.

### 3.1 Sufficient condition for unconditional stability of first-order TASE operators with arbitrary order RK

We provide below the theoretical proof that the combination of first-order TASE operators with arbitrary order RK methods (16) is unconditionally stable if is chosen larger than a critical value.

###### Theorem 3.2 (Unconditional stability of first-order TASE with RK)

The explicit RK time integration (16) of is unconditionally stable when (sufficient condition).

We prove theorem 3.2 in two steps. We first prove the particular case of explicit Euler time integration:

###### Lemma 3.3 (First-order TASE: Explicit Euler)

The explicit Euler time integration of is unconditionally stable, i.e. , if and only if .

The explicit Euler time integration of yields . Using the definition of in Eq. (9), it can be re-arranged as

 yn+1−ynΔt=αλyn+1+(1−α)λyn. (17)

The corresponding stability function is

 σ=1+(1−α)λΔt1−αλΔt. (18)

It is trivial to show that the condition subject to is satisfied if and only if . This concludes the proof of lemma 3.3.

We now prove that unconditional stability for higher-order RK schemes directly follows from unconditional stability with explicit Euler:

###### Lemma 3.4 (First-order TASE: From explicit Euler to arbitrary order RK)

The arbitrary order RK time integration (16) of is unconditionally stable if the explicit Euler integration of where is unconditionally stable.

The explicit time integration of with a arbitrary order RK method (16) yields

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩y(0)=yn,y(i)=i−1∑k=0ai,k(y(k)+ΔtT(1)λ(α)bi,kai,kλy(k)),   ai,k≥0,   i=1,…,syn+1=y(s), (19)

The key is to note from the definition of the first-order TASE operator in Eq. (9) that

 T(1)λ(α)bi,kai,kλ=T(1)λ′(αi,k)λ′ (20)

where and . Importantly, note that since , if and only if . Therefore Eq. (19) becomes

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩y(0)=yn,y(i)=i−1∑k=0ai,k(y(k)+ΔtT(1)λ′(αi,k)λ′y(k)),   αi,k≥0,   i=1,…,syn+1=y(s), (21)

It follows from Eq. (21) that all intermediate stages of arbitrary order RK schemes of the form (16) are convex combinations of explicit Euler operations, with combination factors summing to 1, eigenvalue replaced by , and replaced by . Therefore and for . If we assume that the explicit Euler integration of is unconditionally stable, it trivially follows that , which concludes the proof of lemma 3.4.

The complete proof of theorem 3.2 is obtained by combining both lemmas. From lemma 3.4, we conclude that for the explicit RK time integration (16) of to be unconditionally stable, it is sufficient that the explicit Euler integration of where is unconditionally stable. From lemma 3.2, we conclude that the explicit Euler integration of is unconditionally stable if only if for all and , which is equivalent to . This concludes the proof of theorem 3.2.

In this section, we have proven a sufficient condition for the unconditional stability of the combination of first-order TASE with arbitrary order RK. Given a specific RK method, one can write it in Shu-Osher form (16) and compute the critical value . For example, linear SSP RK methods up to order 4 Gottlieb et al. (2001) are characterized by hence irrespective of the order. Interestingly the resulting sufficient condition is more restrictive than the necessary condition derived in Section 2.4 and shown in Table 1 for RK3 and RK4. Stability diagrams shown in Section 2.6, and numerical examples in Section 4 indicate that in practice the necessary condition provided in provided in Section 2.3 and in Table 1 works well. We further assess this in the proof below in Section 3.2, for first-order TASE but also for higher-order.

### 3.2 General proof

The analysis conducted in Sec. 3.1 becomes cumbersome when the order of the TASE operator is larger than one. Here instead, we proceed with an alternative approach that relies on the following corollary Havasi (2010)

###### Corollary 3.4.1 (A-stability proof from maximum modulus principle)

A numerical method is A-stable if and only if
(a) is analytic in , and
(b) it is stable on the imaginary axis, i.e. for all real values of .

The proof is made in two steps. First, we prove that is an analytic function in , then we assess the stability of the methods for purely imaginary eigenvalues.

(a) We first compute the stability function with . The explicit time integration of with arbitrary order TASE and arbitrary order RK method yields

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩y(0)=yn,y(i)=i−1∑k=0ai,k(y(k)+ΔtT(p)λ(α)bi,kai,kλy(k)),   ai,k≥0,   i=1,…,syn+1=y(s), (22)

which can be recasted in the more explicit form

 yn+1=(1+s−1∑k=0As,k(ΔtT(p)λ(α)λ)k+1)yn, (23)

where can be recursively defined using linear combinations of the coefficients and Gottlieb et al. (2001). The exact definition of is irrelevant to the present proof, hence not provided here. Injecting the explicit definition of the TASE operators given in Eq. 14, with coefficients given in Table 2, we obtain the following stability function

 σ(z)=⎛⎜⎝1+s−1∑k=0As,kzk+1(p−1∑i=0βp,i(2i−αz)−1)k+1⎞⎟⎠. (24)

Using the multinomial theorem, Eq. (24) can be rearranged as

 σ(z)=⎛⎝1+s−1∑k=0As,kzk+1∑ℓ1+ℓ2+⋯+ℓp=k+1(k+1ℓ1,ℓ2,…,ℓp)p−1∏i=0βℓip,i(2i−αz)−ℓi⎞⎠. (25)

It is possible to express the above stability function as where and are defined as

 Q(z)=p−1∏i=0(2i−αz)s, (26)

and

 P(z)=⎛⎝Q(z)+s−1∑k=0As,kzk+1∑l1+l2+⋯+lp=k+1(k+1ℓ1,ℓ2,…,ℓp)p−1∏i=0βℓip,i(2i−αz)s−ℓi⎞⎠. (27)

The key here is to note that hence is a polynomial and, furthermore, its roots are not in . Therefore is an analytic function in . As a polynomial, is also an analytic function in . We can then conclude that as the ratio of two analytic functions where the denominator that does not admit poles in , the stability function is proved to be an analytic function in for arbitrary order TASE and arbitrary order RK scheme.

(b) Analytically investigating whether the stability function defined by Eq. (25) satisfies the condition for purely imaginary is difficult. We have shown previously in Section 2.4 that a necessary condition is given by Eq. (11), which also applies to the case , with real . However, theorem 2.1 does not guarantee the stability of the method on the entire imaginary axis for intermediate values of the time-step . Large values of the time step are supposedly the worst case for stability, yet nothing guarantees that this holds mathematically. In order to prove that for any , the time integration of Eq. (4) remains stable as long as Eq. (11) is satisfied, and complete the present analysis, we use a computer program to plot using a dense sampling of the imaginary axis over very large values of . The results are reported in Fig. 4. The plots show consistent results with the analysis of the stability diagrams made in Section 2.6. In cases where the method is only nearly unconditionally stable, it is worth noting that the detrimental maximal value of the stability function is at most approximately .

Although no analytic proof is provided for the second part of corollary 3.4.1, the numerical proof performed in Fig. 4 can be repeated for another time-marching scheme as needed.

## 4 Numerical verifications

In this section, a set of stiff physical problems, including those with over-resolved spatial resolution, inhomogeneous boundary conditions, multiple stiffphysical process, nonlinear operators, and grid-induced stiffness, will be conducted to assess the stability and accuracy of the proposed TASE framework.

### 4.1 Prototype model examples

Consider the 1D diffusion of a scalar in the infinite domain described by the PDE

 (28)

where denotes the amplitude of the unsteady source term and is its time scale. Assuming the initial condition

 y(x,0)=1−cos(x), (29)

the exact solution is

 y(x,t)=1−cos(x)exp(−t)−Aτs(cos(t/τs)−1). (30)

Since the solution is -periodic, the numerical solution is sought on a spatial grid and periodic boundary conditions are used. The grid resolution , where is the total number of collocation points used to describe the solution, should be chosen small enough to resolve the smallest physical length scale of the problem. In this problem, it is the lengthscale of the initial condition, and is approximated as its wavenumber . The temporal resolution must be smaller than the physical time scales of the problem. There are two physical time scales in this problem, one that describes the fast diffusion of the initial spatial inhomogeneities, , and another that represents the slow uniform oscillation imposed by the source, . Depending on the grid resolution , the final time of interest , and the source amplitude , an alternative to traditional explicit schemes becomes motivated for one of the reasons described in Section 1. The objective of this section is to demonstrate the accuracy and stability properties of TASE operators in all of these scenarios, and remind how the choice of the time step remains relevant in these cases. Over-resolved, quasi-steady and steady cases are examined in Sections 4.1.1, 4.1.2 and 4.1.3, respectively. The semi-discrete version of Eq. (28) resembles that given in Eq. (3), where is the numerical approximation of the second-order derivative operator, but is supplemented with a time dependent source term on the right hand side that is independent from the diffusion term. This second term is not modified in the following, and further discussions and justifications about general treatment of source terms is provided in Section 4.2.

In all cases, the maximum time step for numerical stability is given by , where given in Table 1 depends on the chosen explicit time-marching scheme, and depends on the numerical approximation of the second-order derivative. For second-order central finite-difference, fourth-order central finite-difference, and Fourier differentiation, is respectively equal to , and Moin (2010). The baseline implicit method is the trapezoidal method.

#### 4.1.1 Scenario 1: the physical features are spatially over-resolved

This section treats the example of over-resolving the spatial scales of the problem, in a way that results in an unphysical constraint on the maximum allowed time step for numerical stability. The short-time solution of Eq. (28) at final time is is of interest, the source amplitude is set to 0 and degree of over-resolution is varied with the number of points from 6 to 60 to 600. The spatial and time derivatives are approximated using fourth-order central finite-difference and second-order RK, respectively.

Figure 5(a,b) show the solution profile at five times during the transient decay of the solution for two cases: the least and the most over-resolved case that respectively correspond to resolving the initial profile with or points, respectively, which correspond to and . In both cases, the time step for the explicit method is equal to the maximum time step for linear stability, with the exception of being rounded so that an integer number of time steps reaches exactly the final time of interest, i.e. , which in turn yields and for and , respectively. On the contrary, the TASE solution is obtained with a physics-based time step to resolve the smallest physical time scale of the problem, i.e. , irrespectively of the maximum time step of stability. This yield and for and , respectively. Profiles obtained with the different numerical methods are almost indistinguishable to plotting accuracy. The remarkable agreement is obtained while in the second case the explicit + TASE solver uses a time step larger than that used by the explicit solver by four orders of magnitude.

Figure 5(c) shows the relative error of the final solution profile in function of the time step normalized by the maximum time step allowed for linear stability, , where the different curves correspond to three different levels of over-resolution , and . It confirms the expected order of accuracy of the TASE method. With the addition of the TASE operator, the explicit method has preserved its order of accuracy and becomes unconditionally stable. The error plots shifts downward with increasing over-resolution, which shows that the more over-resolved the initial profile is, the more justified the use of a time step larger than the largest one allowed by stability analysis. The error of the TASE method is slightly superior to that of the explicit method without TASE due to the error introduced by the TASE operator, although it remains a small error. The TASE error is larger than that of the baseline implicit method used in this example. The error saturates for the least two over-resolved cases as it becomes dominated by the spatial differentiation error. For each numerical method, all curves for a given time-marching scheme collapse when plotted versus instead, owing to the fact that the only physics-based requirement to obtaining an accurate solution is the resolution of the physical time scale of the initial solution profile. The plot shows that the error approaches order unity when this ratio becomes of order unity. In order to attain a given error within the bounds of engineering interest, the required time step that can be used with the TASE method when the solution is over-resolved is unusable with traditional explicit methods, that would be unstable in that case. In the inset, the same plot is shown when the Fourier differentiation is used instead of fourth-order central finite difference. Since the spatial differentiation becomes exact in this case, the plot reflect time errors only and the the saturation at small time steps disappears.

In this example, over-resolution of the solution justifies the use of a large time step, that cannot be used with an explicit method, as long as it resolves the smallest physical time scale of interest. The more over-resolution, the more useful TASE becomes. In the limit of a normally resolved feature, TASE remains accurate. This is satisfied by construction of the TASE operator, and is a useful characteristic for multi-scale problems that are inhomogeneous and intermittent.

#### 4.1.2 Scenario 2: the user is interested in the steady-state solution of a stiff problem

Now consider the case where , i.e. , but contrary to before, the user is interested in the final steady-state solution, i.e. . Here we choose . The spatial and time derivatives are approximated using exact Fourier spectral differentiation and fourth-order RK, respectively.

Figure 6(a) show the time decay of the solution at . The time step for the explicit method is approximately equal to the maximum time step for linear stability, i.e. , which corresponds to . On the contrary, the TASE methods need not resolve this time scale but only the time scale of interest, which is here the final time . For illustration, Figure 6(a) shows the profile obtained with the TASE methods for and , which correspond to and , respectively. Profiles obtained with the different numerical methods are almost indistinguishable to plotting accuracy in the limit of small time step. However, the prediction of the solution at the final time of interest is predicted remarkably well by all schemes irrespective of the chosen time step, because in all cases the physics-based constraint is satisfied.

Figure 6(b) shows the relative error of the final solution profile in function of the time step normalized by the maximum time step allowed for linear stability, , and normalized by the final time of interest . It confirms the expected order of accuracy of the TASE method. Compared to the previous example, it shows the flexibility of TASE to handle different order. With the addition of the TASE operator, the explicit method has preserved its order of accuracy and become unconditionally stable. The plot shows that the error approaches order unity when the time step becomes of the order of the final time step of interest. In order to predict a quantity of interest within error bounds of engineering interest, the required time step that can be used with the TASE method is unusable with traditional explicit methods, that would be unstable in that case.

In this example, the separation of scale between the smallest physical time scale that corresponds to the initial decay of the solution and the final time scale of interest justifies the use of a time step much larger than the maximum time step predicted by linear stability analysis. The larger the separation of scale, the more useful TASE becomes. Here, the separation was by a factor fo 10 only, but in problems of practical interest, it can easily be much larger.

#### 4.1.3 Scenario 3: the stiff physics is quasi-steady and irrelevant to the long term transient solution

Now consider the case where , as in Section 4.1.2 but the source term is activated, and the user is not interested in the initial quasi-steady fast decay of the initial profile. Instead the objective is only to accurately capture the slow oscillations encoded by the source term, i.e. . The amplitude source is . The spatial and time derivatives are approximated using exact Fourier spectral differentiation and second-order RK, respectively.

Figure 7(a) show the time development of the solution at . The time step for the explicit method is equal to the maximum time step for linear stability, i.e. , which in turn yields . On the contrary, the TASE methods need not resolve this time scale but only the time scale of interest, which is here the time scale of the source . For illustration, Figure 7(a) shows the profile obtained with the TASE methods for , and , which correspond to , and , respectively. Profiles obtained with the different numerical methods are almost indistinguishable to plotting accuracy in the limit of small time step. However, the prediction of the solution after an initial fast transient decay is predicted remarkably well by all schemes irrespective of the chosen time step, because in all cases the physics-based constraint is satisfied.

Figure 7(b) shows the relative error of the final solution profile in function of the time step normalized by the maximum time step allowed for linear stability, , and normalized by the source time scale . It confirms the expected order of accuracy of the TASE method. With the addition of the TASE operator, the explicit method has preserved its order of accuracy and become unconditionally stable. The plot shows that the error approaches order unity when the time step becomes of the order of the source time scale. In order to predict a given quantity of interest within the error bounds of engineering interest, the required time step that can be used with the TASE method is unusable with traditional explicit methods, that would be unstable in that case. The effect of is also shown. With increasing , the order of convergence is preserved, but the error is increased. Note that in practice it is not a tunable parameter, as theoretically prescribed as described in Section 2.

In this example, the separation of scale between the smallest physical time scale that corresponds to the initial decay of the solution and the source time scale justifies the use of a time step much larger than the maximum time step predicted by linear stability analysis. The larger the separation of scale, the more useful TASE becomes. Here, the separation was by a factor of 100 only, but in problems of practical interest, it can easily be much larger.

### 4.2 Application to problems with inhomogeneous boundary conditions (source terms)

#### 4.2.1 Method

In this section, we treat a more general example where the boundary conditions are implemented as a source term in the semi-discrete equation. In this case, the semi-discrete system typically takes the form

 dYdt=LY+S, (31)

where is a vector that encodes for the physical boundary conditions, and that may depend on and the specific choice of numerical operators used at the boundaries. When the TASE methodology is used, the TASE operator simply needs to be applied to both terms simultaneously, yielding the modified semi-discrete equation

 dYdt=T(p)L(LY+S), (32)

contrary to the inaccurate one that would be obtained if the TASE operator was only applied to the linear operator, . The simultaneous application of the TASE operator to both terms is necessary since both encode for a single physical process and compete near the boundary. In particular, this approach guarantees that the correct quasi-steady or steady solution of the system is captured, as can be understood by substituting in Eq. (32), and solving for . The presence of a source term originating from the boundary condition discretization is a particular case of the more general case of competing operators, which is treated in Section 4.3.1.

#### 4.2.2 Example

In order to demonstrate the method for handling source term, and how the failure manifests if the boundary conditions implemented as a source term is not subject to the TASE operator, the example from Section 4.1.1 is repeated in the subdomain so that the exact boundary Dirichlet boundary conditions yields a non-zero source term in the semi-discrete equation. The example uses collocation points, and the spatial and time derivatives are approximated using fourth-order central differencing and second-order RK schemes, respectively.

Figure 8(a) shows a comparison of the instantaneous solution profile obtained at an intermediate time for different numerical methodologies. The time step for the explicit method is equal to the maximum time step for linear stability, i.e. . On the contrary, the TASE methods use a physics-based time step to resolve the smallest physical time scale of the problem, i.e. , irrespectively of the maximum time step of stability. This yield . The figures compares the profiles obtained with the correct implementation of the TASE methodology, where both linear operator and source term are modified, to the incorrect one where only the linear operator is modified, as described in Sec. 4.2.1, with . The figure shows that similar conclusions as in Sec. 4.1.1 applies when the TASE methodology is correctly implemented. On the contrary, large errors arise near the boundaries when the boundary source term is not modified, and subsequently propagate in the interior of the domain.

The observations described above are confirmed in Fig. 8(b) that shows the relative error of the final solution profile in function of the time step normalized by the maximum time step allowed for linear stability, . Besides conclusions already drawn in previous part of the paper, it is worth noting that the application of TASE operators with wrong handling of boundary source terms preserves its order of accuracy but much larger errors due to the inconsistent numerical treatment of the diffusion physics in the interior and at the edges of the domain.

### 4.3 Application to problems with multiple stiff physical processes (operator splitting)

#### 4.3.1 Method

In this section, we treat a more general example where the problem involves multiple stiff linear operators and

 dYdt=L1Y+L2Y. (33)

The TASE operators can be used to bypass the linear stability restrictions of each operators. Each operator can be modified combined

 dYdt=T(p)L1+L2(L1Y+L2Y), (34)

or separately

 dYdt=T(p)L1L1Y+T(p)L2L2Y. (35)

The choice of one strategy versus the other is problem-dependent, and depends on whether the underlying quasi-steady process to be bypassed by the temporal resolution involves competition of both operators. When the operators compete, i.e. they are both simultaneously active and significant in a (quasi-)steady process at the same location, a combined TASE operator must be constructed and applied to both terms, as in Eq. (34). Alternatively, when they do not compete, it is sufficient to build distinct TASE operators and and apply them to the respective terms in the equation, as in Eq. (35). This physics-based constraint for splitting or grouping the application of TASE to the right-hand side operators is reminiscent of the constraint that underpins the use operator splitting Marchuk (1968); Strang (1968). In general, competing operators should not be split but computed simultaneously in this latter methodology Strang (1968). In the absence of knowledge about the competition among physical operators, the user may safely adopt the strategy described in Eq. (34). It is more general but since its implementation can be more cumbersome and more costly than the alternative one, it is desirable to use Eq. (35) when prior knowledge about the non-competing stiff operators exists.

#### 4.3.2 Example

 ∂y1∂t+U∂y1∂x=D∂2y1∂x2+K(y2−y1), (36)
 ∂y2∂t+U∂y2∂x=D∂2y2∂x2−K(y2−y1), (37)

in the bounded domain . The semi-discrete equation becomes

 dYdt=(Lc+Ld+Lr)Y+(Sc+Sd), (38)

where is the concatenated vector of unknowns, , , and denote the discretized operators that correspond to the convection, diffusion, and reaction terms, respectively. The source terms are denoted and , respectively. The numerical algorithm employs a non-uniform mesh stretched near the right-boundary in order to capture the fine feature that develops over time there, and a finite-volume discretization of the equation that employs second-order spatial derivative operators, with control volumes. The time-marching scheme is third-order RK. The problem parameters , and are chosen in such a way that the reaction time scale is much faster that the convection and diffusion one, i.e. , and . Assuming that the physical length scale is , this renders the following correspond time scales: ,