Higher-order exponential integrators for constrained semi-linear parabolic problems

This paper provides an introduction to exponential integrators for constrained parabolic systems. In addition, building on existing results, schemes with an expected order of convergence of three and four are established and numerically tested on parabolic problems with nonlinear dynamic boundary conditions. The simulations reinforce the subjected error behavior.


Introduction
The modeling of physical processes often leads to timedependent partial differential equations (PDEs). For the numerical solution of such problems, usually weak formulations are derived. This allows the establishment of solution theories and, eventually, the application of discretization methods, leading to the computation of approximate solutions. As the description of the temporal behavior is of particular interest, we will consider the behavior in time and space separately. We will make use of the possibility to express PDEs in a weak form as ordinary differential equations with values in Banach spaces (abstract ODEs). Related solutions are abstract functions and map from a finite time interval to a Banach space, which itself contains functions describing the spatial behavior. When modeling PDEs, which underlie additional constraints, one obtains partial differentialalgebraic equations (PDAEs), which are also called abstract differential-algebraic equations (abstract DAEs). The subject of this article are PDAEs of the forṁ with a nonlinear right-hand side f . They are called semilinear parabolic problems, as the constraints given in eq. (1b) are linear, while eq. (1a) contains a nonlinearity. In the PDAE (1), B is called the constraint operator, which couples eq. (1a) and (1b) with the help of the Lagrange multiplier λ. This article focuses on the derivation of time integration schemes by so-called exponential integrators building on the recent article [2]. These schemes provide a useful tool for the time discretization, as they allow large time steps since they are not restricted by CFL conditions. This means no fixed relations between the spatial and temporal mesh sizes have to be met in order to obtain stable solutions. The main contribution of this article is the establishment of schemes with an expected order of convergence of three and four. This article is structured as follows. In Section 2, we provide further insights into the PDAE (1) and mention the assumptions required for meaningful solutions while providing material for further studies. Then, a This is an open-access article under the CC BY license (http: // creativecommons. org/ licenses/ by/ 4. 0/ ). model example is introduced, which is used for the numerical simulations of this paper. In Section 3, the idea behind the application of exponential integrators is explained. As this does not yield straightforward methods, a practical algorithm is derived in Section 3.1. Building on this, schemes of higher order are introduced in Section 3.2, while providing the tools for the derivation of further algorithms. The results of numerical simulations of the introduced exponential integrators are then presented in Section 4. Finally, a summary and a short outlook are given in Section 5.

Preliminaries
First, we want to give more insight into the PDAE (1) and state the assumptions required to obtain meaningful solutions. We will provide the necessary definitions for the scope of this article and give references for further studies. Then, we will state an example of a PDE that can be weakly formulated as a PDAE of the form (1) and will serve as a model for our numerical simulations in Section 4.
The PDAE (1) is considered on a finite time interval [0, T ], T > 0. Therefore, we are looking for functions u : [0, T ] → V and λ : [0, T ] → Q, such that the PDAE (1) holds true with initial condition u(0) = u 0 . Furthermore, V and Q denote Hilbert spaces with their dual spaces V * and Q * , respectively. As it allows more general right-hand sides, the space V is part of an evolution triple V , H , V * , which is also known as a Gelfand triple. Figuratively speaking, this is a triple of nested spaces and for an introduction of this concept we want to refer to [9,Ch. 23.4]. As a consequence, H is the natural space for the initial data. Following [2], we will make the slightly stronger assumption u(0) = u 0 ∈ V and require the consistency condition Bu 0 = g (0). The operator A ∈ L (V , V * ) is assumed to be linear, continuous and satisfies a Gårding inequality on V ker := ker B ⊆ V , cf. [9,Ch. 23.7]. The linear and continuous operator B ∈ L (V , Q * ) has to satisfy an inf-sup condition, which implies that B has a linear, continuous right-inverse denoted by B − : Q * → V . For an introduction to infsup stability, see, for example, [3,Ch. III.4]. Finally, the right-hand sides are assumed f ∈ L 2 (0, T ; V * ) and g ∈ H 1 (0, T ; Q * ) with some additional assumptions on the nonlinear part, namely the classical Carathéodory condition and a boundedness condition, cf. [2], leading to the existence of a unique solution. By definition of the above spaces, f and g take values in V * and Q * , respectively, almost everywhere in [0, T ] and a formal introduction is given in [9,Ch. 23].
As a model example fitting into the framework of the PDAE (1), we consider the following semi-linear parabolic problem with nonlinear dynamic boundary conditionu Arising from the requirements to obtain a meaningful weak formulation, the following assumptions are made on the parameters: Note, that u and all occurring parameters apart from β may depend on the spatial variable in (2). We emphasize that both, eq. (2a) and the dynamic boundary condition eq. (2b) may include linear or nonlinear reaction terms. Following [1], the PDE (2) can be formulated weakly as a PDAE of the form which we require to hold a.e. in [0, T ] with initial condition Therein, the notation of the time dependence was omitted and p denotes an auxiliary variable modelling the dynamics on the boundary. Clearly, this PDAE is of the form (1) and therefore fits into our considered framework. In the above, V is defined by the product space H 1 (Ω)×H 1/2 (Γ dyn ) for β = 0 and by H 1 (Ω)×H 1 (Γ dyn ) for β > 0. The pivot space is given by H := L 2 (Ω) × L 2 (Γ dyn ) and the space Q by H −1/2 (Γ dyn ). Furthermore, the coupling operator is defined by B Lastly, the operator A is given by where 〈K u, v〉 := Ω (κ∇u) · ∇v dx for all u, v ∈ H 1 (Ω) and, similarly, 〈K Γ p, q〉 := Γ ∇ Γ p · ∇ Γ q dσ for all p, q ∈ H 1 (Γ) denotes the weak form of ∆ Γ . We revisit this model example in the numerical simulations in Section 4.

Exponential Integrators
The PDAE (1) allows the application of exponential integrators for semi-linear parabolic problems with linear constraints, which are introduced in [2]. As it is important for the construction of higher-order schemes, we explain the idea behind it. In the following, we assume without loss of generality that A is elliptic on V ker , cf. [3,Ch. III.4]. This is not a restriction, as we can add the term κu to A and add it to the nonlinearity f accordingly in eq. (1a). The constant κ ≥ 0 has to be chosen appropriate to the Gårding inequality leading to an A that is elliptic on V ker . The derivation presented in [2] makes use of the fact that the space V can be decomposed into As the choice of the complement space V c is not unique in general and allows a certain freedom in the modeling process, it is here chosen as Solutions of the PDAE (1) then decompose into u = u ker + u c and the latter part is fully determined by the constraint (1b), that is by u c (t ) = B − g (t ) ∈ V c . The remaining part u ker is then determined by considering the restriction of eq. (1a) to the test space V * ker . This is meaningful, as the restriction of the operator A yields an operator A ker : V ker → V * ker and the right-hand side is well-defined as a functional in V * ker by restriction of test functions to V ker ⊆ V . By the choice of the complement space V c , the term A u c vanishes and we therefore obtain the abstract ODĖ The solution to (3) can be obtained by applying the variation of constants formula, which requires to interpret the right-hand side as an element of H * ker . This means that the right-hand side is restricted to test functions in H ker , which is the closure of V ker with respect to the norm of H . This yields the solution formula For the approximation of the above formula, we introduce a partition of the time interval [0, T ] by Assuming a uniform partition with step size τ for simplicity, this leads to It remains to discretize the integral term of eq. (4) in a way that yields an explicit scheme. This is achieved by applying an appropriate quadrature rule, which is the idea behind exponential integrators. For more information on exponential integrators, we refer to [7], but also revisit them in Section 3.2 when constructing discretization schemes of higher order. The result of this procedure can be expressed by means of recursively defined ϕ-functions for j ≥ 0. The values for z = 0 are given by ϕ j (0) = 1/ j !. The choice z = −τA ker , for τ > 0, yields ϕ 0 (−τA ker ) = e −τA ker and we therefore obtain for all h ∈ H * ker , where id : V ker → V ker denotes the identity operator, cf. for instance [8,Ch. 11.1]. Note, that the invertibility of A ker follows from the well-known Lax-Milgram lemma as A is elliptic on V ker . These ϕfunctions are beneficial for computational purposes, as their action on a vector can be computed in an efficient manner via Krylov subspace methods as explained in [2, Sec. 5.1].

The Exponential Euler Scheme
For an illustration of the process of obtaining practical algorithmic methods by the procedure mentioned above, we explain this for the exponential Euler scheme using the abbreviation g n := g (t n ) in the following. In addition, u n denotes the approximation of u(t n ). Following [2], this method is obtained by simply choosing the evaluation of f s, u(s) −u c (s) in eq. (4) at the beginning of the interval as a quadrature rule, which yields the scheme u n+1 = B − g n+1 + ϕ 0 (−τA ker )(u n − B − g n ) Note that ϕ 1 (z) = 1 0 e (1−s)z ds. By using the recursion formula (6), this can be traced back to a single evaluation of ϕ 0 and a number of saddle point problems. That is, by applying eq. (6) for j = 0, the exponential Euler scheme equivalently reads u n+1 = B − g n+1 + ϕ 0 (−τA ker )(u n − B − g n − w n ) + w n . Therein, the auxiliary variable w n ∈ V ker is given by the solution of the stationary saddle point problem The evaluation B − g n can be obtained as the solution of the stationary saddle point problem and similarly for B − g n+1 and B −ġ n , guaranteeing For a complete algorithm, it remains to compute the action of ϕ 0 (−τA ker ) on a vector h = u n − B − g n − w n . The corresponding linear PDAE formulation readṡ in [t n , t n+1 ] with initial condition z(t n ) = h. This finally leads to the following algorithm, which converges with order one for sufficiently regular right-hand sides [2].

Algorithm 1 Exponential Euler scheme
Input: step size τ, consistent initial data u 0 ∈ V , right-hand sides f , g 1: for n = 0 to N − 1 do 2: compute B − g n , B − g n+1 and B −ġ n by (9) 3: compute w n by (8) 4: compute z as the solution of (10) on [t n , t n+1 ] with initial data z(t n ) = u n − B − g n − w n 5: set u n+1 = B − g n+1 + z(t n+1 ) + w n 6: end for Remark 1. In each step, Algorithm 1 only requires a single computation of the action of ϕ 0 (−τA ker ) on a vector h and the solution of four saddle point problems (three starting from the second step on). As mentioned earlier, the former can be done efficiently using Krylov subspace methods, as explained in [2, § 5.1]. The latter only requires the solution of linear systems such that no Newton iteration is necessary.

Higher-Order Schemes
In a similar fashion to the derivation of the exponential Euler scheme in the previous subsection, approximation schemes with a higher order of convergence can be obtained by using a better approximation for the integral term in eq. (4). This leads to the application of explicit exponential Runge-Kutta methods for semi-linear parabolic problems, which are more general explicit exponential integrators and have to be adjusted to fit into the constrained framework of the PDAE (1). The convergence of these methods was already analyzed up to order four in [5,6] and, therefore, the usage of the proposed schemes of these papers suggests resulting algorithms of a similar order.
In general, these methods for the unconstrained case, that is forv can be described by Butcher tableaus with the help of the in eq. (5) mentioned ϕ-functions, cf. [5].
where the k i are given recursively by
Remark 3. For the exponential Euler scheme, which was constructed in Section 3.1, the corresponding Butcher tableau reads 0 In [2], convergence of order one was proven for sufficiently regular right-hand sides. Also, a class of twostage methods was transferred to the PDAE case by using the decomposition of solutions and the recursion formula (6) similar to the derivation of the exponential Euler scheme. For a detailed derivation and depiction of the algorithm, we want to refer to the original source. The related Butcher tableau reads 0 c 2 c 2 ϕ 1,c 2 and convergence of order 3/2 was proven for the choice c 2 = 1, which requires the least amount of saddle point problems which need to be solved. Furthermore, convergence of order two was proven under stronger regularity assumptions.
For the derivation of a method with expected order of convergence three, we use the following scheme, which has order three for semi-linear parabolic problems according to [5]: As it leads to the minimal amount of saddle point problems, c 2 = 2/3 is chosen. By repeatedly applying recursion formula (6) for j = 0 and j = 1, the recursively defined k i are then given by For each 1 ≤ i ≤ 3, the auxiliary variable w i ∈ V ker is given by the solution of the stationary saddle point problem (8) with right-hand side k i and w i is given by the solution of the stationary saddle point problem Inserting the above, the related formula reads The resulting algorithm for this three-stage exponential integrator scheme is shown in Algorithm 2.

Algorithm 2 A three-stage exponential integrator
Input: step size τ, consistent initial data u 0 ∈ V , right-hand sides f , g 1: for n = 0 to N − 1 do Step 0 2: compute z as solution of (10) on [t n , t n+ 2 3 ] with initial condition z(t n ) = u n − B − g n − w 1 5: set u n+1 = B − g n+1 + z(t n+1 ) + 1 2 w 1 + 3 2 w 1 + 3 2 w 3 + 3 2 w 3 11: end for A method with expected convergence order four can be constructed similarly. Note that an explicit exponential integrator of at least five stages is required to obtain convergence order four for semi-linear parabolic problems, as pointed out in [5]. Therein, the Butcher scheme  was proposed for a method of order four and is, therefore, after being translated to the constrained PDAE case of this paper, used as a method with an expected order four in the following numerical simulations. As the resulting algorithms become more and more technical with an increasing number of steps, its detailed description is omitted here. The derivation is similar to the just presented three-stage scheme and can be taken from the code used for the numerical simulations in Section 4.

Numerical Simulation
For the numerical simulation using exponential integrators, we want to revisit the example of dynamic boundary conditions, which was presented in Section 2. Therefore, we choose the same setting as in [2], since it satisfies the assumptions needed for the convergence results proven in [2] and extend it by the higher-order methods presented in Section 3. We choose the time interval The initial condition is set to u(0) = u 0 = sin(πx) cos(5πy/2) and p 0 is chosen in a consistent manner, that is For the spatial discretization, bilinear finite elements on a uniform mesh are chosen with spatial mesh size h = 1/32. Figure 1 shows the convergence history of the discrete energy error, which is the (discrete) A -norm of the approximate solution, for the exponential Euler scheme (Algorithm 1), the second-order scheme resulting from the Butcher tableau (11) for c 2 = 1, the threestage scheme of Section 3.2 (Algorithm 2) and the fivestage scheme resulting from tableau (13). It can be seen, that all schemes seem to meet the expected order of convergence. Due to computational errors of the five-stage scheme, its order of convergence is slightly reduced for small time steps. Since this scheme is used to compute the reference solution, it is most likely responsible for the asymptotical behavior for the errors around 1 × 10 −7 . Note, that those errors occur for all of the shown methods, but the error size decreases for smaller stages of the corresponding Butcher tableau (about 1 × 10 −10 for the three-stage and 1 × 10 −13 for the two-stage scheme), and therefore this behavior only results from the five-stage scheme in Figure 1.

Summary and Outlook
Exponential integrators were introduced, and the idea how to obtain higher-order schemes was explained and shown for schemes up to order four. Numerical simulations of a model problem with dynamic boundary conditions fortify the conjectured convergence orders. It remains to prove the convergence of the depicted methods of supposed order three and four, which is likely to be possible in a similar manner to the proofs given in [2]. Furthermore, an optimization of the shown algorithms might be possible to reduce the errors for small step sizes τ, which depend on the number of stages of the related schemes. Another point of interest could be