Adaptive high-order splitting methods for systems of nonlinear evolution equations with periodic boundary conditions

We assess the applicability and efficiency of time-adaptive high-order splitting methods applied for the numerical solution of (systems of) nonlinear parabolic problems under periodic boundary conditions. We discuss in particular several applications generating intricate patterns and displaying nonsmooth solution dynamics. First we give a general error analysis for splitting methods for parabolic problems under periodic boundary conditions and derive the necessary smoothness requirements on the exact solution in particular for the Gray-Scott equation and the Van der Pol equation. Numerical examples demonstrate the convergence of the methods and serve to compare the efficiency of different time-adaptive splitting schemes and of splitting into either two or three operators, based on appropriately constructed a posteriori local error estimators.


Introduction
We are interested in computational methods for nonlinear evolution equations of the type ∂ t u(t) = Au(t) + B(u(t)), t > t 0 , (1.1) on a Banach space B, which in our examples equals L 2 on the d-dimensional torus. Here, A : D ⊆ B → B is an (unbounded) differential operator and B a generally unbounded nonlinear operator whose domain has nonempty intersection with D.
To enable an efficient numerical solution of (1.1) for large-scale applications, adaptive high-order time-discretizations are central. In some applications the promised speed-up will be critical for the feasibility of a simulation. In many realistic models, the stiffness of the operators A and B is different which suggests to use splitting methods which separately propagate the two vector fields. If A is a linear differential operator, effective schemes are known which solve the subproblem efficiently after appropriate space discretization. For the problems discussed in this paper, a Fourier pseudospectral space discretization is the most natural choice as this allows to propagate the linear part by exponentiation of a diagonal matrix.
Parabolic equations often induce high computational demand due to challenging solution dynamics, which suggests to employ adaptive time-stepping in order to accommodate for local variations in the numerical error. However, this is not the only reason for using adaptivity. Typically, the optimal step-size is not known a priori, and an adaptive procedure determines the appropriate value within a few steps, see for example Section 3.3. Moreover, adaptive time-stepping increases the reliability of a computation, see for instance [1].
At the (time-)semi-discrete level, s-stage exponential splitting methods for the integration of (1.1) use multiplicative combinations of the partial flows φ A (t, u) and φ B (t, u). For a single step (0, u 0 ) → (h, u 1 ) with time-step t = h, this reads where the coefficients a j , b j , j = 1 . . . s are determined according to the requirement that a prescribed order of consistency is obtained [2]. Compared to highly implicit methods as for instance implicit Runge-Kutta methods or their exponential counterparts (see [3]), splitting methods are easy to implement and efficient in combination with suitable spatial discretization and appropriate implementations or approximations of the subflows φ A and φ B . This is an important asset of our approach, however we will demonstrate in addition that adaptive choice of the time steps leads to a more efficient solution for problems where the variation in the solution is large. For related work on adaptivity using a pair of lower order methods we refer to [4].
A rigorous error analysis of splitting methods for Schrödinger equations has first been given for the second-order Strang splitting scheme in [5], which has later been extended to higher-order splittings in [6]. The more involved arguments for the nonlinear case have been devised in [7] for the Schrödinger-Poisson and cubic nonlinear Schrödinger equation for second order splitting; higher-order methods are analyzed in [8].
The error analysis relies on an error representation which was first proven in [8]: the local error of a splitting method of order p applied to a nonlinear evolution equation has an error expansion with leading term where C kµ are computable constants and D A , D B represent the Lie derivatives of the two vector fields, respectively. ad µ D A (D B ) denotes the µ -fold commutator. In our subsequent analysis we will make use of this error representation, where the main task will be to compute and estimate the commutators of the vector fields in an appropriate functional analytic setting in the space of periodic functions. To this end, we will resort to a Sobolev theory on the torus, which we review in detail in Appendix A, to which we refer for notations used in the subsequent error analysis.
Detailed understanding and analysis of splitting methods for parabolic problems in particular for the nonlinear case is missing to date. Partial results have been obtained by other authors; recent work for linear problems can be found in [9] and [10]. In particular, in [9], a number of higher order methods with complex coefficients are constructed. In these papers, splitting methods are analyzed in the context of semigroup theory. However, the authors do not exploit the special structure of the local error (as specified in [11] in terms of iterated commutators). Therefore the results in [10] rely on unnaturally restrictive regularity assumptions, and the same is true for the convergence results given in [9]. Section 2 introduces a number of local and global a posteriori error estimators whose performance will subsequently be assessed.
In Section 3, our theoretical framework is applied to analyze the convergence of splitting methods for the Gray-Scott equation, where the regularity requirements on the exact solution are worked out which ensure boundedness of the commutators appearing in the error expansion.
In Section 4, we investigate the Van der Pol system, which has a stiff limit cycle. Adaptive time-stepping is shown to give rise to guaranteed accuracy, and in some cases significantly reduced computation times compared to fixed time steps.
In Section 5 we demonstrate that splitting into three operators can be beneficial computationally if the structure of the vector field enables exact integration of the subproblems, by resorting to computations for the Gray-Scott equations.
The functional analytic framework for the error analysis of splitting methods applied to parabolic problems under periodic boundary conditions is briefly recapitulated in Appendix A, which states the underlying results for the space of periodic functions on the torus. Sobolev embeddings which are used in our error estimates are stated in Appendix B with a brief indication of the proofs.

A posteriori local error estimators
In this section, we briefly describe three classes of computable a posteriori local error estimators which serve as our basis for adaptive time-stepping and which have different advantages depending on the context in which they are applied. Embedded pairs of splitting formulae have been introduced in [14] and are based on reusing a number of evaluations from the basic integrator. For methods of odd order, an asymptotically correct error estimator can be computed at the same cost as for the basic method by employing the adjoint method, see [16], and finally the Milne device relies on the explicit knowledge of the leading error terms of methods of equal order. A collection of splitting coefficients covering also these three types of error estimators has been compiled at the webpage http://www.asc.tuwien.ac.at/˜winfried/splitting/ which we subsequently refer to as [17].

Embedded pairs
In [14], pairs of splitting schemes of orders p and p + 1 are specified. The idea is to select a controllerS of order p + 1 and to construct an integrator S of order p for which a maximal number of compositions coincide with those of the controller. To construct pairs offering an optimal balance between cost and accuracy, we fix a 'good' controller of order p + 1 and wish to adjoin to it a 'good' integrator of order p. Since the number of compositionss in the controller will be higher than the number of compositions s in the integrator, we can select an optimal embedded integrator S from a set of candidates obtained by flexible embedding, where the number of coinciding coefficients is not a priori fixed. The idea is expanded in detail in [16], where optimized methods are determined.

Adjoint pairs and palindromic formulae
For a scheme S of odd order p, the leading local error terms of S and its adjoint S * are identical up to the factor −1, see [2]. Therefore, the averaged additive schemē is a method of order p + 1, and provides an asymptotically correct local error estimate for S(h, u). In this case the additional effort for computing the local error estimate is identical with the effort for the integrator S but not higher as is the case for embedded pairs. This principle is limited to methods of odd order. In particular, in [16] so-called palindromic schemes were constructed which turn out to have small error constants as compared to competing schemes. Therefore, we include palindromic pairs in our investigations.

The Milne device
In the context of multi-step methods for ODEs, the so-called Milne device is a wellestablished technique for constructing pairs of schemes. In our context, one may aim for finding a pair (S,S) of schemes of equal order p such that their local errors L,L are related according to with γ = 1. Then, the additive schemē is a method of order p + 1, and provides an asymptotically correct local error estimate for S(h, u).

Step-size selection
Based on a local error estimator, the step-size is adapted such that a prescribed local error tolerance tol is expected to be satisfied in the subsequent step. If h old denotes the current step-size, the next step-size h new is predicted as (see [18,19]) where we choose α = 0.9, α min = 0.25, α max = 4.0. This simple strategy incorporates safety factors to avoid an oscillating and unstable behavior. The chosen values of α min and α max are commensurable with the recommendations in [2]. The safety factors have not proven critical in our examples, the local changes in the stepsizes are usually smaller from step to step, see for example Figure 8. Only if at the beginning of time propagation the initial stepsize is unsuitable as in Figure 6, where still no instabilities arise in the step-size control, however.

The Gray-Scott equation
As a concrete example, we first study the Gray-Scott system (see [20]) modeling a twocomponent reaction-diffusion process, This system is of the type (1.1), with unknown (u(x, y,t), v(x, y,t)), the vector of concentrations of the two chemical species involved. In many situations this model is closed naturally by periodic boundary conditions. This system is studied as a model for pattern formation with a rich dynamical behavior. For (x, y) ∈ [−4π, 4π] 2 we prescribe the initial condition A visualization of the solution component v at t = 0, 2000 and 4000 is shown in Figure 1. The problem can also naturally be stated in three spatial dimensions and solved by our methods. In Figure 2 we show the component v computed by a complex embedded 4/3 splitting pair from [14] with an underlying spatial discretization with 512 3 basis functions and a tolerance of 10 −5 . The solution is plotted at times t = 2500, t = 3000, t = 4000, and t = 5000. In the following we will only investigate the 2D case, as this does not influence the assessment of the time integrators, but reduces computation time.

Convergence analysis
For the theoretical analysis of the convergence of splitting methods, we use the error representation (1.3). Since the flow induced by the cubic nonlinearity is not unconditionally stable, we have to resort to the three-stage argument first given in [7] for the cubic Schrödinger equation, see also [8]: Along this line, we can prove the following theorem, since for the present situation of a parabolic problem under periodic boundary conditions, the same Sobolev embeddings hold as on the full space R 3 , see Appendix A, so in particular the second order differential operators and the cubic terms and their commutators admit the same bounds. Thus, the following proof strategy can be followed in the same manner, taking into account the commutator bounds given later: Theorem 1 Suppose that the Gray-Scott equation (3.1) possesses a uniquely determined sufficiently regular solution u on the time interval [0, T ]. Then, for any exponential operator splitting method (1.2) of (nonstiff) order p ≥ 2, the following error estimates are valid.
holds true with constant C depending on M 2p−2 .
Proof We work out the analysis in detail for the case p = 2, the general case is proven analogously. For the analysis, we write the Gray-Scott system in the partitioned form Stability is shown in the same manner as for the cubic Schrödinger equation [8], see the outline above. To bound the local error, we compute the commutators of the vector fields. This yields This can be estimated in Sobolev norms by resorting to the embeddings in Appendix B: For the second commutator we compute contains terms of the form uv∆ 2 u and uv∆ 2 v which do not cancel. Consequently, Inductively, the result for higher commutators appearing in estimates for higher-order splitting methods follows.

Numerical results
In this section, we will demonstrate the accuracy of several splitting schemes for the Gray-Scott equation (

Comparisons
After verifying the reliability of the investigated solution methods, we will assess the efficiency of the adaptive time integration methods by giving a comparison to the situation where the same accuracy is achieved with constant time-steps. Moreover, we will compare the efficiency of adaptive time integration based on the second order method in conjunction with the Milne device as compared to the fourth order embedded splitting pair [17,Emb 4(3) A c] and the palindromic scheme [17, PP 3/4 A c]. By construction, the latter also provides an asymptotically correct error estimator, which by its special structure is cheap to evaluate. Runtime was measured on a PC with Intel Core i7-2600 3, 4GHz Quad-Core processor with 16 GB RAM: Table 1 shows the number of steps required in the adaptive integration, the number of equidistant steps with the smallest necessary adaptive time-step, and the computing time for both scenarios. The tolerances were chosen as 10 −5 (top) and 10 −8 (bottom), respectively. We observe that indeed the adaptive methods require fewer steps, but the overall computational cost is higher due to the effort for the evaluation of the error estimator in each step. This suggests an adaptive strategy which does not estimate the error in each step, but only after a certain number of steps with a fixed time-step. This implies that an update of the time-steps every two or three steps should provide a more efficient strategy, but possibly at the cost of reduced numerical stability, since this example shows rather smooth solution dynamics. Indeed, the step-size is adjusted rapidly by exploiting the maximally permitted increase by a factor of 4 from a too small initial guess to the appropriate value, which is assumed throughout the rest of the computation, see Figure 6, which gives the quotient of two consecutive step-sizes over the integration interval. This behavior demonstrates one major advantage of adaptivity, that an unsuitable initial guess of the step-size is automatically adjusted to an optimal value.

The Van der Pol equation
The Van der Pol equation is an ordinary differential equation with limit cycle behavior. It is used as a test of time integration schemes for stiff differential equations. It shares characteristics with simple models for cardiac behavior. The Van der Pol equation is usually considered as an ordinary differential equation, but by adding diffusion terms, one can consider an extension from a set of ordinary differential equations to a pair of coupled partial differential equations with spatial dependence. It is given by It is split into and The convergence result for an order p splitting applied to this system can readily be seen to be the same as Theorem 1. However, the constants in the estimates (3.3)-(3.5) depend on the small parameter ε, (3.4), and C = C(M 2p−2 , ε 2−p ) in (3.5). We must stress that the involved estimates of the exact solution will also be negatively influenced when ε is small. The analysis of the exact solution is not a topic of the present paper, however.
For our comparisons, we solve the problem in one spatial dimension, with x ∈ [−π, π], and choose ε = 10 −3 . The evolution of the solution components with t (on the vertical axis) is illustrated in Figure 7. Results showing the effectiveness of adaptive time stepping for (4.1) are shown in Table 2. For this problem, the lower order method is more efficient. Adaptive step selection yields a speed-up by about a factor 5. Indeed, if we consider the ratio of two consecutive step-sizes, we see some variation in the region of the steep layers in Figure 8, which is obviously sufficiently large to warrant adaptive time-stepping.    The time-steps generated in the course of an adaptive procedure are given in Figure 9. The left plot shows the time-steps to satisfy a tolerance of 10 −5 for the PP 3/4 A c method, and likewise on the right for the [17, PP 5/6 A c] pair.

Splitting into three operators ('ABC-splitting')
Finally, we consider a splitting of the Gray-Scott equations (3.1) into three parts, This has the computational advantage that the flows of the operators B and C can be computed analytically when the other component is frozen. Below we verify the convergence orders for this case for the optimal palindromic splitting PP 3/4 A 3 c. Remark: A formal error analysis for ABC-splitting has not yet been given in the nonlinear case, the linear case has been treated in [15]. However, inspection of the commutators that would critically influence the error shows that a convergence result analogous to Theorem 1 will hold, since commutators of B and C vanish.

Numerical results
The numerical results below were computed by the method [17, PP 3/4 A 3 c]. This is the method of order 3 with the smallest leading error coefficients (see [16]) we could determine and offers the advantage of the cheap error estimator from Section 2.2, see Figure 10. The time-steps generated in the course of an adaptive procedure according to Section 2.4 are given in Figure 11. The plot shows the step-sizes to satisfy a tolerance of 10 −5 for the PP 3/4 A 3 c method.

Comparisons
In order to compare the efficiency of the ABC-splitting approach with the two-operator splitting discussed in Section 3, in Table 3 we give the number of steps required for tolerances 10 −5 and 10 −8 and the resulting computation times. It is observed that the ABC-splitting [17, PP 3/4 A 3 c] requires slightly fewer steps than [17, PP 3/4 A c], but the computation time is higher. The reason is that each individual step is computationally more demanding in the ABC-splitting due to the larger number of required FFT transforms associated with the larger number of compositions. These result from the fact that the number of order conditions is larger in the ABC case and therefore, more free parameters are necessary to construct high-order methods. Indeed, 1000 steps with [

Conclusions and outlook
We have investigated high-order adaptive time-splitting methods for the solution of nonlinear evolution equations of parabolic type under periodic boundary conditions. The theoretical error analysis for the Gray-Scott equations and the Van der Pol equation shows the classical convergence orders under regularity assumptions on the exact solution implied by the Sobolev inequality for functions on the torus. The theory is illustrated by numerical computations showing the established convergence orders. Moreover, adaptive time-stepping strategies have been demonstrated to improve both efficiency and reliability, where high-order methods generally yield a computational advantage for the approximation of regular solutions. Local error estimators based on embedded formulae of splitting coefficients are more efficient than estimators employing the adjoint method, but the former need to be constructed especially by a computationally demanding optimization procedure, while the latter principle can be applied invariantly for methods of odd order.
Indeed, it has been observed that for problems with rapidly varying solutions, an adaptive strategy yields an advantage as compared to uniformly using the smallest time-step required locally. Secondly, a good guess of the time step-size is not commonly available even when the solution is smooth, so adaptive adjustment saves from repeating runs until the optimal step-size is found.
Splitting into three operators promises a computational advantage for the calculation of the individual compositions, but the complexity of high-order integrators of this class implies a significant surplus of necessary compositions which negatively affects the performance.

A Periodic functions and their Fourier transforms
In the following, we recapitulate material from [21] for the convenience of the reader. Consider The space L 2 = L 2 (Q) is a Hilbert space with the inner product Fourier representation of u ∈ L 2 Let k = (k 1 , . . . , k d ) ∈ Z d , and |k| = |k 1 | + · · · + |k d |.
and the inverse transform yields the representation Parseval's identity implies an isometric correspondence Remark 1 Since the torus has finite measure, we have L q ⊆ L p for 1 ≤ p ≤ q ≤ ∞.
We introduce the following notations: H s = H s (Q), α = (α 1 , . . . , α d ) ∈ N d 0 , |α| = α 1 + · · · + α d , α! = α 1 ! · · · α d !. Weak derivatives are denoted by D α u. The norm on H s is H s is a Hilbert space with inner product Fourier representation of D α u. The weak derivative has the Fourier representation and thus as a consequence of Parseval's identity (A.1). Here, k α = k In the following, we will need to resort to the fact that the norms on the Sobolev space H s can equivalently be stated in terms of the Fourier coefficients. The proof of the following lemma is given in [21].
Lemma 1 With computable constants C, C depending on d and s we have Lemma 1 shows that H s is identical to the space Proof The proof is indicated in [21]. In the following we work out the argument in detail. Consider an arbitrary u ∈ H s . With the Cauchy-Schwarz inequality in 2 = 2 d yields , provided that the series is convergent.
is convergent for 2s > 1, i.e., s > 1/2 = d/2. This shows that, for s > d/2, the series (B.2) is convergent and that u ∈ H s satisfies (B.1). Furthermore, the absolute summability of the Fourier coefficients c k implies that the Fourier series for u is uniformly convergent, which in turn implies the continuity of u.
Corollary 1 For s > d/2 + n we have H s ⊆ C n , and the embedding H s → C n is continuous, i.e., u C n ≤ C s,n u H s for all u ∈ H s .

B.2 Integrability
In order to study integrability properties of functions u ∈ H s we need to interrelate them to summability properties of its Fourier transform in q spaces, with û q = ∑ k∈Z d |c k | q 1 q (û = (u k ) k∈Z d ). For the proof of the following result see [22,Theorem 2.1 & 2.2] and also [23].

Remark 2
It can be shown that the assertion of Theorem 3 is also valid for the endpoint case s = d/2 and p < ∞; see [21].
Proof The proof is indicated in [21]. In the following we work out the argument in detail.
Here we have used Hölder's inequality with conjugate exponents 2 q , 2 2−q , and Lemma 1. This estimate makes sense provided the sum in the latter expression is finite, i.e., if We reason as in the proof of Theorem 2: We have for (2sq)/(2 − q) > d, i.e., q > 2d/(2s + d). With 1/p + 1/q = 1 this is equivalent to p < d/( d 2 − s), as asserted.
In the special cases d = 1, 2, 3, which are relevant to our analysis, this means: 2 − s we have H s ⊆ L p . In particular, H 1 ⊆ L 6 . 1 Here, p plays the role of q in (B.3b) and vice versa. We have 1 ≤ q ≤ 2.