Convergence analysis of high-order commutator-free quasi-Magnus exponential integrators for nonautonomous linear evolution equations of parabolic type

The main objective of this work is to provide a stability and error analysis of high-order commutator-free quasi-Magnus (CFQM) exponential integrators. These time integration methods for nonautonomous linear evolution equations are formed by products of exponentials involving linear combinations of the defining operator evaluated at certain times. In comparison with other classes of time integration methods, such as Magnus integrators, an inherent advantage of CFQM exponential integrators is that structural properties of the operator are well preserved by the arising linear combinations. Employing the analytical framework of sectorial operators in Banach spaces, evolution equations of parabolic type and dissipative quantum systems are included in the scope of applications. In this context, however, numerical experiments show that CFQM exponential integrators of nonstiff order five or higher proposed in the literature suffer from poor stability properties. The given analysis delivers insight that CFQM exponential integrators are well defined and stable only if the coefficients occurring in the linear combinations satisfy a positivity condition and that an alternative approach for the design of stable high-order schemes relies on the consideration of complex coefficients. Together with suitable local error expansions, this implies that a high-order CFQM exponential integrator retains its nonstiff order of convergence under appropriate regularity and compatibility requirements on the exact solution. Numerical examples confirm the theoretical result and illustrate the favourable behaviour of novel schemes involving complex coefficients in stability and accuracy.


Class of evolution equations
We consider evolution equations of the form u (t) = A(t) u(t), t ∈ [t 0 , T ], u(t 0 ) given (1.1) defined by a family of time-dependent linear operators (A(t)) t∈[t 0 ,T ] in a Banach space. We are primarily interested in situations where (1.1) is related to a partial differential equation of parabolic type or a dissipative quantum system, respectively. Nonautonomous linear evolution equations arise in sensitivity analysis or optimal control; further relevant applications include driven open quantum systems such as the parametrically driven dissipative Dicke model (see Alvermann et al., 2012 and references given therein).

Nonlinear equations and linearizations
As an elementary illustration, we state a one-dimensional partial differential equation comprising a nonlinear diffusion term, a nonlinear advection term, a nonlinear reaction term and an additional inhomogeneity. The time periodic logistic reaction-diffusion equation given in Pao (2001) can be cast into the above form with and constants c 0 , c 1 , c 2 , c 3 ∈ R. Associated linearized equations involve the Fréchet derivative of the second-order differential operator defining the right-hand side of the equation For instance, the variational equation, which describes the sensitivity of the solution with respect to the reference solution, corresponds to a nonautonomous linear partial differential equation ∂ t u(x, t) = α 2 (x, t) ∂ xx u(x, t) + α 1 (x, t) ∂ x u(x, t) + α 0 (x, t) u(x, t), (1.2b) with space-time-dependent coefficient functions given by . (1.2c) When rewritten as abstract differential equation, the variational equation is of the form (1.1).

Exponential integrators
A variety of contributions confirm that exponential time integration methods are favourable in various respects (see Crouch & Grossman, 1993;Hochbruck & Lubich, 2003;Blanes & Moan, 2006;Thalhammer, 2006;Alvermann & Fehske, 2011;Alvermann et al., 2012;Bader et al., 2016 and references given therein). Exponential integrators for nonautonomous evolution equations commonly rely on the computation of several exponentials, which involve the values of the defining operator at certain nodes and combinations thereof. For the associated spatially semidiscretized equation, the action of the arising matrix exponentials on vectors is often advantageously realized by polynomial approximations such as Chebyshev or Krylov methods (see also Hochbruck & Lubich, 1997;Sidje, 1998;Moler & Van Loan, 2003).

Magnus integrators
A well-established class of exponential time integration methods for nonautonomous linear evolution equations (1.1) is based on a formal solution representation by the Magnus expansion u(t n+1 ) = e Ω(τn,tn) u(t n ), t 0 ≤ t n < t n+1 = t n + τ n ≤ T , see Magnus (1954); the natural approach to truncate the infinite series and to employ quadrature approximations of the arising multiple integrals leads to interpolatory Magnus integrators. However, as has been noticed in Celledoni et al. (2003), the use of commutators may be undesirable for solving stiff systems. In particular, in the context of partial differential equations, Magnus integrators manifest a fundamental difficulty. In consideration of the fact that the commutator in general does not inherit the characteristic properties of the underlying operator, one has to face the issue of well-definedness. Moreover, as relevant applications include partial differential equations in two-and three-space dimensions, the systems resulting from spatial semidiscretization commonly involve large matrices; the computation of the action of discrete counterparts to iterated commutators, such as

Commutator-free quasi-Magnus exponential integrators
In the present work, we focus on a class of exponential time integration methods that circumvents the difficulties of Magnus integrators in the context of partial differential equations, and thus may provide a favourable alternative (see Blanes & Moan, 2006;Thalhammer, 2006;Alvermann & Fehske, 2011). In view of the attempt to avoid the presence of iterated commutators, we employ the notion commutatorfree exponential integrator, and to distinguish them from the commutator-free exponential integrators considered in Celledoni et al. (2003) and Owren (2006), we add the term quasi- Magnus. The basic idea leading to commutator-free quasi-Magnus (CFQM) exponential integrators is to replace the single exponential (1.3) by the composition of several exponentials involving linear combinations of the values of the underlying operator at certain nodes u n+1 = e τn B nJ · · · e τn B n1 ≈ u(t n+1 ) = e Ω(τn,tn) u(t n ), An inherent advantage of CFQM exponential integrators over Magnus integrators is that structural properties of the defining operator family are well preserved by the arising linear combinations. As a consequence, well-definedness and stability of the time-discrete solution can be established, for instance, within the framework of sectorial operators and analytic semigroups under natural (weak) regularity requirements on the defining operator family.

Examples
Henceforth, we denote by p ∈ N the nonstiff order of the method.
(i) As a first example, we mention the exponential midpoint rule of nonstiff order two p = 2 : u n+1 = e τn A(tn+ 1 2 τn) u n ≈ u(t n+1 ) = e Ω(τn,tn) u(t n ), (1.4) which is an instance of a Magnus integrator and likewise fits into the class of CFQM exponential integrators.
(ii) A Magnus integrator of nonstiff order four is based on two Gaussian nodes

Related classes of commutator-free exponential integrators
It is worth mentioning the connection of CFQM exponential integrators to other classes of commutatorfree exponential integrators. In Celledoni et al. (2003), autonomous nonlinear (ordinary) differential equations of the form are studied; the family of vector fields (E ) L =1 is assumed to span, at every point, the tangent space, and (f ) L =1 denotes a family of real or complex-valued functions defined on the underlying manifold. Under the assumption that for a fixed value y the differential equation is more easily solvable than the original one, explicit commutator-free exponential integrators of orders three and four involving three and five exponentials, respectively, are proposed.
A nonautonomous linear equation of the form (1.1) can be recast in this way by adding the time as a coordinate, i.e., by considering the system However, compared with CFQM exponential integrators, this approach leads to more costly schemes, requiring, for instance, five exponentials for order four; besides, we are not aware of schemes of order five or higher that have been proposed in the literature. A simpler alternative consists in considering a pth-order splitting method that is defined by real or complex coefficients (a , b ) L =1 ; for noncommutative matrices X, Y , this in particular implies The application to the autonomous system (1.6), more precisely, to the corresponding subproblems leads to the commutator-free exponential integrator Again, compared with CFQM exponential integrators, the number of exponentials grows considerably with the order; moreover, for order greater than two, schemes with real coefficients necessarily involve negative coefficients, and thus have poor stability properties for evolution equations of parabolic type.

Objective and approach
In this work, our main concern is to provide a stability and error analysis of high-order CFQM exponential integrators for the time integration of nonautonomous linear evolution equations of parabolic type. For this purpose, we employ the analytical framework of sectorial operators in Banach spaces. Preliminary numerical tests for an elementary parabolic equation showed that CFQM exponential integrators of nonstiff order six proposed in the literature suffer from poor stability properties, and a first theoretical analysis delivered insight that the structural quality of CFQM exponential integrators is only preserved under a positivity condition on certain combinations of the coefficients. Additional numerical tests confirmed a conjectured order barrier at order five for schemes involving real coefficients, and the connection to operator splitting method suggested to circumvent this order barrier by the consideration of complex coefficients.
On the basis of these findings, we deduce stability and local error estimates with respect to the norm of the underlying Banach space, which imply that a CFQM exponential integrator retains its nonstiff order of convergence, provided that the exact solution to the considered evolution equation satisfies suitable regularity and compatibility requirements. Numerical results obtained for a variational equation associated with a test equation of the form (1.2) confirm the theoretical convergence estimate and illustrate the favourable behaviour of novel schemes involving complex coefficients in stability and accuracy. The class of methods considered in this article thus has good potential for its use in the time integration of nonautonomous linear evolution equations.

Related work
Our analysis of CFQM exponential integrators extends former work on a fourth-order scheme within this class involving two exponentials; however, in Thalhammer (2006), the main objective was to explain order reductions encountered for parabolic equations under homogeneous Dirichlet boundary conditions, and thus only few expansion steps of the local error were needed. The main original contribution of the present work is the accomplishment of a suitable local error expansion, applicable to the whole class of CFQM exponential integrators. The main novel aspect of the stability analysis is to include complex coefficients and to identify the positivity condition on certain linear combinations of coefficients.
The observation that high-order CFQM exponential integrators, given in the literature involve negative coefficients, and thus suffer from poor stability properties when applied to dissipative equations of parabolic type, motivates the design of novel optimized schemes with complex coefficients. Our related work (Blanes et al., 2016) is concerned with this question; in particular, a set of independent conditions for time-symmetric methods of order six is given there. In addition, the efficiency of the novel schemes is compared with other time integration methods such as Magnus integrators for a dissipative problem in quantum mechanics.

Outline
The present manuscript is organized as follows. The class of CFQM exponential integrators is introduced in Section 2. A rigorous stability and error analysis of high-order CFQM exponential integrators is provided in Section 3. For better readability, the details of the somewhat long-winded local error expansion are included in the Appendix A; as an illustration of the general case, the expansion obtained for a fourthorder CFQM exponential integrator involving two nodes and two exponentials per time step and a Maple implementation of the resulting nonstiff order conditions are included in Appendix B. Numerical examples that confirm and complement the theoretical considerations are given in Section 4.

Notation
Let N = {n ∈ Z : n ≥ 0} denote the set of non-negative integer numbers. For a composition of noncommutative operators, we employ the convenient short notation here, by definition, the empty product is equal to the identity operator. Thoughout, C > 0 denotes a generic constant. For simplicity, we do not distinguish the solutions to a partial differential equation and to the associated abstract differential equation in notation and likewise for the defining operators.

Commutator-free quasi-Magnus exponential integrators
In this section, we introduce the general format of CFQM exponential integrators for the nonautonomous linear evolution equation (1.1) and specify higher-order schemes.

General format
As usual in a time-stepping approach, we consider suitably chosen time grid points t 0 < t 1 < · · · < t N = T and denote by the associated time increments; throughout, we employ the standard assumption that the ratios of subsequent time step sizes are bounded from below and above For a given initial approximation, the time-discrete solution values are determined by recurrence a high-order CFQM exponential integrator can be cast into the format (2.2b) As common, we relate the nodes to quadrature nodes, which we assume to be contained in the unit interval and monotonically increasing 0 ≤ c 1 < · · · < c K ≤ 1; (2.3a) for evolution equations of parabolic type, it is beneficial to permit complex coefficients in the linear combinations. Henceforth, we use the abbreviations and set d 0 = 0 as well as d 0 0 = 1. For a time-independent operator A, we have which explains the basic necessity of the positivity condition indeed, this elemental requirement ensures that a CFQM exponential integrator remains well defined within the analytical framework of sectorial operators and analytic semigroups (see also Section 3). Moreover, we tacitly assume that the consistency condition is satisfied; this relation is a direct consequence of the basic requirement S n (τ n ) = e τnA for a timeindependent operator A. We recall that p ∈ N denotes the nonstiff order of the method.

Examples
As illustration and in view of numerical comparisons described in Section 4, we recall the CFQM exponential integrators introduced before and specify the coefficients of fifth-and sixth-order schemes.
For further examples, among them optimized schemes with small effective error constants, we refer to Blanes & Moan (2006), Blanes et al. (2009), Alvermann & Fehske (2011 and Alvermann et al. (2012). Often, the coefficient matrix a ∈ C J×K is defined by the Gaussian nodes and weights (c i , w i ) K i=1 , which correspond to a quadrature approximation of maximum order 2K. In some situations, however, the use of a different quadrature formula of the same order or higher may be convenient or favourable, providing more accurate approximations without considerably increasing the computational cost. The new coefficient matrix related to quadrature nodes and weights ( c i , We note that the row sums of a are equal to the row sums of a, that is, the validity of condition (2.3c) is independent of the underlying quadrature formula.
(i) Order 2. The exponential midpoint rule (1.4) is based on a single Gaussian node and involves a single exponential at each time step evidently, condition (2.3c) is satisfied. Considering instead the trapezoidal rule, we obtain (ii) Order 4. The fourth-order CFQM exponential integrator (1.5) is based on the Gaussian quadrature formula and requires the evaluation of two exponentials at each time step For instance, for the Simpson rule, we obtain for the sixth-order Gaussian quadrature formula, we have (2.5) (iii) Order 5. A fifth-order CFQM exponential integrator with complex coefficients that satisfies condition (2.3c) is designed in Blanes et al. (2016). Employing a quadrature approximation based on the sixth-order Gaussian quadrature rule, see (2.5), the coefficient matrix reads (iv) Order 6. A nonoptimized sixth-order CFQM exponential integrator results from the coefficients (f jk ) 3 j,k=1 given in Alvermann & Fehske (2011 , Table 3, CF6:6); using again the sixth-order Gaussian quadrature rule yields (2.7) We point out that condition (2.3c) does not hold, since thus, a poor stability behaviour is observed for evolution equations of parabolic type. As an alternative, we introduce a sixth-order CFQM exponential integrator with complex coefficients (see Blanes et al., 2016); when based on the Gaussian nodes (2.5), the scheme is given by

Convergence analysis
In this section, we establish a convergence result for high-order CFQM exponential integrators applied to nonautonomous linear evolution equations of parabolic type. In Section 3.1, we specify our hypotheses on the defining operator and recapitulate basic results related to sectorial operators in Banach spaces; for details on the employed analytical framework, we refer to the standard works Tanabe (1979), Henry (1981), Pazy (1983) and Engel & Nagel (2000) and in particular to the monograph (Lunardi, 1995). An elementary initial-boundary-value problem introduced in Section 3.2 serves as illustration and test problem for numerical comparisons. A fundamental stability estimate and a suitable local error expansion are stated in Sections 3.3 and 3.4. Our main result is finally given in Section 3.5.

Analytical framework
As before, we let (X, · X ) denote the underlying Banach space and assume that the Banach space (3.1) as common, we set X 0 = X and X 1 = D. Examples for intermediate spaces are real interpolation spaces or fractional power spaces (see Henry, 1981;Pazy, 1983;Lunardi, 1995). We employ the following hypotheses on the family (A(t)) t∈[t 0 ,T ] defining the right-hand side of (1.1) (see Lunardi, 1995, Chapter 6.1).
(ii) The graph norm of A(t) and the norm in D are equivalent for any t ∈ [t 0 , T ]; that is, there exists a constant C 2 > 0 such that the relation is valid for all t ∈ [t 0 , T ]. (iii) The regularity requirement A ∈ C ϑ ([t 0 , T ], L(D, X)) holds for some ϑ ∈ (0, 1]; that is, there exists a constant C 3 > 0 such that the bound Let t ∈ [t 0 , T ]. Under the above assumptions, the sectorial operator A(t) : D → X generates an analytic semigroup e σ A(t) σ ≥0 , given by the integral formula of Cauchy for some path Γ surrounding the spectrum of A(t); the resolvent estimate (3.2a) ensures that the operator e σ A(t) : X → X remains bounded, see also relation (3.3) below.
In order to show well-definedness and stability of high-order CFQM exponential integrators, we make use of the following auxiliary considerations (see also Henry, 1981;Pazy, 1983;Lunardi, 1995).
By the required equivalence of the graph norm and the norm in D, this yields the estimate follows. (ii) Let T σ > 0. The above considerations lead to a basic estimate for the analytic semigroup (3.3a) more generally, for exponents 0 ≤ μ ≤ ν ≤ 1 and k ∈ N, we obtain see also (3.1).

Illustration
Let Ω ⊂ R denote a bounded closed interval. In regard to (1.2), we consider the initial value problem for a nonautonomous linear partial differential equation Under the basic assumption that all values of the leading space-time-dependent coefficient function are positive, the equation is of parabolic type. In order to rewrite (3.4) as an abstract Cauchy problem of the form (1.1), we choose the space of continuous functions as underlying Banach space Further restrictions apply when additional boundary conditions are imposed; for instance, in the case of homogeneous Dirichlet boundary conditions, the domain of the defining operator is given by Analogous considerations are valid for alternative choices such as X = L 2 (Ω, R). For a detailed exposition and the natural extension to several space dimensions, we refer to Lunardi (1995, Chapter 3). Provided that the boundary of the spatial domain Ω ⊂ R d is sufficiently regular and that the (complex-valued) space-time-dependent coefficient functions are sufficiently smooth with respect to the spatial variables, it is, for instance, shown that a (strongly) elliptic second-order differential operator under homogeneous boundary conditions generates a sectorial operator in the space of continuous functions (see Lunardi, 1995, Corollary 3.1.21). Furthermore, assuming the coefficient functions to be Hölder continuous with respect to time permits to extend fundamental results for the time-independent case to the time-dependent case (see Lunardi, 1995, Chapter 6). We point out that the construction of the evolution operator associated with a nonautonomous linear equation and in particular the proof of its well-definedness as operator from X to X is a nontrivial task; the stability bound given in Section 3.3 conforms to the first statement in Lunardi (1995, Corollary 6.1.8).

Well-definedness and stability
Employing the analytical framework introduced in Section 3.1, it is ensured that high-order CFQM exponential integrators remain well defined on the underlying Banach space.
Remark 3.3 Let n ∈ {0, 1, . . . , N −1} and j ∈ J . Under Hypothesis 3.1 on the linear operator family and assumption (2.3) on the method coefficients, the operator B nj : D → X defining the CFQM exponential integrator (2.2) is sectorial; indeed, employing the reformulation noting that the operator b j A(t n ) is sectorial if b j > 0 and that the remainder is bounded by the Hölder continuity assumption (3.2c), the statement follows from a perturbation result for sectorial operators (Pazy, 1983, Section 3.2).
A stability bound for high-order CFQM exponential integrators is provided by the following result.
Theorem 3.4 (Stability) Let n 0 , n ∈ {0, 1, . . . , N − 1} be such that n 0 ≤ n. Under Hypothesis 3.1, assumption (2.1) on the sequence of time step sizes, and condition (2.3c) on the method coefficients, the time-discrete evolution operator associated with a high-order CFQM exponential integrator (2.2) satisfies the following bound with a constant C > 0 that depends on t N = T , in general, but is independent of N and the step size sequence Proof. We recall that C > 0 denotes a generic constant which possibly has different values at different occurrences. Without loss of generality, we henceforth assume n > > n 0 ; otherwise, the boundedness of compositions of the time-discrete evolution operator follows at once from a repeated application of (3.3). Our proof of the above stability bound is in the lines of the preceding work (Thalhammer, 2006), see also references given therein. The basic idea is to use the relations which are valid for the analytic semigroup generated by the sectorial operator A(t), due to the consistency condition (2.3) and the bounds in (3.3); for reason that will become clear subsequently, we choose t = t n 0 + c K τ n 0 . Thus, it suffices to estimate the difference We employ the telescopic identity as a consequence, by (3.5), we obtain the estimate In a similar manner, for some i ∈ {n 0 , . . . , n}, the difference S i (τ i ) − e τ i A(t) is rewritten by means of the telescopic identity and estimated as follows X←X .
An application of the integral formula of Cauchy yields the representation and thus implies the bound see also Hypothesis 3.1 and Remarks 3.2 and 3.3. With regard to (3.6), we set t = t n 0 + c K τ n 0 such that this yields the bounds i = n 0 : Altogether by collecting the above bounds and estimating the arising Riemann sums by the corresponding integrals, we obtain A Gronwall-type inequality involving a weakly singular kernel (see, for instance, Brunner & van der Houwen, 1986) finally proves the stated result.

Local error expansion
Our objective is the derivation of a local error expansion for a high-order CFQM exponential integrator that is appropriate for nonautonomous evolution equations involving time-dependent unbounded linear operators. Compared with alternative approaches for nonstiff differential equations, it is essential to capture the remainder terms and to specify the regularity and compatibility requirements on the problem data and the exact solution. As the treatment of the general case entails certain technicalities, detailed calculations are shifted to the appendix; for the less involved case of a fourth-order CFQM exponential integrator comprising two nodes and two exponentials, the basic expansion steps are also recapitulated there.

Local error
We meanwhile fix n ∈ {0, 1, . . . , N − 1}. For the convenience of the reader, we recall useful abbreviations and a basic consistency condition on the coefficients

see (2.2) and (2.3). The local error of a CFQM exponential integrator is defined by
e τnB nj u(t n ). (3.7) Assuming that the coefficients of the method satisfy the pth nonstiff order conditions, our aim is to deduce a local error expansion which implies δ n+1 = O τ p+1 n .

Linearization and solution representation
Let j ∈ J be such that b j = 0. In order to attain a local error representation that is appropriate for further stepwise expansions, we rewrite the right-hand side of the evolution equation (1.1) as follows The variation-of-constants formula together with a linear integral transformation yields and, by induction, this leads to

Local error representation
As a consequence, we obtain the local error representation

Further expansion
Starting from the above integral relation, suitable Taylor series expansions of certain values of A and u lead to a representation of the form with certain coefficients c jk m and expressions F jk m involving the operators B ni for j ≤ i ≤ J; detailed calculations are included in the appendix. A further expansion yields the relation and, in a final step, a suitable expansion of the product involving the operators B ni for j ≤ i ≤ J is employed.

Local error expansion
Altogether this yields the local error expansion where the principal term retains the nonstiff order conditions imposed on a pth order CFQM exponential integrator (2.2); more precisely, a (redundant) set of order conditions results from inspecting the expressions associated with τ q n

Convergence result
By means of the stability bound and the local error expansion given in Sections 3.3 and 3.4, we are able to establish the following convergence result for a CFQM exponential integrator of the form (2.2). We point out once more that the local error expansion (3.10) implies that the stiff order conditions coincide with the nonstiff order conditions. A set of 10 independent order conditions for nonstiff order five and the corresponding seven conditions for time-symmetric schemes of order six are stated in our related work (Blanes et al., 2016), see also references given therein. For notational simplicity, we formulate the global error estimate for the maximum time step size On the basis of the elementary evolution equation introduced in Section 3.2, the regularity and compatibility requirements on the problem data and the exact solution are discussed below.
is valid with a constant C > 0 that does not dependent on n and τ .
Proof. As standard, we employ a telescopic identity to relate the global error to the local errors By means of Theorem 3.4, the estimate follows. Under the given assumptions on the operator family and the exact solution, the local error expansion (3.10) is well defined; in particular, the remainder is bounded in the underlying Banach space Provided that the CFQM exponential integrator (2.2) fulfils the nonstiff pth order conditions, the principal error term in (3.10) vanishes, and we obtain As a consequence, this implies the stated global error bound.
Remark 3.6 We note that any CFQM exponential integrator of order five or higher found in the literature involves merely real coefficients and does not fulfil the positivity condition in (2.3); as a consequence, the above result does not apply to it.

Elementary evolution equation
In order to substantiate the implications of Theorem 3.5, we reconsider the elementary evolution equation (3.4), choosing again the space of continuous functions as underlying function space. Provided that the coefficient functions are sufficiently often continuously differentiable, time derivatives of the defining operator again correspond to second-order differential operators; for instance, the first time derivative is given by Straighforward calculations such as

t) A(t) + A(t) A (t) + A(t) A(t) A(t) u(t)
show that the kth time derivative of the solution involves the k-fold product (A(t)) k for k ∈ N. The regularity requirements of Theorem 3.5, reflected, for instance, in hence correspond to the condition Whenever the solution is subject to boundary conditions, additional restrictions apply; for instance, in the case of homogeneous Dirichlet boundary conditions, the condition u (p−1) (t) ∈ D in particular involves the compatibility requirement

Order reductions
In situations where it is unnatural to assume that certain compositions of the defining operators and the exact solution satisfy the imposed boundary conditions, order reductions are encountered. These order reductions are explained by the fact that the relation (A(t)) k−κ u(t) ∈ X holds for certain exponents k ∈ N and κ ∈ (0, 1), but (A(t)) k u(t) ∈ X. The smoothing property of an analytic semigroup, reflected in relation (3.3), then permits to use estimates such as in order to establish global error estimates involving fractional powers of the time increments. For a detailed error analysis of a fourth-order CFQM exponential integrator and the explanation of significant order reductions, we refer to Thalhammer (2006); a numerical counterexample for an evolution equation subject to homogeneous Dirichlet boundary conditions is also given there.

Full discretization error
Our convergence result for the time-discrete solution is also a basic ingredient for the derivation of an estimate for the full discretization error. We sketch the approach for the exponential midpoint rule applied with constant time step size τ > 0 and exact initial value see also (1.4); the fully discrete solution values are obtained by a composition of the form where the index M reflects the quality of the spatial approximation. For instance, for central finite differences, the number M corresponds to the number of spatial grid points. Together with a suitable interpolation, this procedure defines an element in the underlying function space (and not just a vector comprising real numbers). The triangle inequality implies that the full discretization error is bounded by the time discretization error, see Theorem 3.5, and by the error of the additional space discretization Employing a telescopic identity, the difference u nM − u n takes the form Thus, the crucial point is to estimate products of S M (τ ) and the defect S M (τ ) − S (τ ). Under suitable regularity requirements on the data of the problem, a global error bound of the form with some exponent q > 0 is expected (see also Thalhammer, 2012 and references given therein).

Numerical examples
In this section, we present numerical comparisons for a test equation in a single-space dimension. The definition of the underlying nonlinear equation is somewhat inspired by models arising in pattern formation such as the Kuramoto-Sivashinsky equation and the Gierer-Meinhardt equations, but with a nonlinear diffusion term; for these types of problems, it is also natural to consider periodic boundary conditions. Our test equation is rather elementary from a computational point of view and has been chosen in this way for two reasons. First, in order to enable a comparison with other high-order CFQM exponential integrators given in the literature, a sufficiently low degree of freedom is required; otherwise, the poor stability behaviour of CFQM exponential integrators involving merely real coefficients would entirely spoil the numerical result (not a number). Secondly, our purpose here is to confirm the favourable stability behaviour of schemes involving complex coefficients and our theoretical global error estimate, rather than exhaustively analysing the efficiency of our novel schemes on realistic problems. Simulations for relevant applications are found in Alvermann et al. (2012); for a more detailed assessment of our novel schemes when applied to a dissipative quantum system of higher dimension, we refer to Blanes et al. (2016).

Test equation
Let Ω = [0, 1]. With regard to (1.2) and Section 3.2, we consider the nonlinear diffusion-advection reaction equation defined by the scalar functions f 2 (w) = 1 10 cos(w) + 11 10 , f 1 (w) = 1 10 w, f 0 (w) = w w − 1 2 ; (4.1b) the inhomogeneity is chosen such that the exact solution is equal to The associated variational equation involves the space-time-dependent coefficient functions where all values of the leading coefficient are positive; we, in particular, consider the initial state u(x, 0) = (sin(2 π x)) 2 x ∈ Ω.

Numerical results
For the space discretization, we choose M = 100 equidistant grid points and apply standard symmetric finite differences. The time integration is performed by different CFQM exponential integrators with time step sizes 2 − for ∈ {1, . . . , 10}; the numerical approximation obtained for the finest time step size serves as reference solution.
In Fig. 1, we display the global errors at time T = 1, obtained for the CFQM exponential integrators of nonstiff orders p = 2, 4, 5, 6 introduced in Section 2; the slopes of the lines reflect the orders of convergence. In accordance with Theorem 3.5, the nonstiff orders are retained for smaller time step sizes.
The numerical comparisons also confirm the poor stability behaviour of CFQM exponential integrators that do not verify condition (2.3c). Indeed, for the sixth-order scheme with real coefficients, a satisfactory result is only obtained for sufficiently small time step sizes; for larger values, it fails (not a number). This unstable behaviour changes for the worse when the spatial grid is refined, see Fig. 1.

Order reductions
In situations, where compositions of the involved operators applied to the exact solution values do not satisfy the imposed boundary conditions, an order reduction is expected (see Thalhammer, 2006). In order to illustrate this error behaviour, we consider (4.1) with and impose homogeneous Dirichlet boundary conditions. We determine the global temporal orders at time T = 1. In addition, we measure the global errors in the interior of the domain; that is, we only consider the 80 interior grid points and neglect the contributions of the grid points nearby the boundaries.
The obtained results, displayed in Fig. 1, show order reductions of the fourth-, fifth-and sixth-order schemes to third order, when the error is measured in the maximum norm, whereas higher convergence rates are obtained in the interior of the spatial domain.

Effect of space discretization
For comparison, we include the corresponding results obtained for M = 120 grid points in Fig. 1; due to the fact that the size of the errors are retained for the stable schemes, we conclude that the spatial error is negligible and that the numerical results are independent of the considered space discretization method. by a sectorial operator. For t ∈ [t 0 , T ], σ ≥ 0 and m ∈ N, we set referred as ϕ-functions. Evaluation at zero is understood as right-sided limit σ → 0 + , such that ϕ m (0) = 1 m! I for any m ∈ N. The corresponding result for the analytic semigroup e σ A(t) σ ≥0 ensures that ϕ m (σ A(t)) : X μ → X μ remains bounded for exponents μ ∈ [0, 1] and σ ∈ [0, T σ ], see (3.3). By means of the recurrence relation we obtain the expansion

Taylor series expansions and integration
Let (j, k) ∈ J × K . In order to expand the local error (3.8) with respect to the time increment, we employ the Taylor series expansions see (A.1); clearly, the term involving = 0 vanishes. An elementary calculation yields Inserting the above representation into (3.8) implies The integrals arising in the principal local error term are rewritten by means of the binomial series and the ϕ-functions, which leads to As a consequence, choosing q = p − 1 and r = q − 1 = p − 2, where p denotes the nonstiff order of (2.2), we obtain the local error representation in order to indicate the dependence of the remainder on the time increment and the highest time derivatives, we employ the short notation

Expansion of exponential functions
It remains to suitably expand the principal term in the local error representation (A.3). For this purpose, we apply the following relations for the ϕ-functions see (A.2). This leads to the local error expansion where the remainder is given by (A.3c) and as before, we indicate the highest time-derivatives and write for short. In a further expansion step, we repeatedly employ the recurrence relation .
we point out that the range for the index λ J−1 is chosen in dependence of the preceding index λ J such that the remainder involves the same power of the time increment, namely τ Λ J n . In order to indicate the dependence of the remainder on certain values of A and products thereof, comprising a total of Λ J factors, we henceforth employ the short notation O τ Λ J n , A Λ J ; hence, the above relation rewrites as In an inductive manner, for any j ∈ J , the expansion follows; here, we set |λ| = λ j+1 + · · · + λ J as well as λ! = λ j+1 ! · · · λ J ! for λ = (λ j+1 , . . . , λ J ) ∈ N J−j . As a consequence, we obtain the remainder R (3) n+1 comprises certain compositions of the linear operators B nj for j ∈ J and can be cast into the form O τ p+1 n , A (p) , u (p−1) . Employing the abbreviation we may rewrite δ (3) n+1 as follows

Final expansion
In a final expansion step, we employ the Taylor series expansion here, for j ∈ J , we employ the abbreviations Clearly, the summation over the subindex λ (i) does not arise whenever κ i = 0 and, in particular, the summation over the index λ does not arise whenever |κ| = 0. Altogether we obtain the local error expansion (3.10), where the remainder is of the form A (p) , u (p−1) .
We note that the operators A(t) and A (t) do not commute in general; thus, using the differential equation and inserting the relations u (t) = A(t) u(t) as well as u (t) = A (t) u(t) + A(t) A(t) u(t) does not lead to a further simplification. The requirement δ n+1 = O τ 5 n implies that the quantities C 10 , C 010 , C 20 , C 11 , C 110 , C 0010 , C 020 , C 30 , C 011 , C 21 , C 12 , involving the method coefficients a jk and c k for j = 1, 2 and k = 1, 2 should vanish; this set of (redundant) order conditions possesses a uniquely determined solution, the CFQM exponential integrator (2.4) based on two Gaussian nodes. For completeness, the order conditions are included in a simple Maple-implementation.