AN A POSTERIORI—A PRIORI ANALYSIS OF MULTISCALE OPERATOR SPLITTING D. ESTEP∗ , V. GINTING† , D. ROPP

‡,

J.N. SHADID§ , AND S. TAVENER¶

Abstract. In this paper, we analyze a multiscale operator splitting method for solving systems of ordinary differential equations such as those that result upon space discretization of a reactiondiffusion equation. Our goal is to analyze and accurately estimate the error of the numerical solution, including the effects of any instabilities that can result from multiscale operator splitting. We present both an a priori error analysis and a new type of hybrid a priori – a posteriori error analysis for an operator splitting discontinuous Galerkin finite element method. Both analyses clearly distinguish between the effects of the operator splitting and the discretization of each component of the decomposed problem. The hybrid analysis has the form of a computable a posteriori leading order expression and a provably-higher order a priori expression. The hybrid analysis takes into account the fact that the adjoint problems for the original problem and a multiscale operator splitting discretization differ in significant ways. In particular, this provides the means to monitor global instabilities that can arise from operator splitting. Key words. a posteriori error analysis, adjoint problem, discontinuous Galerkin method, generalized Green’s function, goal oriented error estimates, multiscale method, operator decomposition, operator splitting, reaction-diffusion equations, residual AMS subject classifications. 65N15, 65N30, 65N50

1. Introduction. Operator decomposition is perhaps the most widely used technique for solving multiscale, multiphysics problems. The general approach is to decompose a model into components involving simpler physics over a relatively limited range of scales, and then to seek the solution of the entire system by using numerical solutions of the individual components. This approach has many appealing aspects. For example, it capitalizes on the significant progress that has been made on the stable, accurate and efficient solution of a broad spectrum of single physics problems. It provides an avenue to use highly evolved legacy codes to tackle multiphysics problems. It also provides a natural way to tackle problems encompassing multiple time and length scales. The classic example of operator decomposition is operator splitting for reactiondiffusion equations. The generic picture is a relatively fast, destabilizing reaction component interacting with a relatively slow, stabilizing diffusion component. Accuracy considerations dictate the use of relatively small steps to integrate the reaction component. On the other hand, stability considerations over moderate to long time intervals suggests the use of implicit, dissipative numerical methods for integrating ∗ Department of Mathematics and Department of Statistics, Colorado State University, Fort Collins, CO 80523 ([email protected]). D. Estep’s work is supported in part by the Department of Energy (DE-FG02-04ER25620, DE-FG02-05ER25699), the National Aeronautics and Space Administration (NNG04GH63G), the National Science Foundation (DMS-0107832, DGE0221595003, MSPA-CSE-0434354) and the Sandia Corporation (PO299784) † Department of Mathematics, Colorado State University, Fort Collins, CO 80523 ([email protected]). V. Ginting’s work is supported in part by the Department of Energy (DE-FG02-04ER25620) ‡ SAIC, 4001 Fairfax Dr. Arlington, VA 22203 ([email protected]) § Computational Sciences R&D Group, Sandia National Laboratories, PO Box 5800, MS0316, Albuquerque, NM 87185 ([email protected]) ¶ Department of Mathematics, Colorado State University, Fort Collins, CO 80523 ([email protected]). S. Tavener’s work is supported in part by the Department of Energy (DE-FG02-04ER25620)

1

2

D. ESTEP, V. GINTING, D. ROPP, J.N. SHADID, AND S. TAVENER

diffusion problems. Such methods are expensive to use per step, but relatively large steps can be used on a purely dissipative problem. If the reaction and diffusion components are integrated together, then the small steps required for accurate resolution of the reaction lead to an expensive computation. In a multiscale operator splitting approach, the reaction and diffusion components are integrated independently inside each time interval of a discretization of time and “synchronized” in some fashion only at the nodes of the interval. The reaction component is often integrated by using significantly smaller sub-steps (e.g 10−5 smaller is not uncommon) than those used to integrate the diffusion component, which can lead to a tremendous computational savings. However, operator decomposition presents an entirely new set of accuracy and stability issues, some obvious, some subtle, and all difficult to correct. In the case of operator splitting, the instantaneous interaction between reaction and diffusion is discretized and this has a strong effect on accuracy and stability, even if each component is solved exactly. For example, consider the instability observed in the Brusselator problem [28, 27].

2

10

1.5

10 L2 norm of error

Operator Split Solution

Example 1.1. We illustrate the instability of operator splitting applied to the Brusselator problem in Fig. 1.1. We apply a standard first order splitting scheme to a

1 0.5 0

10 10 10

0

-1

-2

-3

-4

pe

slo

-0.5 0

0.2

0.4 0.6 0.8 Spatial Location

1

10

t = 6.4 t = 16 t = 32 t = 64 t = 80

1

-5 -3

10

-2

10

-1

10

0

10

1

10

Time Step Size

Fig. 1.1. The lefthand plot illustrates typical instability that can arise from multiscale operator splitting applied to Brusselator problem. Solution is shown at time 80. On the right, we show plots of the error in the L2 norm versus time step size at different times.

space discretization of the Brusselator model with 500 discrete points. (see (6.6) in Sec. 6 with α = .6, β = 2, k1 = k2 = .025) We integrate using the trapezoidal rule with time step of .2 for the diffusion and backward Euler with time step of .004 for the reaction. On the left of Fig. 1.1, we show a numerical solution that exhibits nonphysical oscillations that developed after some time. On the right, we show plots of the error versus time steps at different times. There is a critical time step above which the instability develops. Moreover, changing the space discretization does not improve the accuracy. In [28], it is demonstrated that a finer spatial discretization for a constant time step size leads to significantly more error in the long time solution. The observed instability is a direct consequence of the operator splitting, which in effect separates the stabilizing effect of the diffusion component from the destabilizing reaction component. Beginning with the classic work of Marchuk and Strang, there is a highly de-

A POSTERIORI ANALYSIS OF OPERATOR SPLITTING

3

veloped literature devoted to a priori convergence analysis of operator splitting for reaction-diffusion equations, see for example [32, 23, 5, 20, 31, 34, 8] and references therein. Almost all of the classic literature assumes that the time steps for the reaction and diffusion components are comparably sized, though Dawson and Wheeler consider multiscale operator splitting for convection-reaction-diffusion problems in [5]. There are also studies of stability properties of various operator splitting procedures, see for example [21]. All of the existing convergence analysis is carried out under assumptions that, if strictly enforced, would prevent instabilities arising from the operator splitting discretization. However, enforcing such requirements presents serious problems on both theoretical and practical grounds. This provides the motivation to seek an a posteriori error analysis that devises a computational basis for detecting and correcting instabilities and inaccuracies arising from operator splitting as a particular computation proceeds. To simplify the presentation, we consider a simple example of operator splitting applied to a model system of ordinary differential equations of the form: find y ∈ Rl , l ≥ 1 such that ( y˙ = Ay(t) + F (y(t)), 0 < t ≤ T, (1.1) y(0) = y0 , where A is an l × l constant matrix representing a “diffusion component” and F (y) = (F1 (y), F2 (y), · · · , Fl (y))> is a vector of nonlinear functions representing a “reaction component”. Example 1.2. Such a system arises after semi-discretization of a reaction-diffusion governed by the initial boundary value problem,  ∂u   x ∈ Ω, 0 < t,  ∂t = ² ∆u + f (u), (1.2) suitable boundary conditions, x ∈ ∂Ω, 0 < t,    u(·, 0) = u0 (·) where Ω ⊂ Rd is a spatial domain. For example, using a continuous, piecewise linear finite element method for a Dirichlet problem with Ne elements, we obtain the initial value problem (1.1). The unknown y consists of the vector of nodal values of the finite element solution with dimension Ne − 1 and A is a symmetric negative definite sparse matrix. Our results extend in a straightforward way to different operator splitting schemes and differential equations with a nonlinear diffusion component, e.g., y˙ = FD (y) + FR (y), albeit at the cost of complicated notation. Formally, our approach also extends directly to reaction-diffusion problems (1.2) where the splitting is carried out at the continuous differential equation level followed by the discretization independently for each component. A technical issue is dealing with the finite dimensional representation of the reaction component, which is an ordinary differential equation in a Banach space when the reaction is decoupled from the diffusion. Our main goal is to derive a computable a posteriori expression that accurately estimates the error in a specified quantity of interest computed from a multiscale operator splitting approximate solution of (1.1). The a posteriori analysis is based on

4

D. ESTEP, V. GINTING, D. ROPP, J.N. SHADID, AND S. TAVENER

duality, adjoint operators, and variational analysis [15, 13, 10, 11, 16]. To deal with operator splitting, we distinguish the effects of operator splitting from the effects of numerical discretization of the components. The latter can be treated with the standard a posteriori analysis. Estimating the error arising from operator splitting requires a new approach. A main technical issue is the definition of a suitable adjoint problem. The standard approach in error analysis is to form an implicit equation for the error, linearize this problem around an average of the true and approximate solutions using the integral mean value theorem, and then employ the unique adjoint to this linearized error equation [16]. This approach generally works for implicit discretizations, for which the differential operator and its discretization are close in some sense. This approach generally fails or requires significant modification for other discretizations, e.g., explicit, multistep, and implicit/explicit schemes. Indeed, the solution operator for many classes of numerical schemes are associated with different adjoint operators than the solution operator for the original problem. This is particularly true for operator decomposition methods, which are generally associated with radically different adjoint operators than the forward problem (see below and [2, 18]). This is one important reason, perhaps unrecognized, that previous a posteriori analyses of evolution problems have focused on fully implicit time integration methods. Because an adjoint problem carries the global stability information about the quantity of interest computed from the solution, accounting for the differences between adjoint problems associated with the original problem and a numerical discretization are critical for obtaining accurate error estimates. In the estimate derived below, this difference takes the form of “residuals” between certain adjoint operators associated with the fully coupled problem and an analytic operator split version. A practical difficulty with such a result is that solving the adjoint for the fully coupled problem poses the same multiphysics challenges as solving the original forward problem. We therefore develop a new hybrid a priori - a posteriori estimate that combines a computable leading order expression obtained using a posteriori arguments with a provably higher order bound obtained using a priori convergence result. The rest of the paper is organized as follows. In Section 2, we formulate the analytic operator splitting procedure for (1.1) and conduct a preliminary investigation of the instability in operator splitting using an illuminating “blow up” problem. We present a multiscale operator splitting Galerkin finite element method for (1.1) in Section 3. We begin the analysis by presenting the results of an a priori convergence analysis of the finite element method in Section 4. The main result on the a posteriori analysis of the operator splitting finite element method is presented in Section 5 followed by several numerical examples in Section 6. In Section 7, we give the details of proof of the a priori result. Finally, in Section 8, we present a conclusion. 2. Analytic operator splitting. In this section, we define an analytic operator splitting version of (1.1). We first discretize [0, T ] into 0 = t0 < t1 < t2 < · · · < tN = N T with diffusion time steps {∆tn }n=1 , ∆tn = tn − tn−1 , and ∆t = max1≤n≤N (∆tn ). For the first discretization, we introduce a theoretical discretization in which each component is solved exactly. We define a piecewise continuous approximate solution y˜(t) =

t − tn−1 tn − t y˜n−1 + y˜n , ∆tn ∆tn

tn−1 ≤ t ≤ tn ,

with nodal values y˜n obtained from the following procedure.

(2.1)

5

A POSTERIORI ANALYSIS OF OPERATOR SPLITTING

Algorithm 2.1 [Analytic Operator Splitting] • Set y˜0 = y0 . • For n = 1 to N – Compute y r (t− n ) satisfying the reaction component ( y˙r = f (y r (t)), tn−1 < t ≤ tn , y r (t+ ˜n−1 n−1 ) = y – Compute y d (t− n ) satisfying the diffusion component ( y˙d = Ay d (t), tn−1 < t ≤ tn , r − y d (t+ n−1 ) = y (tn )

(2.2)

(2.3)

– Set y˜n = y d (t− n) end for Example 2.2. To provide motivation, we consider a problem in which the reaction component exhibits finite time blow up when undamped by the diffusion component. The problem is ( y˙ + λy = y 2 , t > 0, (2.4) y(0) = y0 ∈ R, which has exact solution y(t) =

λy0 , y0 − (y0 − λ) eλt

(2.5)

when λ 6= 0. The exact solution exists for all time and tends to zero as t → ∞ when λ > y0 . On the other hand, there is finite time blow up if λ < y0 . Applying the analytic operator splitting (2.2), (2.3) to (2.4), the solutions of the two components and the true solution are, y r (t) =

d− yn−1 d− 1 − yn−1 (t − tn−1 )

,

y d (t) = e−λ(t−tn−1 ) ynr− ,

y˜n =

e−λ∆tn y˜n−1 , 1 − ∆tn y˜n−1

when the reaction component is defined. We see that splitting off the smoothing effect provided by instantaneous interaction with the diffusion component means that the reaction component can blow up in finite time. The different behavior of the reaction and diffusion components introduces a difference in scales. Consider λ = 2 with initial condition y0 = 1, so y(t) = 2/(1+e2t ). This solution has an asymptotic decay rate proportional to e−2t , i.e., an asymptotic decay time scale of 1/2. This same asymptotic behavior is observed for any y0 < 2. The solution to the diffusion equation, y d (t) = e−2(t−tn−1 ) ynr− , has a fixed decay rate proportional to e−2(t−tn−1 ) , i.e., a fixed decay time scale τD = 1/2. However, d− d− the solution to the reaction equation, y r (t) = yn−1 /(1 − yn−1 (t − tn−1 )), becomes d− d− unbounded at t = tn−1 + 1/yn−1 for all yn−1 > 0, and has a doubling time equal d− d− d− to tn−1 + 1/(2yn−1 ), suggesting a reaction time scale τR = 1/yn−1 . If yn−1 τD . If yn−1 ≈ λ, τR ≈ τD , but the reaction time scale is not constant and

6

D. ESTEP, V. GINTING, D. ROPP, J.N. SHADID, AND S. TAVENER

decreases rapidly as the reaction progresses independently undamped by diffusion. Operator splitting has created two uncoupled problems with very different stability properties and with very different time scales. Furthermore, the time scales are not clearly separated nor do they remain in constant proportion. The time scale τD is related to the time scale of the coupled process in the limit of small solutions, but the time scale τR is completely unrelated. On the left side of Fig. 2.1, we plot the true solution and the nodal values of the approximation y˜ for N = 50 diffusion steps and M = 1 reaction step per diffusion step. The approximation is reasonably accurate.

Solution

Solution

0.8 0.6 0.4

2

2

1.5

1.5 Solution

1

1

0.5

0.5

0.2 0

0 0

0.5

1 time

1.5

2

1

0 0

0.5

1 time

1.5

2

0

0.5

1 time

1.5

2

Fig. 2.1. Plots of the solution y˜ and the true solution. Left: N = 50, M = 1. Middle: N = 10, M = 5. Right: N = 5, M = 10. The nodal values of y˜ are denoted by the larger points while the smaller points denote node values of the reaction components on the mesh used for the reaction.

Next, we increase the diffusion step by choosing N = 10 and, in order to maintain the same resolution, we correspondingly increase to M = 5. In Fig. 2.1, we plot y˜, the reaction component y r , and the true solution. The node values of y˜ are relatively close to those of y. The subsequent nodal values of the reaction component solution y r inside each step move away from the true solution. This large departure is somewhat counteracted by application of the diffusion operator. The multiscale reaction components exhibit significant growth inside each diffusion step, which severely affects accuracy. If we increase the diffusion step by taking N = 5 and maintain resolution in the reaction component by taking M = 10, the approximation becomes even less accurate. In Fig. 2.1, we plot y˜, the reaction component y r , and the true solution. If we increase the diffusion step further, then the reaction component actually blows up inside a diffusion step. 3. A multiscale operator splitting finite element method. We first discretize [0, T ] into 0 = t0 < t1 < t2 < · · · < tN = T with diffusion time steps N {∆tn }n=1 , ∆tn = tn − tn−1 , and ∆t = max1≤n≤N (∆tn ). For each diffusion step, we choose a (small) time step ∆sn = ∆tn /Mn with ∆s = max1≤n≤N (∆sn ), and the nodes tn−1 = s0,n < s1,n < · · · < sMn ,n = tn (see Fig. 3.1). We associate the time intervals In = [tn−1 , tn ] and Im,n = [sm−1,n , sm,n ] with these discretizations. We observe that an operator splitting discretization is actually a consistent discretization of the analytic operator split problem, that is a consistent discretization of (2.2) followed by a consistent discretization of (2.3). Since the a posteriori analysis is based on the framework of variational formulation and adjoint problems, we formulate the discretization as a discontinuous Galerkin (dG) finite element method in time ([6, 7, 13, 17, 16]).

7

A POSTERIORI ANALYSIS OF OPERATOR SPLITTING t0 Diffusion Integration:

t1

1

∆t

s ∆s1

Reaction Integration: s

0,1

...

s

0,2

M1,1

∆t2

... ∆s2

s

t2

∆t

M2,2

∆s3

s

...

0,3

t3

3

s s

0,4

t4

4

∆t

...

s

5

∆t

t5

M4,4

M3,3

Fig. 3.1. Discretization of time used for multiscale operator splitting

The variational formulation of (1.1) reads, for n = 1, . . . , N , find y d ∈ C 1 (In ) such that Z Z  (y˙d , v) dt = (Ay d , v) dt for all v ∈ C 1 (In ), (3.1) I In  dn + r − y (tn−1 ) = y (tn ), where y r ∈ C 1 (In ) satisfies Z Z  (y˙r , w) dt = (F (y r ), w) dt for all w ∈ C 1 (In ), In In  r + y (tn−1 ) = y˜n−1 ,

(3.2)

−

and we set y˜n = y d n . The finite element approximate solutions are sought in a piecewise polynomial spaces, n o V (qd ) = U : U |In ∈ P (qd ) (In ), 1 ≤ n ≤ N , n o V (qr ) (In ) = U : U |Im,n ∈ P (qr ) (Im,n ), 1 ≤ m ≤ Mn , for n = 1, · · · , N , and In = [tn−1 , tn ] and Im,n = [sm−1,n , sm,n ]. P (qd ) (In ) denotes the space of polynomials in Rl of degree qd on In . A similar definition holds for P (qr ) (Im,n ). We let Un+,− denote the left- and right-hand limits of U at tn and [U ]n = Un+ − Un− the jump value of U at tn . Let Y˜ (t) be the piecewise continuous finite element approximation of the operator splitting with tn − t ˜ t − tn−1 ˜ Y˜ (t) = Yn−1 + Yn , ∆tn ∆tn

tn−1 ≤ t ≤ tn .

The nodal values Y˜n are obtained from the following procedure: Algorithm 3.1 [Multiscale Operator Splitting Finite Element Method] • Set Y˜0 = y0 . • For n = 1 to N ˜ – Set Y r − 0,n = Yn−1 .

8

D. ESTEP, V. GINTING, D. ROPP, J.N. SHADID, AND S. TAVENER

– For m = 1 to Mn compute Y r |Im,n ∈ P (qr ) (Im,n ) satisfying Z

³

Y˙ r , W

´

Im,n

¡ ¢ + dt + [Y r ]m−1,n , Wm−1

Z (F (Y r ), W ) dt for all W ∈ P (qr ) (Im,n ). (3.3)

= Im,n

end for − – Set Y d n−1 = Y r − Mn ,n . – Compute Y d |In ∈ P (qd ) (In ) satisfying Z

¡ ¢ + (Y˙d , V ) dt + [Y d ]n−1 , Vn−1 In Z = (AY d , V ) dt for all V ∈ P (qd ) (In ). (3.4) In

− – Set Y˜n = Y d n . end for

In general, the dG method using polynomials of degree q converges with order up to q + 1 at all points t while its nodal values from the left converge with order up to 2q + 1 under certain conditions, see [11, 13]. Example 3.2. For qr = qd = 0, Z F (Y r ) dt → F (Y r − m−1,n ) ∆sn , Im,n

Z

−

In

AY d dt → AY d n ∆tn .

The dG approximations for the two components are r− r− Y r− m,n = Y m−1,n + F (Y m−1,n ) ∆sn ,

−

−

−

Y d n = Y d n−1 + AY d n ∆tn .

The former is equivalent to the forward Euler scheme, while the latter to the backward Euler scheme. By employing quadrature formulas to evaluate the integrals in (3.3) and (3.4), many popular finite difference schemes can be described using this variational framework, see [6, 7, 17, 16]. Using quadrature requires a straightforward extension of the a posteriori analysis presented below. The results also extend easily to higher order dG methods as well as the continuous Galerkin (cG) method [15, 16]. Example 3.3. We can use the dG qd = 0 for the diffusion component and the cG qr = 1 implemented with the trapezoidal rule to obtain the reaction component approximation Y r m,n = Y r m−1,n +

1 ∆sn (F (Y r m−1,n ) + F (Y r m,n )) , 2

which is equivalent to the Crank-Nicholson scheme. 4. An a priori convergence analysis. In this section, we carry out an a priori convergence analysis for the multiscale operator splitting dG finite element method.

A POSTERIORI ANALYSIS OF OPERATOR SPLITTING

9

We require this convergence result for the hybrid a posteriori–a priori error analysis. The convergence analysis uses familiar tools from standard results for difference schemes. However, we carry out the analysis in an unusual way by using the analytic operator split problem (2.2) and (2.3) directly. The goal is to estimate the global error (y − Y˜ ) = (y − y˜) + (˜ y − Y˜ ), which we decompose as the sum of the error due to the analytical operator splitting and the error due to the numerical discretization of the components of the problem. We give estimates of these two errors in Theorem 4.2 and Theorem 4.4 respectively and combine these in Theorem 4.5. Pl For y ∈ Rl , we define the Euclidean norm as |y|2 = i=1 |yi |2 . 4.1. A priori error analysis of analytic operator splitting. With yn = y(tn ), we write the Taylor expansion of y(t) solving (1.1) around tn−1 as, y(t) = yn−1 + (t − tn−1 ) (Ayn−1 + F (yn−1 )) 1 + (t − tn−1 )2 (A + F 0 (yn−1 ))(Ayn−1 + F (yn−1 )) + O(∆t3n ), (4.1) 2 where F 0 (y) denotes the Jacobian of F (y). Similarly, r+ r+ yr − n = y n−1 + ∆tn F (y n−1 ) +

1 2 0 r+ 3 ∆t F (y n−1 ) F (y r + n−1 ) + O(∆tn ), 2 n

and − yd n

µ ¶ 1 2 2 + 3 = I + ∆tn A + ∆tn A + O(∆tn ) y d n−1 . 2

+

−

r+ Since y d n−1 = y r − ˜n−1 , and y˜n = y d n , we have n , y n−1 = y

y˜n = y˜n−1 + ∆tn (A˜ yn−1 + F (˜ yn−1 )) ¢ 1 2 ¡ 2 yn−1 ) + F 0 (˜ yn−1 ) F (˜ yn−1 ) + O(∆t3n ). (4.2) + ∆tn A y˜n−1 + 2AF (˜ 2 Furthermore, substitution of (4.2) in (2.1) yields a representation of y˜(t) with tn−1 ≤ t ≤ tn : y˜(t) = y˜n−1 + (t − tn−1 ) (A˜ yn−1 + F (˜ yn−1 )) ¡ ¢ 1 + ∆tn (t − tn−1 ) A2 y˜n−1 + F 0 (˜ yn−1 ) F (˜ yn−1 ) + 2AF (˜ yn−1 ) + O(∆t3n ). (4.3) 2 Lemma 4.1. In the case of y˜n−1 = yn−1 , Z 2 |yn − y˜n | ≈ O(∆tn ) and |y(t) − y˜(t)| dt ≈ O(∆t3n ). In

Proof. With y˜n−1 = yn−1 , (4.2) yields yn − y˜n =

¢ 1 2¡ 0 ∆tn F (yn−1 )Ayn−1 − AF (yn−1 )) + O(∆t3n , 2

from which we get the first part of the lemma. Moreover, subtraction of (4.3) from

10

D. ESTEP, V. GINTING, D. ROPP, J.N. SHADID, AND S. TAVENER

(4.1) with y˜n−1 = yn−1 yields 1 (t − tn−1 )2 (A + F 0 (yn−1 ))(Ayn−1 + F (yn−1 )) 2 ¡ ¢ 1 − ∆tn (t − tn−1 ) A2 yn−1 + F 0 (yn−1 ) F (yn−1 ) + 2AF (yn−1 ) + O(∆t3n ). 2

y(t) − y˜(t) =

and the second estimate in the lemma results from integration over In . Theorem 4.2 (Global analytic operator splitting error). Let y be the solution of (1.1) and y˜ be the solution of the analytical operator splitting (2.2) and (2.3). Then |yN − y˜N | ≈ O(∆t). Remark 4.1. If we use higher order splitting, analogous analysis yields a higher order accuracy. Proof. Subtraction of (4.3) from (4.1) gives yn − y˜n = (yn−1 − y˜n−1 ) + ∆tn (F (yn−1 ) − F (˜ yn−1 ) + A(yn−1 − y˜n−1 )) + O(∆t2n ), ≤ C |yn−1 − y˜n−1 | + O(∆t2n ), where C is a generic constant independent of ∆t. Summing, |yN − y˜N | ≤ C

N X

O(∆t2n ) ≤ C ∆t.

n=1

Example 4.3. We consider the blow up example in Section 2. Assuming a uniform diffusion time step, induction yields y˜N =

y0 e−λtN (1 − e−λ∆t ) . (1 − e−λ∆t ) − y0 ∆t (1 − e−λtN )

Using (2.5), we obtain an estimate of the analytic operator splitting error ¡ ¢¡ ¢ 1 − e−2tN 2∆t − (1 − e−2∆t ) |yN − y˜N | = ≤ C∆t. (1 + e2tN ) ((1 − e−2∆t ) − (1 − e−2tN ) ∆t) Table 4.1 shows the error versus the time step. Table 4.1 Analytic operator splitting errors at T=2.0 for the blow up example in Section 2

∆t 0.050000 0.025000 0.012500 0.006250 0.003125

Error 0.00185 0.00089 0.00044 0.00021 0.00010

order 1.102 1.048 1.023 1.011 1.005

4.2. Analysis of the multiscale operator splitting dG finite element method. We now turn to the analysis of the numerical solution of the multiscale operator splitting problem (2.3) and (2.2). The unusual feature of this problem is the

11

A POSTERIORI ANALYSIS OF OPERATOR SPLITTING

numerical solution of an alternating sequence of diffusion and reaction components. The analysis is carried out using the analog of the standard local error analysis for a finite difference scheme. For each component solve on each interval, we decompose the error as a sum of the error in the initial condition inherited from the previous component solve and the error of the numerical solution of the component assuming exact initial conditions on the current interval. We describe the main results below and give the detailed proof in Section 7. Theorem 4.4 (Numerical discretization error). Let y˜ be the solution of analytic operator splitting and Y˜ be the multiscale operator splitting dG finite element solution, which approximates y˜. Assume that there exists a positive constant L such that |F (u)− F (v)| ≤ L |u − v| for all u, v. Then for qd = 0, 1 and qr = 0, 1, |˜ yN − Y˜N | ≈ O(∆tqd +1 ) + O(∆sqr +1 ),

(4.4)

Proof. We set e = y˜ − Y˜ , and write − − − − en = y˜n − Y˜n = y d n − Y d n = (y d n − Xn− ) + (Xn− − Y d n ) = θn− + ζn− ,

(4.5)

where X and θ are defined in Lemma 7.5. Subtracting (3.4) from (7.6) yields Z Z + + + ˙ V ) dt − (ζ, (Aζ, V ) dt + (ζn−1 , Vn−1 ) = (er − n , Vn−1 ), In

In

r− r− (qd ) for all V ∈ P (qd ) (In ), where er − . Arguing as in the n = y n − Y n . Now, ζ ∈ V proof of Lemma 7.1, we find

|ζn− | ≤ exp(8CA ∆tn ) |er − n |, where CA = |A|. By Lemma 7.2, this inequality together with (4.5) and Lemma 7.5 give the recursive relation, ¡ ¢ |en | ≤ C ∆tqnd +2 + exp(τn )∆sqnr +1 ∆tn + exp(2τn ) |en−1 | , where τn = C∆tn . From this recursive inequality, we have |eN | ≤ C ∆tqd +1

N X n=1

exp(2(n − 1)τn )∆tn + C ∆sqr +1

N X

exp((2n − 1)τn ) ∆tn .

n=1

yielding the desired result. Combining Theorems 4.2 and 4.4 we have the following global a priori error bound. Theorem 4.5 (Global a priori error bound). Let Y˜ be multiscale operator splitting dG finite element solution. Assume that there exists a positive constant L such that |F (u) − F (v)| ≤ L |u − v| for all u, v. Then for qd = 0, 1 and qr = 0, 1, there exists constants C1 , C2 , C3 such that, |yN − Y˜N | ≤ C1 ∆t + C2 ∆tqd +1 + C3 ∆sqr +1 .

12

D. ESTEP, V. GINTING, D. ROPP, J.N. SHADID, AND S. TAVENER

5. An a posteriori–a priori analysis. We now turn to the hybrid a posteriori– a priori analysis of the multiscale operator splitting dG finite element method. As mentioned in Sec. 1, the analysis differs from the standard a posteriori analysis for an evolution problem in several ways. Mainly, we have to account for the fact that different adjoint operators are associated with the fully coupled problem (1.1) and the analytic operator split version (2.2) - (2.3). The difference between these adjoints provides the means to estimate the effects of operator splitting on the stability properties of a quantity of interest computed from a solution.

5.1. Defining adjoint problems. To define adjoint operators for the forward problems, we seek analogs of the classic representation formula involving the Green’s function of a linear elliptic problem. In order to obtain unique adjoint problems from this condition, we use linearization. We assume that y = 0 is a steady state solution of (1.1), which can be achieved by assuming that Homogeneity Assumption:

F (0) = 0,

and we linearize in a region around 0. In terms of applications to reaction-diffusion problems, there are mathematical reasons for making the homogeneity assumption [30] and it is satisfied in a great many cases. We can modify the analysis to allow for linearization around any other constant steady-state solution and more generally around any known function of time. The selection of a suitable “point” for linearization requires some insight into the behavior of the solution being approximated. We illustrate below in the bistable example presented in Sec. 6. On time interval (tn−1 , tn ), we consider the linearized problem ( y˙ = A y(t) + F 0 (y) y(t), tn−1 < t ≤ tn , y(tn−1 ) = yn−1 , where F 0 (y) is a matrix whose ij-th entry is expressed as Z F 0 (y)ij

1

= 0

∂Fi (sy) ds. ∂yj

We note that F 0 (y)y = F (y) because F (0) = 0. The generalized Green’s function ϕ satisfies the adjoint problem ( > −ϕ˙ = A> ϕ(t) + F 0 (y) ϕ(t), tn > t ≥ tn−1 , (5.1) ϕ(tn ) = ψn , >

where ψn determines the quantity of interest (y(tn ), ψn ), and A> and F 0 (y) the transpose of A and F 0 (y), respectively. RT 0

denote

Remark 5.1. The analysis extends to cover quantities of interest of the form (y, ψ) dt in a straightforward way.

We choose ψn = ϕ(t+ n ), which couples the local adjoint problems (5.1) to form a global adjoint problem. This definition yields a simple representation of the solution value over one time step. Taking the Euclidean inner product of (5.1) with y and

13

A POSTERIORI ANALYSIS OF OPERATOR SPLITTING

integrating over In gives Z ³ ´ > 0= −ϕ˙ − A> ϕ(t) − F 0 (y) ϕ(t), y dt In Z Z d =− (ϕ, y) dt + (y˙ − Ay − F 0 (y)y, ϕ)dt. In dt In

(5.2)

The second term on the right hand side vanishes by (1.1), and thus integrating the first term yields the (yn , ψn ) = (yn−1 , ϕn−1 ).

(5.3)

Using the analogous approach for each component of the analytic operator splitting (2.2) and (2.3), we define the two adjoint problems for n = 1, · · · , N , ( > −ϕ˙r = F 0 (y r ) ϕr (t), tn > t ≥ tn−1 , (5.4) r ϕr (t− n ) = ψn , ( −ϕ˙d = A> ϕd (t), d ϕd (t− n ) = ψn .

tn > t ≥ tn−1 ,

(5.5)

Note that (5.4) is linearized around y r . We obtain the following representations for the component solutions r r+ r+ (y r − yn−1 , ϕr + n , ψn ) = (y n−1 , ϕ n−1 ) = (˜ n−1 ),

−

+

+

(5.6)

+

d (y d n , ψnd ) = (y d n−1 , ϕd n−1 ) = (y r − n , ϕ n−1 ).

(5.7)

5.2. A local representation of the analytic operator splitting error. As with the a priori analysis, we begin with the decomposition y − Y˜ = (y − y˜) + (˜ y − Y˜ )

(5.8) −

In order to construct a representation of the nodal value yn − y˜n , where y˜n = y d n , it is + natural to set ψnr = ϕd n−1 . Thus, taken together, the adjoint components (5.4)-(5.5) are analogous to an analytic operator split version of the adjoint problem (5.1), except that the linearization is taken around y r instead of y. Note that the adjoint component problems (5.4)-(5.5) are solved “in the reverse order”, which mirrors the fact that the adjoint of a composition of linear transformations is equal to the composition of the adjoints of the transformations computed in the reverse order. Now we combine (5.6) and (5.7) to get −

(˜ yn , ψnd ) = (y d n , ψnd ) = (˜ yn−1 , ϕr + n−1 ).

(5.9)

To obtain a representation of the local error of operator splitting, we assume that the adjoint problems (5.1), and (5.5) and (5.4) have the same data ψnd = ψn at time tn .

14

D. ESTEP, V. GINTING, D. ROPP, J.N. SHADID, AND S. TAVENER

Using (5.3) and (5.9) we obtain the following error representation (yn − y˜n , ψn ) = (yn−1 , ϕn−1 ) − (˜ yn−1 , ϕr + n−1 ).

(5.10)

At this point, it is convenient to change notation and introduce the solution operators of the various adjoint problems. Let Φn (y) be the solution operator for the full adjoint problem (5.1) so Φn (y)ψn = ϕn−1 . For a problem of the form (1.1), µ ¶ Z > > 0 Φn (y) = exp ∆tn A + F (y) dt . In

Similarly, we set Φrn (y r ) and Φdn to be the solution operators for (5.4) and (5.5), respectively, so r+ Φrn (y r )ϕr − n = ϕ n−1

and

+

Φdn ψnd = ϕd n−1 .

For (1.1), µZ Φrn (y r ) = exp

>

F 0 (y r ) dt In

¶

¡ ¢ and Φdn = exp ∆tn A> .

r r d Under the assumptions above, ϕr + n−1 = Φn (y )Φn ψn , and (5.10) reads

(yn − y˜n , ψn ) = (yn−1 , Φn (y)ψn ) − (˜ yn−1 , Φrn (y r )Φdn ψn ). To conduct a local analysis, we assume that the forward problems have the correct data y˜n−1 = yn−1 . Thus, Theorem 5.1 (Local analytical splitting error). The splitting error over a single diffusion time step is represented by (yn − y˜n , ψn ) = (˜ yn−1 , ∆Φn ψn ),

(5.11)

with ∆Φn = Φn (y) − Φrn (y r )Φdn . We note that the first term on the right hand side of (5.11) includes the effects associated with the linearization of the nonlinear reaction function around y and y r , respectively. The goal now is to use (5.11) to derive a computable estimate. 5.3. An error representation for the multiscale operator splitting dG finite element solution. We now derive error estimate for the numerical discretization of the analytic operator splitting. Let ϑd define the adjoint solution associated with (3.4) (diffusion component) satisfying ( −ϑ˙d = A> ϑd (t), tn > t ≥ tn−1 , ϑd (t− n ) = ψn . Furthermore, let ϑr define the adjoint solution associated with (3.3) (reaction component) satisfying ( −ϑ˙r = (Fˆ 0 (y r , Y r ))> ϑr (t), sm,n > t ≥ sm−1,n , r ϑr (sm,n ) = ψm,n ,

15

A POSTERIORI ANALYSIS OF OPERATOR SPLITTING +

for m = Mn , · · · , 1, with ψ r Mn ,n = ϑd n−1 and ψ r m,n = ϑrm,n for m < Mn . Thus ϑr is continuous across the internal reaction time nodes sm,n , m = 1, · · · , Mn − 1. Here Fˆ 0 (y r , Y r ) is the Jacobian matrix linearized along the trajectory zy r + (1 − z)Y r , with entries Z 1 ∂Fi Fˆij0 = (zy r + (1 − z)Y r ) dz. 0 ∂yj

We set ed = y d − Y d and er = y r − Y r . Obviously on time interval In , Z Z d ˙d > d (e , ϑ + A ϑ ) dt = 0 and (y˙d − Ay d , ϑd ) dt = 0. In

In

Using integration by parts Z Z − + + (ed , ϑ˙d ) dt = (ed n , ψn ) − (ed n−1 , ϑd n−1 ) − In

(e˙d , ϑd ) dt, In

we obtain − (ed n , ψn )

+ + (ed n−1 , ϑd n−1 )

=

Z −

(Y˙d − AY d , ϑd ) dt.

In +

Furthermore, taking V = Πϑd ∈ V (qd ) in (3.4) and using the fact that ed n−1 = d er − Mn ,n − [Y ]n−1 , yields − (ed n , ψn )

=

d+ (er − n , ϑ n−1 ) −

Z In

+ + (Y˙d − AY d , ϑd − Πϑd ) dt − ([Y d ]n−1 , ϑd n−1 − Πϑd n−1 ).

To get a representation for er − n , we use similar arguments to obtain à ! Z Mn X r− r− r+ r+ 0 r r r r r ˆ ˙ 0= (e m,n , ϑ m,n ) − (e m−1,n , ϑ m−1,n ) − (e − (F (y , Y ))e , ϑ ) dt Im,n

m=1

Since (Fˆ 0 (y r , Y r ))er = F (y r ) − F (Y r ), and y r satisfies (2.2), this reduces to à ! Z Mn X r− r− r+ r+ r r r ˙ 0= (e m,n , ϑ m,n ) − (e m−1,n , ϑ m−1,n ) + (Y − F (Y ), ϑ ) dt . Im,n

m=1

Taking V = Πϑr ∈ V (qr ) (In ) in (3.3), we then have 0=

Mn X m=1

Z r− (er − m,n , ϑ m,n )

−

r+ (er − m−1,n , ϑ m−1,n )

(Y˙ r − F (Y r ), ϑr − Πϑr ) dt

+ Im,n

r+ + ([Y r ]m−1,n , ϑr + m−1,n − Πϑ m−1,n ).

16

D. ESTEP, V. GINTING, D. ROPP, J.N. SHADID, AND S. TAVENER

r+ r− Now since ϑr − ˜n−1 − Y˜n−1 , we get m,n = ϑ m,n , e 0,n = y +

d (er − yn−1 − Y˜n−1 , ϑr + 0,n ) Mn ,n , ϑ n−1 ) = (˜ Ã ! Z Mn X r r r r r+ r+ r ˙ − (Y − f (Y ), ϑ − Πϑ ) dt + ([Y ]m−1,n , ϑ m−1,n − Πϑ m−1,n ) . Im,n

m=1

− Note that ed n = y˜n − Y˜n , and thus we can combine all these representations into the following result. Theorem 5.2 (Local discretization error). The discretization error over a single diffusion step is represented by

(˜ yn − Y˜n , ψn ) = (˜ yn−1 − Y˜n−1 , ϑr + 0,n ) Ã ! Z Mn X + + (Y˙ r − F (Y r ), ϑr − Πϑr ) dt + ([Y r ]m−1,n , ϑr m−1,n − Πϑr m−1,n ) − Im,n

m=1

Z

+ + (Y˙d − AY d , ϑd − Πϑd ) dt − ([Y d ]n−1 , ϑd n−1 − Πϑd n−1 ).

− In

(5.12)

5.4. A computable error representation. In light of (5.12) and (5.11), which represent the errors over a single diffusion step, and using y˜n−1 = yn−1 , we obtain the recursive relation (yn − Y˜n , ψn ) = (yn−1 − Y˜n−1 , ϑr + 0,n ) + (yn−1 , ∆Φn ψn ) ÃZ ! M n X + + r r r r r r − (Y˙ r − F (Y ), ϑ − Πϑ ) dt + ([Y ]m−1,n , ϑ m−1,n − Πϑ m−1,n ) m=1

Im,n

Z

− In

+ + (Y˙d − AY d , ϑd − Πϑd ) dt − ([Y d ]n−1 , ϑd n−1 − Πϑd n−1 ). (5.13)

Undoing (5.13) after choosing the adjoint data ψn−1 = ϑr + 0,n for n = 2, · · · , N − 1, and summing the resulting expressions up to N , yields the following result. Theorem 5.3. Given an adjoint data ψN , the operator splitting error at the final time can be represented as (yN − Y˜N , ψN ) = −

Mn µ Z N X X n=1 m=1

−

N X

³

Y˙ r − F (Y r ), ϑr − Πϑr

Im,n N µZ X n=1

(yn−1 , ∆Φn ψn )

n=1

³ In

´

¡ ¢ r+ dt + [Y r ]m−1,n , ϑr + m−1,n − Πϑ m−1,n

¶

´ ³ ´¶ d d+ d+ d d d ˙ d Y − AY , ϑ − Πϑ dt + [Y ]n−1 , ϑ n−1 − Πϑ n−1 ˜ 1 + Q2 + Q3 , (5.14) =Q

17

A POSTERIORI ANALYSIS OF OPERATOR SPLITTING

where ψn−1 = ϑr + 0,n for n = 2, · · · , N − 1. ˜ 1 represents the accumulated effects of operator splitting measured using the Q adjoint residual ∆Φn ψn weighted by the solution values yn−1 . Note that if y tends to zero as time passes, then the accumulated effects of operator splitting decrease correspondingly. On the other hand, if y is increasing, then the effects of operator splitting also increase. The quantities Q2 and Q3 are typical a posteriori error representations for the dG finite element approximate solutions of the component problems and are computable up to the linearization of the error equation (which typically does not cause difficulty [16, 9]). Note that we expect that the dG approximations to be much closer to the analytic operator split solutions, O(∆sqr +1 ) and O(∆tqd +1 ) respectively, than to the solution of the full problem. ˜ 1 contains quantities that are not computable, namely the exact Unfortunately, Q solution of the fully coupled problem (1.1) and the associated adjoint operator. Even if we have the true solution, we cannot expect to solve the adjoint of the full problem since this nominally presents the same difficulties as the original forward problem. ˜ 1 into a computable leading order expression and a remainder We now decompose Q that is provably higher order. As a first step, (yN − Y˜N , ψN ) =

N ³ N ³ ´ ´ X X ˜ Yn−1 , ∆Φn ψn + Q2 + Q3 + yn−1 − Y˜n−1 , ∆Φn ψn n=1

n=1

= Q1 + Q2 + Q3 + Q4 , (5.15) where Q2 and Q3 are the same as for (5.14), and the new expression Q1 does not depend directly on the solution values. Next, we derive a computable estimate for the factor ∆Φn ψn in Q1 and then use asymptotic analysis to prove that the remainder from Q1 after substituting the estimate as well as Q4 are higher order. Example 5.4. To motivate (5.15) as a natural representation, we consider the blow up example in Sec. 2 once again. Since there is no numerical solution involved, the expressions Q2 and Q3 drop out. The adjoint associated with the blow up problem is ϕ(t) = Φn (t) ψn ,

tn > t ≥ tn−1 ,

Φn (t) =

yn−1 − (yn−1 − λ)eλ(t−tn−1 ) . yn−1 − (yn−1 − λ)eλ∆tn

Since the diffusion component is linear, the associated adjoint is readily seen to be ϕd (t) = Φdn (t) ψnd ,

Φdn (t) = eλ(t−tn ) .

tn > t ≥ tn−1 ,

For the reaction component, the associated adjoint is ϕ(t) = Φrn (t) ψn ,

tn > t ≥ tn−1 ,

Φrn (t) =

d− 1 − yn−1 (t − tn−1 ) d− 1 − yn−1 ∆tn

,

for t ∈ In . At the final time level tN , (5.15) reads, yN − y˜N =

N X n=1

d− yn−1 ∆Φn (tn−1 ) ψn +

N X ¡

¢ d− yn−1 − yn−1 ∆Φn (tn−1 )ψn ,

n=1

where ψN = 1 and ψn = Φrn+1 (tn ) Φdn+1 (tn ) ψn+1 for 1 ≤ n < N .

18

D. ESTEP, V. GINTING, D. ROPP, J.N. SHADID, AND S. TAVENER

Table 5.1 illustrates the performance of the a posteriori error representation. We use λ = 0.9 and initial condition y0 = 1. The column with title “order” shows the computed order of convergence of Q4 to zero. Table 5.1 Operator splitting estimated errors at T=2.0 for blow up example in Section 2

∆t 1/80 1/160 1/320 1/640 1/1280

Exact -0.176863 -0.085044 -0.041723 -0.020667 -0.010285

Q1 -0.183464 -0.086602 -0.042101 -0.020760 -0.010309

Q4 0.006600 0.001557 0.000378 0.000093 0.000023

order 2.176 2.083 2.040 2.020 2.010

This example suggests that ∆Φn = Φn (y) − Φrn (y r )Φdn ≈ O(∆t2n ).

(5.16)

If this is true, then Q4 is O(∆tqd +2 ) and can be ignored in the asymptotic limit when computing an estimate. It helps to decompose (5.16) into two separate quantities, ¡ ¢ ¡ ¢ Φn (y) − Φrn (y r )Φdn = Φn (y) − Φrn (y)Φdn + (Φrn (y) − Φrn (y r )) Φdn = E1 + E2 . E1 measures the splitting error around the true solution, while E2 measures the effect of switching the linearization from around the true solution y to around y r in the nonlinear reaction. We next derive an asymptotic representation for E1 . Lemma 5.5. E1 = E 1 + O(∆t3n ), with ¡ ¢ 1 y ) − R(˜ y ) A> and R(y) = E 1 = ∆tn A> R(˜ 2 Proof. Using the expansion exp(B) = be expressed as: E1 =

P∞

1 i i=0 i! B

Z

>

F 0 (y) dt. In

for a square matrix B, E1 can

¢ 1¡ ∆tn A> R(y) − R(y) ∆tn A> + O(∆t3n ). 2

(5.17)

E1 can be approximated by the first term in (5.17) up to O(∆t3n ). For practical reasons, the dependence on y has to be avoided. A natural candidate to replace y is the operator splitting solution y˜, which is approximated to high order by the operator splitting dG finite element method. By adding and subtracting appropriate terms, we write ¡ ¢ 1 1 1 ∆tn A> R(˜ y ) − R(˜ y ) A> + ∆tn A> (R(y) − R(˜ y )) + ∆tn (R(˜ y ) − R(y)) A> . 2 2 2 (5.18) It is obvious that the first term is computable. The task is to demonstrate that the E1 ≈

19

A POSTERIORI ANALYSIS OF OPERATOR SPLITTING

remaining terms in (5.18) are higher order. Specifically, we show that ³ ¡ ¢ ´ Y˜n−1 , A> (R(y) − R(˜ y ) ψn ≈ O(∆t2n ),

(5.19)

and ³

´ Y˜n−1 , (R(˜ y ) − R(y)) A> ψn ≈ O(∆t2n ).

(5.20)

For (5.19), we may write ´ ³ ´ ³ Y˜n−1 , A> (R(y) − R(˜ y ))ψn = AY˜n−1 , (R(y) − R(˜ y ))ψn ≤ |A|l∞ |Y˜n−1 |l∞ |R(y) − R(˜ y )|l1 |ψn |l1 , where we have used the usual vector norms and the corresponding subordinate matrix norms. In particular, since R(y) involves matrix transpose, l ¯Z X ¯ ¯ |R(y) − R(˜ y )|l1 = max ¯ 1≤i≤l

(F 0

ij (y)

−

F0

In

j=1

¯ Z ¯ ¯ y )) dt¯ ≤ C ij (˜

|y − y˜| dt, In

where we have used the Lipschitz continuity of F 0 ij . By Lemma 4.1 this yields (5.19). Similar arguments hold for (5.20) and the proof is complete. Next, we seek an approximate estimate for E2 that does not depend on the unknown solution y. Lemma 5.6. Assume that O(∆t3n ),

E2 + with E 2 = Proof. We write

∂Fi ∂yj

(Φrn (˜ y)

−

∂ 2 Fi ∂yj ∂yk are Φrn (y r )) Φdn .

and

Lipschitz continuous. Then, E2 =

E2 = (Φrn (˜ y ) − Φrn (y r )) Φdn + (Φrn (y) − Φrn (˜ y )) Φdn . We want to show that (Φrn (y) − Φrn (˜ y )) Φdn ≈ O(∆t3n ). This requires estimating the difference between the solution operators of the reaction component adjoints corresponding to linearization around y and y˜. The solution operator Φrn (y) is associated with the problem ( > −ω˙ = F 0 (y) ω(t), tn > t ≥ tn−1 , (5.21) ω(tn ) = Φdn ψn , while the solution operator Φrn (˜ y ) is associated with ( > −ω ˜˙ = F 0 (˜ y) ω ˜ (t), tn > t ≥ tn−1 , d ω ˜ (tn ) = Φn ψn .

(5.22)

We expect that if y˜ closely approximates y, then the solutions of these differential equations are close as well. In fact, Lemma 4.1 implies that |yn − y˜n | ≈ O(∆t2n ) for y˜n−1 = yn−1 .

20

D. ESTEP, V. GINTING, D. ROPP, J.N. SHADID, AND S. TAVENER

Using a Taylor expansion, we write dω(tn ) 1 2 d2 ω(tn ) + ∆tn + O(∆t3n ) dt 2 dt2 Ã !! > 0 (y ) > > 2 1 dF n 2 I + ∆tn F 0 (yn ) − ∆tn + [F 0 (yn ) ] Φdn ψn + O(∆t3n ). 2 dt

ω(tn−1 ) = ω(tn ) − ∆tn à = Similarly, à ω ˜ (tn−1 ) =

I+

∆tn F 0 (y˜n )

>

1 − ∆t2n 2

Ã

dF 0 (y˜n ) dt

!!

>

+

> [F 0 (y˜n ) ]2

Φdn ψn + O(∆t3n ).

In these two expansions, we have applied the initial condition ω(tn ) = ω ˜ (tn ) = Φdn ψn . Taking the difference between the expansions, ω(tn−1 ) − ω ˜ (tn−1 ) µ ´ > > > > 1 d ³ 0 = ∆tn (F 0 (yn ) − F 0 (˜ yn ) ) + ∆t2n F (˜ yn ) − F 0 (yn ) 2 dt ¶ > 2 > 2 1 2 0 0 + ∆tn ([F (y˜n ) ] − [F (yn ) ] ) Φdn ψn + O(∆t3n ). (5.23) 2 >

>

>

>

Now, we show that F 0 (yn ) −´F 0 (˜ yn ) is O(∆t2n ), and that ([F 0 (y˜n ) ]2 − [F 0 (yn ) ]2 ) ³ > > and d/dt F 0 (˜ yn ) − F 0 (yn ) are at least O(∆tn ). In elemental form, we may write Z F 0 (yn )ij

−

F 0 (˜ yn )ij

1

= 0

µ

∂Fi (syn ) ∂Fi (s˜ yn ) − ∂yj ∂yj

¶ ds ≤ C|yn − y˜n | ≈ O(∆t2n ),

where we have used the assumption that ∂Fi (y)/∂yj is Lipschitz continuous. Moreover, we may write µ ¶ µ ¶ > > > > [F 0 (y˜n ) ]2 − [F 0 (yn )]2 = F 0 (y˜n ) − F 0 (yn ) F 0 (y˜n ) µ ¶ > > + F 0 (yn ) F 0 (y˜n ) − F 0 (yn ) . Arguing as above shows this term is O(∆t2n ), and thus the contribution of the third term in (5.23) is O(∆t4n ). The same conclusion about the second term in (5.23) follows 2 Fi using the assumption that ∂y∂j ∂y is Lipschitz continuous. This finally leads to the k conclusion that ω(tn−1 ) − ω ˜ (tn−1 ) = O(∆t3n ). Finally, Lemmas 5.5 and 5.6 involve y˜ and y r , which we do not have in practice. We use the same kinds of arguments to replace them by the approximations Y˜ and Y r at the cost of further higher order expressions. The result is contained in Theorem 5.7. Theorem 5.7 (Computable a posteriori error estimate). A computable error

21

A POSTERIORI ANALYSIS OF OPERATOR SPLITTING

representation for the multiscale operator splitting dG finite element method is (yN − Y˜N , ψN ) = −

Mn N X X

ÃZ

n=1

!

(Y˙ r − F (Y r ), ϑr − Πϑr ) dt + ([Y

N µZ X

n=1

(Y˜n−1 , (E 1 + E 2 )ψn )

Im,n

n=1 m=1

−

N X

In

r ]m−1,n , ϑr + m−1,n

−

Πϑr + m−1,n )

¶ + + (Y˙d − AY d , ϑd − Πϑd ) dt + ([Y d ]n−1 , ϑd n−1 − Πϑd n−1 )

+ O(∆tqd +2 ) + O(∆t ∆sqr +1 ), where E1 =

³ ´ 1 ∆tn A> R(Y˜ ) − R(Y˜ ) A> , 2

³ ´ and E 2 = Φrn (Y˜ ) − Φrn (Y r ) Φdn .

We have finally obtained the hybrid a posteriori - a priori error estimate, where the leading order expressions are computable and the remainder is provably higher order. Recall that the term E 1 is a matrix that consists of the diffusion matrix A> and R the Jacobian of the nonlinear reaction, R(Y˜ ) = In F 0 (Y˜ )dt. These are computable and do not require the solution of an adjoint problem. On the other hand, the term E2 involves the computation of three adjoint solutions, namely one adjoint solution of the diffusion component, and the adjoint solutions of the reaction components linearized about Y˜ and Y r , respectively.

5.5. Remark on the definition of the adjoint problem. The homogeneity assumption F (y) = 0 implies that y = 0 is a steady state solution of (1.1). We can also define adjoint problems using other steady state solutions or a given function of time. We illustrate by assuming that c is a vector such that General Homogeneity Assumption:

F (c) = 0.

We define the variable z = y − c and set Z 1 0 F (z) = F 0 (sz + c) ds, 0

so that F 0 (z)z = F (y) − F (c) = F (y). We can define a differential equation governing z, z˙ = A z(t) + F (y) + Ac = A z(t) + F 0 (z) z(t) + Ac.

(5.24)

The associated generalized Green’s function ϕ satisfies the adjoint differential equation >

−ϕ˙ = A> ϕ(t) + F 0 (z) ϕ(t).

(5.25)

22

D. ESTEP, V. GINTING, D. ROPP, J.N. SHADID, AND S. TAVENER

Arguing as in Sec. 5.1 yields Z ¡ ¢ 0= −ϕ˙ − A> ϕ(t) − (f 0 (z))> ϕ(t), z dt In Z Z Z d 0 =− (ϕ, z) dt + (z˙ − Az − F (z)z − Ac, ϕ)dt + (Ac, ϕ)dt. In dt In In

(5.26)

The second term on the right hand side vanishes by (5.24), and thus integrating the first term yields the Z (zn , ψn ) = (zn−1 , ϕn−1 ) + (Ac, ϕ)dt. (5.27) In

Using the analogous approach for each component of the analytic operator splitting (2.2) and (2.3), we define respectively, z r = y r − c and z d = y d − c, and the associated differential equations z˙r = F (y r ) = F 0 (z r )z r ,

z˙d = Az d + Ac.

To these differential equations, we define the respective adjoint problems, >

−ϕ˙r = F 0 (z r ) ϕr (t),

−ϕ˙d = A> ϕd (t),

Again, arguing as in Sec. 5.1, we obtain the following representations for the component solutions r r+ r+ (z r − n , ψn ) = (z n−1 , ϕ n−1 ), Z ¡ ¢ − + + (z d n , ψnd ) = (z d n−1 , ϕd n−1 ) + Ac, ϕd dt. In

Next we can follow the analysis in Sec. 5.2 to get an error representation analogous to Theorem 5.1. Theorem 5.8. Let c be a vector such that F (c) = 0. Then the splitting error over a single diffusion time step can be represented by Z ¡ ¢ (yn − y˜n , ψn ) = (˜ yn−1 − c, ∆Φn ψn ) + Ac, ϕ − ϕd dt In

with ∆Φn = Φn (y) − Φrn (y r )Φdn . We note that when c = 0, the original error representation in Theorem 5.1 is recovered. Also, we showed above that the first term involving ∆Φn ψn can be cast into expressions that are computable. Similar analysis can be carried out for the last term involving ϕ − ϕd . To see this, we write ϕ − ϕd = (ϕ − ω) + (ω − ω ˜ ) + (˜ ω − ϕd ),

(5.28)

where ω and ω ˜ are the functions governed by differential equations (5.21) and (5.22). Using the expansion of ϕ and ω gives ¡ ¢ ϕ(t) − ω(t) = (tn − t) A> ψn − ∆tn A> ψn + O (tn − t)2 + ∆t2n (5.29) ¡ ¢ = (tn−1 − t) A> ψn + O (tn − t)2 + ∆t2n ,

23

A POSTERIORI ANALYSIS OF OPERATOR SPLITTING

from which we obtain Z ¡ ¢ 1 (Ac, ϕ − ω) dt = − ∆t2n Ac, A> ψn + O(∆t3n ). 2 In Using similar argument as in the proof of Lemma 5.6, we can deduce that Z (Ac, ω − ω ˜ ) dt = O(∆t3n ).

(5.30)

(5.31)

In

Finally, the function ω ˜ is actually computed in the error estimator and thus we may write Z Z ¢ ¡ ¢ ¡ ¢ ¡ 1 Ac, ϕ − ϕd dt = − ∆t2n Ac, A> ψn + Ac, ω ˜ − ϕd dt + O(∆t3n ). (5.32) 2 In In We note that in this setting, the rest of the error estimator components remains as before. 6. Numerical examples. In this section, we present several numerical examples that show the performance of the error estimates. All of these problems originate as reaction-diffusion initial boundary value problems. Unless otherwise noted, we use a continuous, piecewise linear Galerkin finite element method with Ne elements to discretize space. All forward problems are solving using the lowest order, piecewise constant dG method, which is equivalent to backward Euler scheme, while the adjoint solutions are computed using a second order, piecewise linear, continuous Galerkin method, which is equivalent to Crank Nicholson scheme. 6.1. Assessment of the operator splitting exact error. To gain some insights on the splitting procedure, we solve a nonlinear initial boundary value problem described in (1.2) with the following specification: f (x, u) =

u2 , sin(πx)

ud = 0,

u0 (x) = sin(πx).

After spatial discretization, the vector of initial conditions g has entries gj = sin(πxj ) = sin(jπ/Ne ). while using the trapezoidal rule to approximate the integrals involving the forcing terms gives Fi (y) ≈

1 yi2 1 yi2 yi2 (0 + )+ ( + 0) = . 2 sin(πxi ) 2 sin(πxi ) sin(πxi )

(6.1)

The eigenvalues of the matrix A are λi =

2k 2k − 2 cos(iπ/Ne ), h2 h

(6.2)

with the associated eigenvectors Vi whose entries are Vij = sin(ijπ/Ne ).

(6.3)

We note that g = V1 and we assume that the exact solution of (1.1) has entries, y(t) = ξ(t)V1

(6.4)

24

D. ESTEP, V. GINTING, D. ROPP, J.N. SHADID, AND S. TAVENER

0

0

N=20, M=100 N=40, M=50 N=80, M=25

-1

-1.5

-2

M=25 M=50 M=100 M=200

-0.5 Exact Error

Exact Error

-0.5

-1

-1.5

0

0.2

0.4 0.6 Component

0.8

-2

1

0

.2 0

0.4 0.6 Component

0.8

1

Fig. 6.1. Comparison of exact errors at T = 2 against spatial location: reaction time step is kept constant (left), diffusion time step is kept constant (right)

where ξ(t) is a scalar function to be determined later and V1 is the eigenvector of matrix D associated with the eigenvalue λ1 . Consequently, we have Fi (y) = Fi (ξV1 ) =

(ξV1i )2 = ξ 2 V1i , sin(πxi )

(6.5)

using fact that V1i = sin(iπ/Ne ) = sin(πxi ). Substitution of (6.5) and (6.4) into (1.1) yields the ordinary differential equation (noting that AV1 = λ1 V1 ), ˙ + λ1 ξ(t) = ξ 2 (t), ξ(t) with ξi (0) = 1. This has analytical solution ξ(t) =

λ1 . 1 − (1 − λ1 )eλ1 t

Fig. 6.1 shows exact errors of the splitting procedure for this problem. The plot on the left shows the exact error at T = 2 plotted against the spatial location while the reaction time step is kept constant. This demonstrates that the error behaves like O(∆t). The plot on the right shows the exact error at T = 2 while the diffusion time step is kept constant. As the reaction time step is decreased, the error also decreases until it reaches a condition where the dominating error comes from the diffusion time step. After this, there is no improvement in accuracy. This result confirms the a priori analysis established in the previous section that the error in the operator splitting is dominated by the component that is of O(∆t). 6.2. Performance of the a posteriori error representation. The main purpose of the next four examples is to test the accuracy of the computable parts of the hybrid a posteriori - a priori error estimate where we drop the uncomputable higher order terms. Since we do not have true solutions for most of these problems, we approximate exact errors by computing the difference between the operator splitting solution with numerical solution of the fully coupled problems computed on a very fine mesh. We refer to this approximated error as the “error”. The first three examples fall into the category of stable diffusion interacting with unstable reaction. The reaction components in the three examples present a range

25

0 0 Error and Estimate

Error, Estimate, and Contributions

A POSTERIORI ANALYSIS OF OPERATOR SPLITTING

-0.05 -0.1 Estimate O/S First Component O/S Second Component Diffusion component Reaction component

-0.15 -0.2 -0.25

0

0.2

0.4 0.6 Component

-0.2 -0.4 -0.6 -0.8

0.8

1

0

1

0.5

1.5

Time

Fig. 6.2. Blow Up results: comparison of errors against the spatial location (left), time history of errors at the midpoint location (right), the dotted line is the exact error and the (+) is the estimated error

of instability. The last example is radically different in that it represent competition between stable reaction and diffusion components. 6.2.1. A blow up problem. The first tained from  ∂2u ∂u 2   ∂t − 0.05 ∂x2 = u , u(0, t) = u(1, t) = 0,   u(x, 0) = 4x(1 − x),

example is the “blow up” problem obx ∈ (0, 1), t > 0, t > 0, x ∈ (0, 1).

As mentioned in Sec. 2, the solution of the reaction component exhibits finite time blow up when undamped by the diffusion component. This is perhaps the most extreme form of instability. For this computation, we use 20 spatial finite elements. Table 6.1 shows the effectivity index of the error estimate computed at the final time T = 1. In this computation, we keep the reaction time step constant and varied the diffusion time step and number of reaction time steps. Fig. 6.2 shows a comparison of the errors computed using ∆t = 0.05 and M = 50 reaction time steps. The left plot in Fig. 6.2 compares the errors at T = 1 plotted against the spatial location. In this plot, we plot each of the expressions contributing to the overall estimate. In the plot on the right, the time history up to T = 1.3 of the error at the mid-point location is shown. Table 6.1 Operator splitting error estimate for the blow up problem at T = 1, reaction time step = 10−3

∆t 10−1 10−2 10−3

M 100 10 1

Exact Err (%) 11.07 1.35 0.45

Effectivity 1.0286 1.0067 1.0020

6.2.2. Chemical dynamics: the Brusselator problem. The Brusselator problem is represented by a coupled set of equations first introduced by Prigogine

26

D. ESTEP, V. GINTING, D. ROPP, J.N. SHADID, AND S. TAVENER .002

.001 0 -.001

.03

.001

Error and Estimate

Species 1

Error and Estimate

Error and Estimate

.002

Species 1

0 Species 2

-.001

Species 1 .01 0 -.02

Species 2 -.002 0

0.2 0.4 0.6 0.8 Component

1

-.002 0

Species 2 0.5

1 1.5 Time

2

-.04 0

0.2 0.4 0.6 0.8 Component

1

Fig. 6.3. Brusselator results. Left: comparison of errors against the spatial location at T = 2. Middle:time history of errors at the midpoint location on [0, 2]. Right: comparison of errors against the spatial location at T = 40. The dotted line is the exact error and the (+) is the estimated error

and Lefever [26] as a model of chemical dynamics,  ∂ 2 u1 ∂u1   − k = α − (β + 1)u1 + u21 u2 , x ∈ (0, 1), t > 0,  1  2  ∂t ∂x  2  ∂u2 ∂ u2 x ∈ (0, 1), t > 0, − k2 = βu1 − u21 u2 , 2 ∂t ∂x     u1 (0, t) = u1 (1, t) = α, u2 (0, t) = u2 (1, t) = β/α, t > 0,  u (x, 0) = u (x), u (x, 0) = u (x), x ∈ (0, 1), 1 1,0 2 2,0

(6.6)

where u1 and u2 are the concentration of species 1 and 2, respectively. We use α = 2, β = 5.45, k1 = 0.008, k2 = 0.004 and initial conditions u1 (x, 0) = α + 0.1 sin(πx) and u2 (x, 0) = β/α + 0.1 sin(πx), which yields an oscillatory solution. In this case, the reaction is very mildly unstable, with at most polynomial rate accumulation of perturbations as time passes. We use a 32 node spatial finite element discretization, resulting in an ODE system with dimension 62. We note that in original form, the reaction terms do not satisfy the requirement F (0) = 0. Instead, we use a vector c as discussed in Sec. 5.5 with ci = α for i = 1, · · · , Ne − 1 and ci = β/α for i = Ne , · · · , 2Ne − 2, so that F (c) = 0. Fig. 6.3 compares the errors computed using ∆t = 0.01 and M = 10 reaction time steps to the hybrid a posteriori error estimates. We show results for [0, 2], when the solution is still in a transient stage, and at T = 40 when the solution has become periodic. All the results show that the exact and estimated errors are in remarkable agreement. 6.2.3. Chaotic dynamics: the Lorenz problem. In this example, we treat the chaotic Lorenz equations as a reaction component of a reaction-diffusion problem,  ∂ 2 u1 ∂u1   − 0.5 = 10(u2 − u1 ), x ∈ (0, 1), t > 0,  2  ∂t ∂x   2  ∂ u ∂u 2 2   − 0.5 = u1 (28 − u3 ) − u2 x ∈ (0, 1), t > 0,  2 ∂t ∂x 2 (6.7) ∂u3 ∂ u3 8  − 0.5 = u1 u2 − u3 x ∈ (0, 1), t > 0,   2  ∂t ∂x 3    u (0, t) = u (1, t) = 0, t > 0, i = 1, 2, 3, i i    ui (x, 0) = .1 sin(πx), x ∈ (0, 1), i = 1, 2, 3,

27

0.4

0.15 0.1

0.05 0

0.2 0.1

Error and Estimate

0.4 0.2 0 -0.2

Species 1

-0.4 0

0.5

1 1.5 Time

1 0.5 0 -0.5 0

2

0.1 0.05 0

0.5

1 1.5 Time

0 0.2 0.4 0.6 0.8 1 Component 1

Species 2

1.5

Species 3

0.15

0 0.2 0.4 0.6 0.8 1 Component

2

0.6 Error and Estimate

0.3

0

0 0.2 0.4 0.6 0.8 1 Component

0.2

Species 2

Error and Estimate

Species 1

Error and Estimate

0.2

Error and Estimate

Error and Estimate

A POSTERIORI ANALYSIS OF OPERATOR SPLITTING

0.5 0 -0.5 -1 -1.5 -2

2

Species 3 0

0.5

1 1.5 Time

2

Error and Estimate

Fig. 6.4. Lorenz results. Top row: comparison of errors against the spatial location. Bottom row: time history of errors at the midpoint location. The dotted line is the exact error and the (+) is the estimated error

0.2

0.25

0.1

0.125

0.2

0.1

Species 2

Species 1 0

0

.5

1

0

0

.5

Species 3

0 1

0

.5

1

Component Fig. 6.5. Accuracy of the error estimate for the Lorenz problem at time T = 16, the dotted line is the exact error and the (+) is the estimated error

The chaotic nature of the Lorenz equations implies that errors accumulate locally at an exponential rate in an average sense, placing this instability between that of the blow up and the Brusselator problems. Here, we use 20 spatial finite elements, resulting in system of ODEs with dimension 57. Fig. 6.4 shows typical error estimates for the three species. Across the top row, we compare estimates and errors at the final time T = 2, using ∆t = 5 × 10−3 and M = 10 reaction time steps. We show the time history at the midpoint across the bottom row. Again the figures show a good agreement between exact and estimated errors. We illustrate the accuracy at a longer time of T = 16 in Fig. 6.5.

28

D. ESTEP, V. GINTING, D. ROPP, J.N. SHADID, AND S. TAVENER 0 Error Estimate

Error and Estimate

0

-0.005

-0.01

“error”

“error” adjoint using 0 adjoint using 1

-0.005

-0.01

adjoint using 0 adjoint using 1 -0.015 0

0.2

0.4 0.6 Component

0.8

1

-0.015 0.2

0.25 0.3 0.35 Component

0.4

Fig. 6.6. Bistable problem error estimate at T = 41, the right plot is the zoomed view of the boxed portion on the left plot.

6.2.4. The Bistable Problem. The bistable problem,  ∂u ∂2u   − ² 2 = u − u3 , x ∈ (0, 1), t > 0,  ∂t ∂x ux (0, t) = ux (1, t) = 0, t > 0,    u(x, 0) = u0 (x), x ∈ (0, 1),

(6.8)

is a well studied example of nonlinear relaxation to equilibrium in the presence of competing stable steady states [3, 19, 1, 12, 16]. The stable steady states are u = 1 and u = −1 while u = 0 is an unstable steady state. Here, both the diffusion and the reaction are stabilizing while instability arises because of the competition between these two components. Non-equilibrium solutions are characterized by long periods of “metastability” during which the solution is nearly stationary over periods of time on the order of exp(1/²) that are punctuated by rapid transients. The spatial profile of a metastable solution consists of narrow layers between regions where the solution has values of 1 and −1. We emphasize that the interaction between reaction and diffusion is very delicate in a metastable solution, and it is by no means clear that operator splitting is a reasonable approach. For the computation we use 100 spatial finite elements. We compute the solution until final time level T = 41. Fig. 6.6 shows comparison of the error. Here, we use c = 0 and c = 1 as the linearization points. The figure shows an observable improvement when computing the error using c = 1 than c = 0. Note that the solution eventually converges to 1. 7. Details of the a priori analysis. The following section contains the sequence of lemmas required for the proof of Theorem 4.4. Lemma 7.1 determines the error in the numerical solution of the reaction component arising from errors in the initial conditions. Lemma 7.1. Let Z ∈ V (qr ) (In ) satisfy Z Z + + − + ˙ W ) dt − (Z, (F (Z), W ) dt + (Zm−1,n , Wm−1 ) = (Zm−1,n , Wm−1 ) (7.1) Im,n

Im,n

− r for all W ∈ P (qr ) (Im,n ), m = 1, 2, · · · , Mn , and Z0,n = yr + n−1 . With ξ = Z − Y , − |ξM |2 ≤ exp(16L∆tn ) |˜ yn−1 − Y˜n−1 |2 . n ,n

A POSTERIORI ANALYSIS OF OPERATOR SPLITTING

Proof. By construction, ξ ∈ V (qr ) (In ). It is obvious that for qr = 0, 1, Z 2 − ˙2 . |ξ|2 dt ≤ 2∆sn |ξm,n |2 + ∆s3n |ξ| Im,n 3 Im,n

29

(7.2)

Subtracting (3.3) from (7.1) yields Z Z + + − + ˙ (ξ, W ) dt − (F (Z) − F (Y r ), W ) dt + (ξm−1,n , Wm−1 ) = (ξm−1,n , Wm−1 ) Im,n

Im,n

(7.3) for every W ∈ P (qr ) (Im,n ). Setting W = (t − sm−1,n )ξ˙ in (7.3) and estimating yields 1 2 ˙2 ²2 2 ∆sn |ξ|Im,n ≤ L ∆sn 2 2

Z |ξ|2 dt + Im,n

1 ˙2 , ∆s2n |ξ| Im,n 2²2

for some ² > 0. Thus, Z ˙2 ∆s2n |ξ| Im,n

2

|ξ|2 dt,

≤ C² L ∆sn

(7.4)

Im,n

where C² = ²4 /(²2 − 1). Substitution of (7.4) into (7.2) gives ¶Z µ 2 − |ξ|2 dt ≤ 2∆sn |ξm,n |2 , 1 − C² ∆s2n L2 3 Im,n from which we obtain Z Im,n

− |ξ|2 dt ≤ 4∆sn |ξm,n |2 ,

(7.5)

provided 1 − 23 C² ∆s2n L2 > 1/2. Next, we choose W = ξ in (7.3) and use the Lipschitz continuity of F to get Z 1 − 2 1 + − + 2 |ξ | + |ξm−1,n | ≤ L |ξ|2 dt + |ξm−1,n | |ξm−1,n |. 2 m,n 2 Im,n This gives Z − |ξm,n |2

≤

− |ξm−1,n |2

|ξ|2 dt.

+ 2L Im,n

Using (7.5), this implies − |ξm,n |2 ≤

1 − |ξ − |2 ≤ exp(16∆sn L) |ξm−1,n |2 , 1 − 8∆sn L m−1,n

provided 1 − 8∆sn L > 1/2. Applying this recursive relation Mn times yields the desired result. Lemma 7.2 determines the error in the numerical solution of the reaction component due to discretization and the inherited initial conditions.

30

D. ESTEP, V. GINTING, D. ROPP, J.N. SHADID, AND S. TAVENER

Lemma 7.2. Let y r be the solution of analytical reaction component and Y r be the dG numerical solution approximating y r . For n = 1, · · · , N , r− qr +1 |y r − + exp(8L∆tn ) |˜ yn−1 − Y˜n−1 |. n − Y M,n | ≤ C ∆tn ∆sn

Proof. We write − − r− r− r− yr − n − Y M,n = (y n − ZM,n ) + (ZM,n − Y M,n ),

where Z ∈ V (qr ) (In ) is defined as in Lemma 7.1. Moreover, it has been established (cf. Delfour et al. [7], Johnson [22], and Estep [13]) that − qr +1 . |y r − n − ZM,n | ≤ C ∆tn ∆sn

Using Lemma 7.1, the proof is complete. We also use the following result from Delfour et al. [6]. Lemma 7.3. Let w ∈ H 1 (In ) satisfying w˙ + A> w = 0 in In . Then there exists a constant C independent of s ∈ In , such that kwkr+1,n ≤ C|w(s)|, where the norm k · kr,n is associated with the Sobolev space ( ) r ° i °2 X ° ° d w ° ° H r (In ) = w : kwk2r,n = Πw° ≤ C ∆tqn |w(s)|. ° dt ° 2 L (In ) Proof. Since w in Lemma 7.3 is homogeneous, interpolation property (cf. Strang et al [33], and Ciarlet [4]) and Lemma 7.3 imply that for tn−1 ≤ s ≤ tn , ° ° ° ° °d ° ° ° ° Πw + A> Πw° = ° d (Πw − w) + A> (Πw − w)° ° dt ° ° dt ° ≤ C kw − Πwk1,n ≤ C ∆tqn kwkq+1,n ≤ C ∆tqn |w(s)|, Lemma 7.5 bounds the error in the numerical solution to a constant coefficient differential equation. Lemma 7.5. Let X ∈ P qd (In ) satisfy Z Z + + + ˙ (X, V ) dt − (AX, V ) dt + (Xn−1 , Vn−1 ) = (y r − (7.6) n , Vn−1 ), In

In

− d − qd +2 = yr − , for for all V ∈ P qd (In ), where Xn−1 n . With θ = y − X, |θn | ≤ C ∆tn qd = 0, 1.

A POSTERIORI ANALYSIS OF OPERATOR SPLITTING

31

Proof. Obviously y d defined in (2.3) satisfies (7.6). Let U ∈ P qd (In ) satisfy Z Z + + + qd ˙ (U , V ) dt − (Ay d , V ) dt + (Un−1 , Vn−1 ) = (y r − n , Vn−1 ), for all V ∈ P (In ). In

In

(7.7)

For sufficiently small ∆tn , the following estimates holds (see [13]): |θ|In ≤ |y d − U |In + |U − X|In ≤ C∆tqnd +1 |y (qd +1) |In + C∆tnqd +3/2 |y (qd +1) |In . (7.8) Subtracting (7.6) from the variational formulation of y d gives Z Z + + ˙ V ) dt − (θ, (Aθ, V ) dt + (θn−1 , Vn−1 ) = 0. In

(7.9)

In

Integration by parts gives Z (θn− , Vn− )

(θ, V˙ + A> V ) dt.

=

(7.10)

In

For qd = 0, we choose V = Vn− = θn− in (7.9) and use (7.8) to get |θn− | ≤ C∆t2n . For qd = 1, we choose V = Πw in (7.9), with w defined as in Lemma 7.3, and w(tn ) = θn− . Then using (7.8) and Lemma 7.4 Z p − 2 ˙ + A> Πw)| dt ≤ C∆t3n ∆tn |θ− |, |θ| |Πw (7.11) |θn | ≤ n In

which gives the desired estimate. 8. Conclusion. In this paper, we derive both an a priori convergence estimate and a hybrid a posteriori - a priori estimate for the multiscale operator splitting discontinuous Galerkin finite element method for a system of ordinary differential equations. The system of equations is assumed to have a form that typically arises from space discretization of a reaction-diffusion problem. The a priori analysis uses the fact that an operator splitting approximation is a consistent approximation of an analytic operator splitting problem. The a posteriori analysis takes into account the fact that the original problem and an analytic operator split version are associated with different adjoint problems. The differences in the adjoints provides the means to accurately understand the effects of operator splitting on the stability properties of a quantity of interest computed from the solution. To produce a computable a posteriori estimate of the error, we manipulate the original a posteriori error representation to obtain a computable expression plus uncomputable quantities that are provably higher order in an asymptotic sense. We conclude with some examples that demonstrate the accuracy of the computable parts of the hybrid a posteriori - a priori estimate. Our results can be extended in a straightforward way to differential equations with a more general nonlinear form, higher order splitting schemes, different time finite element methods. The results also extend formally, and rigorously for the most part, to operator splitting finite element schemes for reaction-diffusion problems where the splitting occurs at the pde level. We will pursue the complete extension to the latter in future work. The general approach to analyzing operator splitting can be used to analyze operator decomposition in a wide variety of settings, see e.g. [2, 18].

32

D. ESTEP, V. GINTING, D. ROPP, J.N. SHADID, AND S. TAVENER REFERENCES

[1] L. Bronsard and R. Kohn, On the slowness of phase-boundary motion in one space dimension, Comm. Pure Appl. Math., 43 (1990), pp. 983–997. [2] V. Carey, D. Estep, and S. Tavener, A posteriori analysis and adaptive error control for operator decomposition methods for coupled elliptic systems I: one way coupled systems, SINUM, (2006). in revision. [3] J. Carr and R. Pego, Metastable patterns in ut = ²2 ∆u − f (u), Comm. Pure Appl. Math., 42 (1989), pp. 523–576. [4] P. G. Ciarlet, The finite element method for elliptic problems, North-Holland Publishing Co., Amsterdam, 1978. Studies in Mathematics and its Applications, Vol. 4. [5] C. N. Dawson and M. F. Wheeler, Time-splitting methods for advection-diffusion-reaction equations arising in contaminant transport, in ICIAM 91 (Washington, DC, 1991), SIAM, Philadelphia, PA, 1992, pp. 71–82. [6] M. Delfour, W. Hager, and F. Trochu, Discontinuous Galerkin methods for ordinary differential equations, Math. Comp., 36 (1981), pp. 455–473. [7] M. C. Delfour and F. Dubeau, Discontinuous polynomial approximations in the theory of one-step, hybrid and multistep methods for nonlinear ordinary differential equations, Math. Comp., 47 (1986), pp. 169–189, S1–S8. With a supplement. [8] S. Descombes, Convergence of a splitting method of high order for reaction-diffusion systems, Math. Comp., 70 (2001), pp. 1481–1501 (electronic). [9] S. Eastman, Analysis and Application of the Nonlinear Power Method, PhD thesis, Colorado State University, 2005. [10] K. Eriksson, D. Estep, P. Hansbo, and C. Johnson, Introduction to adaptive methods for differential equations, in Acta numerica, 1995, Acta Numer., Cambridge Univ. Press, Cambridge, 1995, pp. 105–158. , Computational differential equations, Cambridge University Press, Cambridge, 1996. [11] [12] D. Estep, An analysis of numerical approximations of metastable solutions of the bistable equation, Nonlinearity, 7 (1994), pp. 1445–1462. [13] D. Estep, A posteriori error bounds and global error control for approximation of ordinary differential equations, SIAM J. Numer. Anal., 32 (1995), pp. 1–48. [14] D. Estep, A short course on duality, adjoint operators, Green’s functions, and a posteriori error analysis, 2004. Sandia National Laboratories, Albuquerque, New Mexico. Notes can be downloaded from http://math.colostate.edu/~estep. [15] D. Estep and D. French, Global error control for the continuous Galerkin finite element method for ordinary differential equations, RAIRO Mod´ el. Math. Anal. Num´ er., 28 (1994), pp. 815–852. [16] D. Estep, M. G. Larson, and R. D. Williams, Estimating the error of numerical solutions of systems of reaction-diffusion equations, Mem. Amer. Math. Soc., 146 (2000), pp. viii+109. [17] D. Estep and A. M. Stuart, The dynamical behavior of the discontinuous Galerkin method and related difference schemes, Math. Comp., 71 (2002), pp. 1075–1103 (electronic). [18] D. Estep, S. Tavener, and T. Wildey, A posteriori analysis and improved accuracy for an operator decomposition solution of a conjugate heat transfer problem. submitted to SINUM, 2006. [19] G. Fusco and J. Hale, Slow-motion manifolds, dormant instability, and singular perturbations, J. Dynam. Differ. Equa., 1 (1989), pp. 75–94. [20] W. Hundsdorfer and J. Verwer, Numerical solution of time-dependent advection-diffusionreaction equations, vol. 33 of Springer Series in Computational Mathematics, SpringerVerlag, Berlin, 2003. [21] W. H. Hundsdorfer and J. G. Verwer, Stability and convergence of the Peaceman-Rachford ADI method for initial-boundary value problems, Math. Comp., 53 (1989), pp. 81–101. [22] C. Johnson, Error estimates and adaptive time step control for a class of one step methods for stiff ordinary differential equations, SIAM J. Numer. Anal., 25 (1988), pp. 908–926. [23] G.I. Marchuk, On the theory of the splitting-up method, in In: Proceedings of the Second Symposium on Numerical Solution of Partial Differential Equations, SVNSPADE, 1970, pp. 469–500. [24] G. I. Marchuk, V. I. Agoshkov, and V. P. Shutyaev, Adjoint equations and perturbation algorithms in nonlinear problems, CRC Press, Boca Raton, FL, 1996. [25] C. C. Ober and J. N. Shadid, Studies on the accuracy of time-integration methods for the radiation-diffusion equations, J. Comput. Phys., 195 (2004), pp. 743–772. [26] I. Prigogine and R. Lefever, Symmetry breaking instabilities in dissipative systems, J. Chem. Phys, 48 (1968), pp. 1695–1700.

A POSTERIORI ANALYSIS OF OPERATOR SPLITTING

33

[27] D. L. Ropp and J. N. Shadid, Stability of operator splitting methods for systems with indefinite operators: reaction-diffusion systems, J. Comput. Phys., 203 (2005), pp. 449–466. [28] D. L. Ropp, J. N. Shadid, and C. C. Ober, Studies of the accuracy of time integration methods for reaction-diffusion equations, J. Comput. Phys., 194 (2004), pp. 544–574. [29] J. Sandelin, Global Estimate and Control of Model, Numerical, and Parameter Error, PhD thesis, Department of Mathematics, Colorado State University, Fort Collins, CO 80523, 2006. [30] J. A. Smoller, Shock Waves and Reaction-Diffusion Equations, Springer-Verlag, New York, 2 ed., 1983. [31] B. Sportisse, An analysis of operator splitting techniques in the stiff case, J. Comput. Phys., 161 (2000), pp. 140–168. [32] G. Strang, On the construction and comparison of difference schemes, SIAM J. Numer. Anal., 5 (1968), pp. 506–517. [33] G. Strang and G. J. Fix, An analysis of the finite element method, Prentice-Hall Inc., Englewood Cliffs, N. J., 1973. Prentice-Hall Series in Automatic Computation. [34] J. G. Verwer and B. Sportisse, A note on operator splitting analysis in the stiff linear case, (1998).