# Exponential Time Differencing for the Tracer Equations Appearing in Primitive Equation Ocean Models

The tracer equations are part of the primitive equations used in ocean modeling and describe the transport of tracers, such as temperature, salinity or chemicals, in the ocean. Depending on the number of tracers considered, several equations may be added to and coupled to the dynamics system. In many relevant situations, the time-step requirements of explicit methods imposed by the transport and mixing in the vertical direction are more restrictive than those for the horizontal, and this may cause the need to use very small time steps if a fully explicit method is employed. To overcome this issue, we propose an exponential time differencing (ETD) solver where the vertical terms (transport and diffusion) are treated with a matrix exponential, whereas the horizontal terms are dealt with in an explicit way. We investigate numerically the computational speed-ups that can be obtained over other semi-implicit methods, and we analyze the advantages of the method in the case of multiple tracers.

There are no comments yet.

## Authors

• 2 publications
• 2 publications
• 13 publications
10/25/2021

### Partially Explicit Time Discretization for Nonlinear Time Fractional Diffusion Equations

Nonlinear time fractional partial differential equations are widely used...
03/24/2021

### Hall Effects on Casson Fluid Flow along a Vertical Plate

The Hall effects on Casson fluid flow along a vertical plate has been in...
08/12/2019

### An energy consistent discretization of the nonhydrostatic equations in primitive variables

We derive a formulation of the nonhydrostatic equations in spherical geo...
05/27/2021

### High-Order Multirate Explicit Time-Stepping Schemes for the Baroclinic-Barotropic Split Dynamics in Primitive Equations

In order to treat the multiple time scales of ocean dynamics in an effic...
07/25/2017

### Kinetic Simulation of Collisional Magnetized Plasmas with Semi-Implicit Time Integration

Plasmas with varying collisionalities occur in many applications, such a...
04/01/2022

### Nonlocal transport equations in multiscale media. Modeling, dememorization, and discretizations

In this paper, we consider a class of convection-diffusion equations wit...
04/29/2020

### An energetically balanced, quasi-Newton integrator for non-hydrostatic vertical atmospheric dynamics

An energetically balanced, implicit integrator for non-hydrostatic verti...
##### 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

The primitive equations are the constitutive system of equations in ocean modeling. This system is composed of a momentum equation, an equation for the ocean thickness, equations for the transport of tracers, as well as an equation of state. The momentum and thickness equations describe the dynamics, i.e. the change in time of the velocity and thickness of the water. The tracer equation describes the transport of tracers, such as temperature, salinity or chemicals. Depending on the number of tracers considered, several equations may be added to and coupled to the dynamics system. The coupling between the tracers and the dynamics depends on the nature of tracers. Temperature and salinity impact the density of the layers of ocean water, influencing the dynamics, and for this reason they are called active tracers. Hence, between the dynamics and active tracers there is a two-way coupling. However, for most tracers the coupling is only one-way, with such tracers being called passive [1].

To this day, the numerical solution of the primitive equations remains a challenging task. One of the main issues is the presence of multiple time-scales, where different processes (e.g., external and internal gravity waves, eddies, biochemical reactions) take different times to be completed. Since the primitive equations arise from hyperbolic conservation laws ([2, 3]), explicit time integrators and Runge-Kutta schemes would hypothetically be good choices for solving it. However, these schemes are not capable of efficiently handling multiple time-scales, because the time-steps restrictions are too severe, resulting in a significant degradation of performance. To better handle multiple time-scales, methods have been developed with the help of mathematical and algorithmic techniques such as splitting strategies and semi-implicit approaches. Due to the different properties of the dynamics and the tracer equations, schemes have been created to separately deal with the two subsystems ([4, 5]). Methods for the dynamics need to take into account the different wave speeds of the model, with the goal of having a model for which time step sizes are governed by the slow wave speeds or the speed of advection and not by the fast wave speeds as is the case for standard explicit schemes. Examples of time-stepping schemes for the dynamics include implicit ([6]) and split-explicit ([7, 4]) methods. More recently, exponential time differencing (ETD) methods, also known as exponential integrators, have gained attention in the ocean modeling community due to their stability properties that allow time steps considerably larger than those dictated by the CFL condition. In [8], an ETD scheme has been developed for the rotating shallow water equations with multiple horizontal layers, which correspond to a vertical discretization of the primitive equations dynamics in an isopycnal vertical coordinate system. The main idea behind exponential integrators is a splitting of the right-hand side term of an equation into a linear part and a remainder, i.e.

 ∂tθ=F(θ)=Aθ+R(θ),

with an appropriate choice of the linear operator . For a review of exponential integrators we refer to [9].

In this work, we devise an ETD method for the tracer equation. A tracer is supposed to satisfy a conventional advection-diffusion equation of the form,

 ∂tθ+∇x⋅(uθ)+Dxθ=q(θ) (1)

where is the tracer studied, is the velocity of water, which is usually split in to the horizontal velocity and the vertical velocity , is a diffusion term, and represents interior sources or sinks. The equation (1) has to be solved on a three-dimensional domain, split into a two-dimensional horizontal and a vertical coordinate , and . We note that, the tracer equation can usually also be written as

 ∂tθ+∇x⋅(uθ)+Dxθ+∂z(wθ)−∂z(κz∂zθ)=q(θ), (2)

with the horizontal diffusion , and the vertical diffusion coefficient . To avoid problems with shocks, the velocity field is usually chosen to be divergence free and tangential to the boundary.

When dealing with the tracer equation, a key point is appropriately including vertical mixing. Tracer vertical mixing usually occurs on small time-scales and can be induced by density differences and/or by turbulent motions. For explicit time stepping schemes, the time step requirement is usually set by the horizontal advective CFL condition, hence very small time steps may need to be used with explicit time stepping methods to include realistic vertical mixing of tracers. To avoid this issue, in the time-stepping schemes used by popular ocean models, the vertical diffusion term is treated implicitly. In POP [5], the vertical tracer diffusion term is treated with an implicit Euler algorithm, whereas the remaining terms of the equation are treated with a leapfrog algorithm. In MPAS-Ocean [4], tracer equations are stepped forward with the mid-time velocity values and this process is repeated in a predictor-corrector way. Implicit vertical mixing of tracers completes each time-step, where, as in POP, the vertical tracer diffusion term is treated with an implicit Euler algorithm.

Another important phenomenon to deal with, besides vertical mixing, is when an inflow of cold water near the coast leads to cold water on top, scenario that creates density and pressure variations and downward motion. In an isopycnal configuration, there is no vertical transport, however, if a mostly Eulerian coordinate system is employed (z-star, z-level), rapid variations in the pressure/density may induce a fast vertical flow. This phenomenon is amplified by the usage of fine meshes in the vertical, with much smaller vertical than horizontal spacing, for example [4] employs km horizontal, m vertical. As we said, many ocean models treat the vertical advection explicitly and so, when a fast transport of water among different layers occurs, the model may be unable to appropriately capture this behavior. Consequently, instabilities in the simulation may occur, causing the need to decrease the time-step. To resolve this issue, we propose an ETD solver where all the vertical terms, i.e. vertical advection and diffusion, are treated with a matrix exponential, whereas the horizontal are dealt with in an explicit way. This means that we are splitting the linear operator into two parts: that accounts for all vertical terms, and that contains all horizontal terms. is then incorporated in the remainder , so the actual linear operator we are working with is . This operator splitting has two advantages. First, by treating exponentially terms related to fast time-scales, bigger time steps can be taken, and so computational speed-ups that can be obtained over other explicit methods. Compared to semi-implicit methods, we expect higher accuracy due to an exact treatment of the fast scales. Second, by including only the operator in the exponential, not the whole operator , the computational cost is reduced, effectively lowering the total computational cost and time of the whole method.

Another important challenge to face when dealing with tracer equations, and in general with primitive equations, is in the numbers of tracers considered. The amount of tracers in an ocean simulation is usually around 40, but may increase up to 70, causing a significant computational load. Hence, efficiently solving multiple tracers equations is an important task in an ocean model. When an ETD scheme is used, efficiency corresponds to a low-cost evaluation of -functions. For a given matrix , assembling is generally prohibitive in the large scale context; iterative methods are used instead, i.e. Krylov subspace algorithms ([10, 11]). Scaling and squaring is used usually for dense matrices ([12, 13, 14, 15, 16]) However, it has also been used in the context of multiwavelet Galerkin methods, where a certain level of sparsity can be maintained throughout the iterations of the method and the approximation of ([17]). In this work, together with the usage of a Krylov method, we pursue a second approach based on scaling and squaring relations, focusing on the fact that there is no communication in ocean due to the vertical exponential, since the domain-decomposition is horizontal, and the exponential is vertical. The proposed approximation scheme is based on polynomials of moderate degree that result in a consistent and stable approximation of with low bandwidth. To efficiently solve multiple tracer equations, we intend to preassemble instead of computing for every right hand side . Additionally, since the current ocean models use only 40-100 layers, even a full storage of exponential is feasible.

Finally, our solver is compared with existing ocean models, to make sure that it is able to reproduce similar results under the same physical conditions. To do so, the whole primitive equations are solved and two tests from [18] are performed. The tracer equations is coupled with the dynamics, and to solve this latter system a second ETD solver (Exponential Rosenbrock Euler) is used. The results obtained with our ETD solvers are compared with those obtained with three other codes: MPAS-Ocean, MITgcm and MOM.

The paper is organized as follows. Section 2 describes the discretization of the tracer equation both in the vertical and in the horizontal. Section 3 focuses on exponential integrators and their implementation. For the computation of -functions, a restarted Krylov subspace method is described, and a scheme based on scaling and squaring relations is proposed. For the latter scheme, an error analysis is provided. In section 4, the proposed ETD solver for the tracer equation is presented, and the properties of the operator splitting are described. The numerical tests are shown and discussed in section 5, while the conclusions follow in section 7.

## 2 Discretization of the Tracer Equation

In this section, the discretization of the tracer equation (1) in the vertical and in the horizontal is described. In most ocean models, an hydrostatic condition is assumed, leading to describing the primitive equation by the incompressible Boussinesq equations in hydrostatic balance. Following this assumption, the tracer equation in continuous form can be re-written as

 ∂(˜ρT)∂t=−∇⋅(˜ρuT)−∂(˜ρTw)∂z+DTx+DTz+FT, (3)

where is the tracer, is the horizontal velocity, is the vertical velocity and is the pseudo-density. The variable represents the vertical coordinate and it defined positive upward. and indicate the horizontal and vertical diffusion term, respectively. These terms are defined as

 DTx =∇⋅(hκx∇T), (4) DTz =h∂∂z(κz∂T∂z), (5)

where and are the horizontal and vertical diffusion, respectively. Without loss of generality, we assume that no forcing term is present, i.e. . This assumption is equivalent to consider that no external factors have an influence on the tracer behavior.

To discretize equation (3), we employ z-level coordinates in the vertical [19] and a finite-volume method using a C-grid staggering in the horizontal [20]. For more details about the discretization of all primitive equations, please refer to [4].

### 2.1 Tracer Equation with vertical discretization

As in MPAS-Ocean [4, 19], the vertical coordinate we use is Arbitrary Lagrangian-Eulerian (ALE). With ALE, several coordinate systems can be specified depending on the application. Common choices for the vertical coordinates include z-level, where all layers have a fixed thickness except for the top layer, z-star, where all layer thicknesses vary in proportion to the sea surface height, and isopycnal, where there is no vertical transport between layers.

The tracer equation with vertical discretization is

 ∂(hkTk)∂t=−∇⋅(hkukTk)−¯¯¯¯¯Tkwk+¯¯¯¯¯¯¯¯¯¯Tk+1wk+1+[DTx]k+[DTz]k, (6)

where indicates the vertical layer, is the top layer and increases downward up to ; is the mean elevation of the free surface, and the coordinate is positive upward. The variable indicates the transport of fluid from layer to , i.e. across the top interface of layer . The pseudo-density has been replaced by the , which is the layer-thickness. The operator , on a generic variable , is the vertical average between the layer and the above layer , i.e.,

 ¯¯¯¯¯¯ψk=ψk−1+ψk2. (7)

Finally, and indicate the discretized horizontal and vertical tracer diffusion terms, respectively, and are defined as

 [DTx]k=∇⋅(hkκx∇Tk),[DTz]k=hkδzk––––(κz¯¯¯¯¯¯¯δzk(Tk)). (8)

The discrete operators and , on a generic variable , are defined as

 ¯¯¯¯¯¯¯δzk(ψk)=ψk−1−ψk¯¯¯¯¯hk, (9) δzk––––(ψk)=ψk−ψk+1hk. (10)

The vertical tracer diffusion term can actually be rewritten without introducing the operator as .

The choice of the vertical coordinate system is enforced in the computation of the vertical transport. can be found by solving the thickness equation for . The thickness equation discretized in the vertical has the form

 ∂hk∂t+∇⋅(hkuk)+wk−wk+1=0. (11)

For a Boussinesq fluid, this equation represents the continuity equation for the pseudo-density . To obtain from (11), all variables at the previous time step must be known, in particular the time derivative of at layer , , must be known. For this purpose, a new quantity, named , is introduced. represents the desired thickness for the new time, and it is used to compute using a first-order finite difference approximation. In this way, can be found as

 wk=wk+1−∇⋅(hkuk)−hALEk−hkΔt. (12)

For isopycnal simulations, the vertical transport is set to zero in (12). For other coordinate systems, like z-level and z-star, the way is computed determines the type of coordinates chosen. For example, for z-level vertical coordinates

 hALE1 =hrest1+ζ, hALEk =hrestk, for k>1,

where is the layer thickness when the ocean is at rest, and is the sea surface height defined as . For z-level coordinates (and for z-type coordinates in general) the resting thickness is considered constant in each horizontal layer, but for other coordinate systems, like sigma coordinates, varies horizontally in proportion to the column’s total depth. The simulations presented in section 5 use z-level vertical coordinates. Please refer to [19] for more details about the computation of for other coordinate systems.

### 2.2 Tracer Equation with horizontal discretization

The horizontal discretization is a C-grid, finite-volume method applied to a spherical centroidal Voronoi tessellation (SCVT) mesh. Height, tracers, pressure and kinetic energy are defined at centers of the convex polygons, and the velocity is located at cell edges. Vorticity (curl of velocity) is defined at cell vertices. In the following, the subscripts and indicate the discretized variables through cell centers and edges, respectively. Since we are focusing on the discretization of the tracer equation only, we will not work with variables and operators defined at cell vertices.

The tracer equation with horizontal discretization is

 ∂(hk,iTk,i)∂t= −[∇⋅(ˆ(hk,:)euk,:ˆ(Tk,:)e)]i−¯¯¯¯¯¯¯¯Tk,iwk,i+¯¯¯¯¯¯¯¯¯¯¯¯¯Tk+1,iwk+1,i +[DTx]k,i+[DTz]k,i, (13) k,i =[∇⋅(ˆ(hk,:)eκx[∇Tk,:]e)]i,[DTz]k,i=hk,iδzk––––(κz¯¯¯¯¯¯¯δzk(Tk,i)). (14)

Each variable now has two sub-scripted indices, the first indicating the vertical layer, and the second indicating its position on the horizontal grid, namely either or . Colons in subscripts may be places as second index to indicate that multiple edges or cell centers are used in computing the horizontal operator. For a generic variable , the symbol represents the averaging of the variable from two adjacent centers to the corresponding edge. We would like to point out that the vertical transport through the sea surface and at the bottom surface is zero, i.e. and . Moreover, we consider on all boundary edges.

For a generic vector field

and variable , the discrete horizontal operators and are defined as

 [∇⋅Yk,:]i =1Ai∑e∈E(i)ne,iYk,ele, (15) e =1de∑i∈C(e)−ne,iψk,i. (16)

indicates the Voronoi cell area, is the distance between cell centers, is edge length and represents the sign of the vector at edge with respect to cell . The sets are the edges about cell , and the sets are the cells neighboring edge . Thus, the divergence moves from edges to cell-centered quantity, while the gradient moves from cell centers to edges.

## 3 Exponential Time Integration

This section describes exponential time differencing methods, that later will be used to solve the tracer equation. ETD methods have already been employed to solve the single layer ([21, 22, 23]) and multi-layer ([8]) shallow water equations.

### 3.1 Exponential Integrators

Let

be a system of partial differential equations (PDEs), where

denotes the vector of the solution variables for , and is the right-hand-side term. The interval refers to one time step. The main idea behind exponential integrators is a splitting of the right-hand-side term into a linear part and a remainder, i.e.

 ∂tT=F(T)=AnT+R(T), (17)

where represents a linear operator, and denotes the remainder, which in general is nonlinear. Applying the variation of constants formula to equation (17), the solution at time , i.e. , is obtained as

 Tn+1=exp(ΔtAn)Tn+∫Δt0exp((Δt−τ)An)R(T(tn+τ))dτ. (18)

At this point, to build a concrete exponential integrator, an approximation of must be considered. By substituting with its Taylor expansion truncated at (), namely

the solution can be approximated by

 Tn+1≈exp(ΔtAn)Tn+s∑k=11(k−1)![∫Δt0exp((Δt−τ)An)τk−1dτ]dk−1R(T(tn+τ))dτk−1∣∣ ∣∣τ=0.

The above expression can be rewritten as

 Tn+1≈exp(ΔtAn)Tn+s∑k=1Δtkφk(ΔtAn)dk−1R(T(tn+τ))dτk−1∣∣ ∣∣τ=0, (19)

by introducing the, so called, -functions defined by

 φk(ΔtAn) =1Δtk(k−1)!∫Δt0exp((Δt−τ)An)τk−1dτ (20) =1(k−1)!∫10exp((1−σ)ΔtAn)σk−1dσ,k=1,2,…,s. (21)

By performing the change of variable , (21) can be obtained from (20). For , we have that .

The parameter in (19) indicates the number of stages of the method. By taking , the exponential Euler method is obtained as

 Tn+1≈exp(ΔtAn)Tn+Δtφ1(ΔtAn)R(Tn)=Tn+Δtφ1(ΔtAn)F(Tn), (22)

where , with

indicating the identity matrix. For a generic linear operator

, this method is first-order accurate, but if is the Jacobian matrix of the system evaluated at , namely , then the method becomes second-order accurate [9, 23]. In this case, the scheme (22) is called Exponential Rosenbrock Euler. By taking , a second-order single-step method with two stages is obtained as

 T(1ststage)n=Tn+Δtφ1(ΔtAn)F(Tn), (23) Tn+1=T(1ststage)n+Δtφ2(ΔtAn)(R(T(1ststage)n)−R(Tn)), (24)

where the first derivative of is approximated with a first-order finite-difference approximation

 dR(T(tn+τ))dτ∣∣∣τ=0≈R(T(1ststage)n)−R(Tn)Δt,

and . This scheme is known as ETD2-RK ([24, 9]).

ETD2-RK fulfills the stiff order conditions of order two described in [24], section 5.1. However, by relaxing some of these conditions, different schemes can be obtained that are more convenient for computations. For instance, the following two-stage predictor-corrector scheme makes use only of the function.

 T(1ststage)n=Tn+Δtφ1(ΔtAn)F(Tn), (25) Tn+1=T(1ststage)n+12Δtφ1(ΔtAn)(R(T(1ststage)n)−R(Tn)). (26)

This scheme fulfills the nonstiff order conditions up to order two, and the stiff order conditions up to order one. From an implementational point of view, it is simpler than ETD2-RK, since does not need to be assembled, and a further advantage in terms of computational time can be obtained if the matrix can be precomputed and stored efficiently.

### 3.2 Computation of the φ-functions

The computation of the -functions is of major relevance for ETD methods. Let us assume that is in , where and is the index of the -function. Fist of all, let us consider the special case , for which we have . The most popular method for computing the matrix exponential is the scaling and squaring algorithm [25, 15]. In [14], the scaling and squaring algorithm used by the MATLAB function expm is described, and a variation of this method that alleviates overscaling problems is presented in [12]. The Expokit package [26] uses a scaling and squaring method as well for computing matrix exponentials. An alternative to scaling and squaring methods is given by Krylov subspace projections [27, 11]. Krylov schemes can be used for evaluating any -function with index . With Krylov schemes, matrix-vector products of the form are computed, without the actual construction of the matrix . To overcome some memory issues associated with standard Krylov methods, restarted schemes have been developed [10]. The C/C++ library SLEPc [28] uses the restarted algorithm presented in [10].

In the following, a brief description of the algorithm presented in [10] is given, since this method is used in the numerical results section, and scaling and squaring algorithm for indexes

is presented, for which an error estimate is derived.

#### 3.2.1 Krylov subspace approximation

The Krylov subspace approximation has been found very efficient to compute matrix-vector products due to the optimality of the matrix polynomials produced by Krylov methods. The idea behind a Krylov subspace approach is to approximate the vector , which lives in , in a smaller space of dimension . The Krylov approximation of is based on an Arnoldi decomposition of

 AVm=Vm+1˜Hm=VmHm+ηm+1,mvm+1eTm, (27)

where is a matrix whose columns form an orthogonal basis for the Krylov subspace of dimension

 Km(A,b)=span{b,Ab,…,Am−1b},

is an upper Hessenberg matrix, and denotes the th vector of the standard basis of . The matrices and are such that

 Hm=VTmAVm, (28)

therefore can be seen as the projection of the action of to the Krylov subspace .

Standard Krylov methods requires the storage of , i.e. the storage of vectors of size , which may be costly for moderate to large values of . Therefore, a standard Krylov subspace method may be impractical because of the storage requirements associated with . To overcome this issue, the Arnoldi approximation could be modified in a way that allows the construction of successively better approximations of based on a sequence of Krylov spaces of small dimension. Methods based on this approach are called restarted Krylov subspace methods, and they use an Arnoldi-like decomposition where a sequence of ascending (not necessarily orthonormal) basis vectors is introduced. We now focus on the restarted Krylov subspace algorithm presented in [10], which can be summarized as follows. An Arnoldi-like decomposition of is constructed

 AˆVp=ˆVpˆHp+np+1vpm+1eTpm, (29)

where , , and
, . Since (29) is an Arnoldi-like decomposition, the columns of are only blockwise orthonormal. The matrices , and the scalars are obtained from proper Arnoldi decompositions. Setting

 φk(ˆHp)=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣φ1,1kφ2,1kφ2,2k⋮⋮⋱φp,1kφp,2k…φp,pk⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦,where φj,jk=φk(Hj),j=1,2,…,p,

the approximation of after restart cycles is given by

 φpk=ˆVpφk(ˆHp)e1=[V1V2⋯Vp]φk(ˆHp)e1=p∑j=1Vjφj,1ke1=φp−1k+Vpφp,1ke1. (30)

Therefore, the approximation is obtained from the previous approximation plus a correction term. Only has to be stored from the previous cycle of the algorithm, and the matrix (together with ) can be discarded after computing . An efficient implementation of this algorithm can be found in [29], where it is also shown how to stably compute the coefficient vector .

#### 3.2.2 Scaling and Squaring

We now present a scaling and squaring method for the computation of . We are going to base the scaling and squaring method on the recursive relations for the -functions; cf. [30, 31, 32]. They are given as (see [32, Lemma 3]):

 2kφk(2A)=exp(A)φk(A)+k−1∑j=0φk−j(A)j! (31) = ⎧⎪⎨⎪⎩φk/2(A)2+2∑k/2j=01j!φk−j(A)for k even,φ(k−1)/2(A)φ(k+1)/2(A)+2∑(k−1)/2j=0φk−j(A)j!+φ((k+1)/2)(A)((k+1)/2)!for k odd. (32)

Note that in both expressions the matrix function of the matrix is expressed as a product of two matrix functions of plus some correction terms. The first identity is simpler, but the second one reduces the number of correction terms, which is slightly more efficient for computations; however [32] mentions that the first form is more stable in numerical experiments.

The first formula follows from the definition (20) using , and a splitting of the integral into the intervals and

 2kφk(2A) =∫10exp((2−τ)A)τk−1(k−1)!dτ+∫21exp((2−τ)A)τk−1(k−1)!dτ =exp(A)∫10exp((1−τ)A)τk−1(k−1)!+∫10exp((1−τ)A)(1+τ)k−1(k−1)!dτ,

where the second integral was shifted to . Now, using the binomial identity in the second term and the definition of and shows (31). Now, to derive the second identity, we use in a first step

 exp(A)φk(A)=(I+φ1(A)A)φk(A)=φk(A)+φ1(A)(Aφk(A))=φk(A)+φ1(A)(φk−1(A)−I/(k−1)!)=φ1(A)φk−1(A)+φk(A)−φ1(A)/(k−1)!.

Proceeding iteratively times with the first term, we obtain

 exp(A)φk(A)=φm(A)φk−m(A)+m∑j=1φk−j(A)j!−k−1∑j=k−m−1φk−j(A)j!.

using this identity for and respectively, we obtain (32).

In order to evaluate the matrix functions, we are going exploit these recursive relations to reduce the computation to matrix functions of a scaled matrix . Then, we use a polynomial approximation for this matrix. In order to ensure accuracy and stability, we do this in the following way: First, let

 p00(z)=Tr(z)+zr+1q(z)=exp(z)+O(|z|r+1)

be a polynomial of degree that approximates the exponential function up to order . Here, is the Taylor approximation to , and is a remainder. Then, for all we define the consistent polynomial approximations to the -functions as

 p0k(z)=z−k(pM0(z)−Tk(z))=r−k∑j=0zj(j−k)!+zr+1−kq(z) (33)

Now, we define the higher order recursive approximations for to using (31) as

 pMk(z)=2−k⎛⎝pM−10(z/2)pM−1k(z/2)+k−1∑j=0pM−1k−j(z/2)j!⎞⎠. (34)

This definition ensures that the resulting functions have similar properties as the original functions.

###### Proposition 1.

For any and there holds

 pMk(z) =z−1(pMk−1(z)−1/(k−1)!), (35) pMk(z) =z−k(pM0(z)−Tk−1(z)). (36)
###### Proof.

For this holds according to definition. For higher we follow an induction argument.

First, we show (35). We use the recursive definition (34) to obtain

 zpMk(z)=12k−1⎛⎝pM−10(z/2)(z/2)pM−1k(z/2)+k−1∑j=0(z/2)pM−1k−j(z/2)j!⎞⎠

and use the induction hypothesis to obtain

 2k−1zpMk(z) =pM−10(z/2)(pM−1k−1(z/2)−1(k−1)!)+k−1∑j=0pM−1k−j−1(z/2)−1/(k−j−1)!j! =pM−10(z/2)pM−1k−1(z/2)+k−2∑j=0pM−1k−1−j(z/2)j!−k−1∑j=01j!(k−j−1)! =2k−1pMk−1(z)−2k−1/(k−1)!,

using again the recursive definition (34) for and the well-known summation formula of binomial coefficients. Dividing by yields (35). Concerning (36), we note that is suffices to repeatedly apply (35) for , , …, . ∎

###### Remark 1.

The above result also implies that the recursion from (32) is equivalent for the construction of , since the equality (35) can be used to convert between both versions.

###### Remark 2.

In (45)-(46), which is the ETD method we actually use to later solve the tracer equation, the only -function needed is . For this specific case, the recursion formula (31) looks like

 φ1(A)=12(exp(12A)+I)φ1(12A).

and

 pM1(z)=z−1(pM0(z)−1).

Using the above relations, the following algorithm computes from .

Finally, we provide an error estimate for the approximation . First, we consider the polynomial approximation on a subset of the complex plane. To prepare for the general case, we let be some compact subset of the negative complex plane shifted by , and assume that the underlying polynomial fulfills the stability assumption

 |p00(z)|≤exp(τρ0)for all z∈τΣ, where 0≤τ≤1. (37)

Moreover, since is compact, there exists such that

 |p00(z)−exp(z)|≤cq|z|r+1for % all z∈Σ. (38)

This will be the basis of the error estimates.

###### Remark 3.

There are many situations where this assumption is fulfilled. In the case that is the Taylor polynomial and , can be chosen as the intersection of the negative half-plane and the well-known stability region of a Runge-Kutta scheme of order . In this case, due to the relation and we have . For an overview over known results on polynomials with optimal stability properties for various forms of and a computational approach to determine them, we refer to Ketcheson [33].

We first investigate the case of the matrix exponential .

###### Proposition 2.

Let the polynomial fulfill (37). Then for all we have the stability and approximation properties:

 |pM0(z)| ≤exp(τρ0), (39) |pM0(z)−exp(z)| ≤cqexp(τρ0)|z|r+12−Mr, (40)

which hold for all , for the enlarged stability radii .

###### Proof.

The first statement follows in a straightforward way from (37) since . For the second statement, we use an induction argument, noting that the case follows directly from (38). For , we let be arbitrary we use (34) to obtain

 ∣∣exp(z)−pM0(z)∣∣ =∣∣(exp(z/2)+pM−10(z/2))(exp(z/2)−pM−10(z/2))∣∣ ≤2exp(τρ0/2)cqexp(τρ0/2)|z/2|r+12−(M−1)r