Exponential integrators for semi-linear parabolic problems with linear constraints

TL;DR

This paper develops new exponential integrators for semi-linear parabolic problems with linear constraints, combining existing methods with saddle point solutions to ensure constraints are maintained during time discretization.

Contribution

It introduces a novel class of semi-explicit exponential integrators that handle constraints effectively, with proven convergence and demonstrated numerical performance.

Findings

Proven convergence rates for the proposed integrators.

Successful application to parabolic equations with nonlinear boundary conditions.

Enhanced stability and accuracy in constrained systems.

Abstract

This paper is devoted to the construction of exponential integrators of first and second order for the time discretization of constrained parabolic systems. For this extend, we combine well-known exponential integrators for unconstrained systems with the solution of certain saddle point problems in order to meet the constraints throughout the integration process. The result is a novel class of semi-explicit time integration schemes. We prove the expected convergence rates and illustrate the performance on two numerical examples including a parabolic equation with nonlinear dynamic boundary conditions.

Equations

$$\dot{u}(t)+\mathcal{A}u(t)=f(t,u)$$

$${\phi }_0(z):=e^z,{\phi }_{k+1}(z):=\frac{{\phi }_k(z)-{\phi }_k(0)}{z}.$$

$${\phi }_k(z)={\int }_0^1{e}^{(1-s)z}\frac{s^{k-1}}{(k-1)!}\mathrm{d}s.$$

$$\dot{u}(t)+\mathcal{A}u(t)=\sum _{k=1}^n\frac{f_k}{(k-1)!}{t}^{k-1}\in X$$

$$u(t)={\phi }_0(-t\mathcal{A})u_0+\sum _{k=1}^n{\phi }_k(-t\mathcal{A})f_k{t}^k.$$

$$\dot{u}(t)$$

$$\mathcal{B}u(t)$$

$$\text{adjustlimits}{\inf }_{q\in \mathcal{Q}\setminus \{0\}}{\sup }_{v\in \mathcal{V}\setminus \{0\}}\frac{⟨\mathcal{B}v,q⟩}{\Vert v{\Vert }_{\mathcal{V}}\Vert q{\Vert }_{\mathcal{Q}}}\ge \beta .$$

$${\mathcal{A}}_1\leftarrow {\mathcal{A}}_1+\frac{c_{{\mathcal{A}}_2}^2}{2{\mu }_{\mathcal{A}}}{\mathrm{id}}_{\mathcal{H}}\text{and}{\mathcal{A}}_2\leftarrow {\mathcal{A}}_2-\frac{c_{{\mathcal{A}}_2}^2}{2{\mu }_{\mathcal{A}}}{\mathrm{id}}_{\mathcal{H}}.$$

$$⟨{\mathcal{A}}_1{v}_{\ker },v_{\ker }⟩\ge {\mu }_{\mathcal{A}}\Vert v_{\ker }{\Vert }_{\mathcal{V}}^2-c_{{\mathcal{A}}_2}\Vert v_{\ker }{\Vert }_{\mathcal{V}}\Vert v_{\ker }{\Vert }_{\mathcal{H}}+\frac{c_{{\mathcal{A}}_2}^2}{2{\mu }_{\mathcal{A}}}\Vert v_{\ker }{\Vert }_{\mathcal{H}}^2\ge \frac{{\mu }_{\mathcal{A}}}{2}\Vert v_{\ker }{\Vert }_{\mathcal{V}}^2$$

$$\mu \Vert v_{\ker }{\Vert }_{\mathcal{V}}^2\le \Vert v_{\ker }{\Vert }_{{\mathcal{A}}_1}^2\le C\Vert v_{\ker }{\Vert }_{\mathcal{V}}^2.$$

$$\mathcal{V}={\mathcal{V}}_{\ker }\oplus {\mathcal{V}}_{\text{c}}\text{with}{\mathcal{V}}_{\ker }=\ker \mathcal{B},{\mathcal{V}}_{\text{c}}=\mathrm{im}{\mathcal{B}}^{-}.$$

$${\mathcal{V}}_{\text{c}}:=\{v\in \mathcal{V}|\mathcal{A}v\in {\mathcal{V}}_{\ker }^0\}=\{v\in \mathcal{V}|⟨\mathcal{A}v,w⟩=0\text{ for all }w\in {\mathcal{V}}_{\ker }\}.$$

$${\mathcal{A}}_{\ker }:=\mathcal{A}{|}_{{\mathcal{V}}_{\ker }}:{\mathcal{V}}_{\ker }\to {\mathcal{V}}_{\ker }^{\ast }:=({\mathcal{V}}_{\ker }{)}^{\ast }.$$

$$u$$

$$\Vert u(t)-{\mathcal{B}}^{-}g(t){\Vert }_{\mathcal{A}}^2\le \Vert u_0-{\mathcal{B}}^{-}g(0){\Vert }_{\mathcal{A}}^2+{\int }_0^t\Vert f(s)-{\mathcal{B}}^{-}\dot{g}(s){\Vert }_{\mathcal{H}}^2\mathrm{d}s.$$

$$\Vert f(t,v){\Vert }_{\mathcal{H}}\le k(t)(1+\Vert v{\Vert }_{\mathcal{V}})$$

$$\Vert f(t,v_1)-f(t,v_2){\Vert }_{\mathcal{H}}\le L\Vert v_1-v_2{\Vert }_{\mathcal{V}}$$

$$u$$

$$\dot{u}(t)$$

$$\mathcal{B}u(t)$$

$$\Vert \tilde{u}(t)-u_0{\Vert }_{\mathcal{V}}$$

$${\int }_0^t|k{|}^2\mathrm{d}s$$

$$L^2{t}_1$$

$${\int }_0^t\frac{3}{\mu }|k{|}^2(1+\Vert \tilde{u}{\Vert }_{\mathcal{V}}^2)\mathrm{d}s$$

$$D:=\big{\{}y\in C(0,t_{1};\mathcal{V})\,|\,\|y-\tilde{u}\|_{C([0,t_{1}],\mathcal{V})}\leq r/2\big{\}}$$

$$\mu \Vert (S_{t_1}{y}_i-\tilde{u})(t){\Vert }_{\mathcal{V}}^2$$

$$\displaystyle\overset{\eqref{eqn:f_bound}}{\leq}\!\int_{0}^{t}|k(s)|^{2}\big{(}1+\|y_{i}(s)-u_{0}\|_{\mathcal{V}}+\|u_{0}\|_{\mathcal{V}}\big{)}^{2}\,\mathrm{d}s$$

$$\displaystyle\overset{\hphantom{\eqref{eqn:f_bound}}}{\leq}\!\big{(}1+r+\|u_{0}\|_{\mathcal{V}}\big{)}^{2}\int_{0}^{t}|k(s)|^{2}\,\mathrm{d}s$$

$$\displaystyle\mu\,\|(S_{t_{1}}y_{1}-S_{t_{1}}y_{2})(t)\|^{2}_{\mathcal{V}}\overset{\eqref{eqn:PDAE_f_independent}}{\leq}\int_{0}^{t}\|f(s,y_{1}(s))-f(s,y_{2}(s))\|^{2}_{\mathcal{H}}\,\mathrm{d}s\overset{\eqref{eqn:f_lipschitz}}{\leq}L^{2}t_{1}\|y_{1}-y_{2}\|^{2}_{C(0,t_{1};\mathcal{V})}$$

$$\mu\,\|(u^{\star}-\widetilde{u})(t)\|^{2}_{\mathcal{V}}=\mu\,\|(S_{t_{1}}u^{\star}-\widetilde{u})(t)\|^{2}_{\mathcal{V}}\leq\int_{0}^{t}|k(s)|^{2}\big{(}1+\|(u^{\star}-\tilde{u})(s)\|_{\mathcal{V}}+\|\tilde{u}(s)\|_{\mathcal{V}}\big{)}^{2}\,\mathrm{d}s.$$

$$\|(u^{\star}-\widetilde{u})(t)\|^{2}_{\mathcal{V}}\leq\int_{0}^{t}\tfrac{3}{\mu}\,|k(s)|^{2}\big{(}1+\|\tilde{u}(s)\|^{2}_{\mathcal{V}}\big{)}\,\mathrm{d}s\cdot\exp\Big{(}\tfrac{3}{\mu}\int_{0}^{t}|k(s)|^{2}\,\mathrm{d}s\Big{)}$$

$${\dot{u}}_{\text{ker}}+{\mathcal{A}}_{\ker }{u}_{\text{ker}}={\dot{u}}_{\text{ker}}+\mathcal{A}{u}_{\text{ker}}=f(t,u_{\text{ker}}+u_{\text{c}})-{\dot{u}}_{\text{c}}\text{in }{\mathcal{V}}_{\ker }^{\ast }.$$

$${\iota }_0:\mathcal{H}\equiv {\mathcal{H}}^{\ast }\to {\mathcal{H}}_{\ker }^{\ast }\equiv {\mathcal{H}}_{\ker }.$$

Full text

C. Zimmer acknowledges support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) within the SFB 910, project number 163436311.

**Key words. PDAE, exponential integrator, parabolic equations, time discretization

AMS subject classifications. 65M12, 65J15, 65L80**

1. Introduction

Exponential integrators provide a powerful tool for the time integration of stiff ordinary differential equations as well as parabolic partial differential equations (PDE), cf. [Cer60, Law67, HO10]. Such integrators are based on the possibility to solve the linear part – which is responsible for the stiffness of the system – in an exact manner. As a result, large time steps are possible which makes the method well-suited for time stepping, especially for parabolic systems where CFL conditions may be very restrictive. For semi-linear ODEs and parabolic PDEs exponential integrators are well-studied in the literature. This includes explicit and implicit exponential Runge-Kutta methods [CM02, HO05a, HO05b], exponential Runge-Kutta methods of high order [LO14], exponential Rosenbrock-type methods [HOS09], and multistep exponential integrators [CP06].

In this paper, we construct and analyze exponential integrators for parabolic PDEs which underlie an additional (linear) constraint. This means that we aim to approximate the solution to

[TABLE]

which at the same time satisfies a constraint of the form $\mathcal{B}u(t)=g(t)$. Such systems can be considered as differential-algebraic equations (DAEs) in Banach spaces, also called partial differential-algebraic equations (PDAEs), cf. [EM13, LMT13, Alt15]. PDAEs of parabolic type include the transient Stokes problem (where $\mathcal{B}$ equals the divergence operator) as well as problems with nontrivial boundary conditions (with $\mathcal{B}$ being the trace operator). On the other hand, PDAEs of hyperbolic type appear in the modeling of gas and water networks [JT14, EKLS*+*18, AZ18b] and in elastic multibody modeling [Sim98].

To the best of our knowledge, exponential integrators have not been considered for PDAEs so far. In the finite-dimensional case, however, exponential integrators have been analyzed for DAEs of (differential) index 1 [HLS98]. We emphasize that the parabolic PDAEs within this paper generalize index-2 DAEs in the sense that a standard spatial discretization by finite elements leads to DAEs of index 2. Known time stepping methods for the here considered parabolic PDAEs include splitting methods [AO17], algebraically stable Runge-Kutta methods [AZ18a], and discontinuous Galerkin methods [VR18].

In the first part of the paper we discuss the existence and uniqueness of solutions for semi-linear PDAEs of parabolic type with linear constraints. Afterwards, we propose two exponential integrators of first and second order for such systems. The construction of this novel class of time integration schemes benefits from the interplay of well-known time integration schemes for unconstrained systems and stationary saddle point problems in order to meet the constraints. Since the latter is done in an implicit manner, the combination with explicit schemes for the dynamical part of the system leads to so-called semi-explicit time integration schemes. As exponential integrators are based on the exact evaluation of semigroups, we need to extend this to the constrained case. The proper equivalent is the solution of a homogeneous but transient saddle point problem, which is a linear PDAE.

The resulting exponential Euler scheme requires the solution of three stationary and a single transient saddle point problem in each time step. All these systems are linear, require in total only one evaluation of the nonlinear function, and do not call for another linearization step. Further, the transient system is homogeneous such that it can be solved without an additional regularization (or index reduction in the finite-dimensional case). The corresponding second-order scheme requires the solution of additional saddle point problems. Nevertheless, all these systems are linear and easy to solve. In a similar manner – but under additional regularity assumptions – one may translate more general exponential Runge-Kutta schemes to the constrained case. Here, however, we restrict ourselves to schemes of first and second order.

The paper is organized as follows. In Section 2 we recall the most important properties of exponential integrators for parabolic problems in the unconstrained case. Further, we introduce the here considered parabolic PDAEs, summarize all needed assumptions, and analyze the existence of solutions. The exponential Euler method is then subject of Section 3. Here we discuss two approaches to tackle the occurrence of constraints and prove first-order convergence. An exponential integrator of second order is then introduced and analyzed in Section 4. Depending on the nonlinearity, this scheme converges with order $2$ or reduced order $\nicefrac{{3}}{{2}}$. Comments on the efficient computation and numerical experiments for semi-linear parabolic systems illustrating the obtained convergence results are presented in Section 5.

2. Preliminaries

In this preliminary section we recall basic properties of exponential integrators when applied to PDEs of parabolic type. For this (and the later analysis) we consider the well-known $\varphi$ functions. Further, we introduce the precise setting for the here considered parabolic systems with constraints and discuss their solvability.

2.1. Exponential integrators for parabolic problems

As exponential integrators are based on the exact solution of linear homogeneous problems, we consider the recursively defined $\varphi$-functions, see, e.g. [SWP12, Ch. 11.1],

[TABLE]

For $z=0$ the values are given by $\varphi_{k}(0)=1/k!$. The importance of the $\varphi$-functions comes from the fact that they can be equivalently written as integrals of certain exponentials. More precisely, we have for $k\geq 1$ that

[TABLE]

We will consider these functions in combination with differential operators. For a bounded and invertible operator $\mathcal{A}\colon X\to X$ where $e^{t\mathcal{A}}:=\exp(t\mathcal{A})$ is well-defined, we can directly use the formulae in (2.1) using the notion $\frac{1}{\mathcal{A}}=\mathcal{A}^{-1}$. As a result, the exact solution of a linear abstract ODE with polynomial right-hand side can be expressed in terms of $\varphi_{k}$. More precisely, the solution of

[TABLE]

with initial condition $u(0)=u_{0}$ and coefficients $f_{k}\in X$ is given by

[TABLE]

If $-\mathcal{A}\colon D(\mathcal{A})\subset X\to X$ is an unbounded differential operator which generates a strongly continuous semigroup, then we obtain the following major property for the corresponding $\varphi$-functions.

Theorem 2.1** (cf. [HO10, Lem. 2.4]).**

Assume that the linear operator $-\mathcal{A}$ is the infinitesimal generator of a strongly continuous semigroup $e^{-t\mathcal{A}}$. Then, the operators $\varphi_{k}(-\tau\mathcal{A})$ are linear and bounded in $X$.

With the interpretation of the exponential as the corresponding semigroup, the solution formula for bounded operators (2.3) directly translates to linear parabolic PDEs of the form (2.2) with an unbounded differential operator $\mathcal{A}$, cf. [HO10].

The construction of exponential integrators for $\dot{u}(t)+\mathcal{A}u(t)=f(t,u)$ is now based on the replacement of the nonlinearity $f$ by a polynomial and (2.3). Considering the interpolation polynomial of degree 0, i.e., evaluating the nonlinearity only in the starting value of $u$, we obtain the exponential Euler scheme. The corresponding scheme for constrained systems is discussed in Section 3 and a second-order scheme in Section 4.

2.2. Parabolic problems with constraints

In this subsection, we introduce the constrained parabolic systems of interest and gather assumptions on the involved operators. Throughout this paper we consider semi-explicit and semi-linear systems meaning that the constraints are linear and that the nonlinearity only appears in the low-order terms of the dynamic equation. Thus, we consider the following parabolic PDAE: find $u\colon[0,T]\to\mathcal{V}$ and $\lambda\colon[0,T]\to\mathcal{Q}$ such that

[TABLE]

Therein, $\mathcal{V}$ and $\mathcal{Q}$ denote Hilbert spaces with respective duals $\mathcal{V}^{*}$ and $\mathcal{Q}^{*}$. The space $\mathcal{V}$ is part of a Gelfand triple $\mathcal{V},\mathcal{H},\mathcal{V}^{*}$, cf. [Zei90, Ch. 23.4]. This means that $\mathcal{V}$ is continuously (and densely) embedded in the pivot space $\mathcal{H}$ which implies $\mathcal{H}^{*}\hookrightarrow\mathcal{V}^{*}$, i.e., the continuous embedding of the corresponding dual spaces. In this setting, the Hilbert space $\mathcal{H}$ is the natural space for the initial data. Note, however, that the initial condition may underlie a consistency condition due to the constraint (2.4b), cf. [EM13]. For the here considered analysis we assume slightly more regularity, namely $u(0)=u_{0}\in\mathcal{V}$, and consistency of the form $\mathcal{B}u_{0}=g(0)$.

The assumptions on the operators $\mathcal{A}\in\mathcal{L}(\mathcal{V},\mathcal{V}^{*})$ and $\mathcal{B}\in\mathcal{L}(\mathcal{V},\mathcal{Q}^{*})$ are summarized in the following.

Assumption 2.2* (constraint operator $\mathcal{B}$).*

The operator $\mathcal{B}\colon\mathcal{V}\to\mathcal{Q}^{*}$ is linear, continuous, and satisfies an inf-sup condition, i.e., there exists a constant $\beta>0$ such that

[TABLE]

Assumption 2.3* (differential operator $\mathcal{A}$).*

The linear and continuous operator $\mathcal{A}\colon\mathcal{V}\to\mathcal{V}^{*}$ has the form $\mathcal{A}=\mathcal{A}_{1}+\mathcal{A}_{2}$ with $\mathcal{A}_{1}\in\mathcal{L}(\mathcal{V},\mathcal{V}^{*})$ being self-adjoint and $\mathcal{A}_{2}\in\mathcal{L}(\mathcal{V},\mathcal{H})$. Further, we assume that $\mathcal{A}$ is elliptic on $\mathcal{V}_{\ker}:=\ker\mathcal{B}$, i.e., on the kernel of the constraint operator.

Without loss of generality, we may assume under Assumption 2.3 that $\mathcal{A}_{1}$ is elliptic on $\mathcal{V}_{\ker}$. This can be seen as follows: With $\mu_{\mathcal{A}}$ denoting the ellipticity constant of $\mathcal{A}$ and $c_{\mathcal{A}_{2}}$ the continuity constant of $\mathcal{A}_{2}$, we set

[TABLE]

This then implies

[TABLE]

for all $v_{\ker}\in\mathcal{V}_{\ker}$. Hence, we assume throughout this paper that, given Assumption 2.3, $\mathcal{A}_{1}$ is elliptic on $\mathcal{V}_{\ker}$. As a result, $\mathcal{A}_{1}$ induces a norm which is equivalent to the $\mathcal{V}$-norm on $\mathcal{V}_{\ker}$, i.e.,

[TABLE]

Remark 2.4*.*

The results of this paper can be extended to the case where $\mathcal{A}$ only satisfies a Gårding inequality on $\mathcal{V}_{\ker}$. In this case, we add to $\mathcal{A}$ the term $\kappa\,\operatorname{id}_{\mathcal{H}}$ such that $\mathcal{A}+\kappa\,\operatorname{id}_{\mathcal{H}}$ is elliptic on $\mathcal{V}_{\ker}$ and add it accordingly to the nonlinearity $f$.

Assumption 2.2 implies that $\mathcal{B}$ is onto such that there exists a right-inverse denoted by $\mathcal{B}^{-}\colon\mathcal{Q}^{*}\to\mathcal{V}$. This in turn motivates the decomposition

[TABLE]

We emphasize that the choice of the right-inverse (and respectively $\mathcal{V}_{\text{c}}$) is, in general, not unique and allows a certain freedom in the modeling process. Within this paper, we define the complementary space $\mathcal{V}_{\text{c}}$ as in [AZ18a] in terms of the annihilator of $\mathcal{V}_{\ker}$, i.e.,

[TABLE]

The analysis of constrained systems such as (2.4) is heavily based on the mentioned decomposition of $\mathcal{V}$. Furthermore, we need the restriction of the differential operator to the kernel of $\mathcal{B}$, i.e.,

[TABLE]

Note that we use here the fact that functionals in $\mathcal{V}^{*}$ define functionals in $\mathcal{V}_{\ker}^{*}$ simply through the restriction to $\mathcal{V}_{\ker}$. The closure of $\mathcal{V}_{\ker}$ in the $\mathcal{H}$-norm is denoted by $\mathcal{H}_{\ker}:=\overline{\mathcal{V}_{\ker}}^{\mathcal{H}}$. Assumption 2.3 now states that $\mathcal{A}_{\ker}$ is an elliptic operator. This in turn implies that $-\mathcal{A}_{\ker}$ generates an analytic semigroup on $\mathcal{H}_{\ker}$, see [Paz83, Ch. 7, Th. 2.7].

Finally, we need assumptions on the nonlinearity $f$. Here, we require certain smoothness properties such as local Lipschitz continuity in the second component. The precise assumptions will be given in the respective theorems.

Example 2.5**.**

The (weak) formulation of semi-linear parabolic equations with dynamical (or Wentzell) boundary conditions [SW10] fit into the given framework. For this, the system needs to be formulated as a coupled system which leads to the PDAE structure (2.4), cf. [Alt19]. We emphasize that also the boundary condition may include nonlinear reaction terms. We will consider this example in the numerical experiments of Section 5.

2.3. Existence of solutions

In this section we discuss the existence of solutions to (2.4), where we use the notion of Sobolev-Bochner spaces $L^{2}(0,T;X)$ and $H^{1}(0,T;X)$ for a Banach space $X$, cf. [Zei90, Ch. 23]. For the case that $f$ is independent of $u$, the existence of solutions is well-studied, see [Tar06, EM13, Alt15]. We recall the corresponding result in the special of $\mathcal{A}$ being self-adjoint, which is needed in later proofs.

Lemma 2.6**.**

Let $\mathcal{A}\in\mathcal{L}(\mathcal{V},\mathcal{V}^{\ast})$ be self-adoint and elliptic on $\mathcal{V}_{\ker}$ and let $\mathcal{B}$ satisfy Assumption 2.2. Further, assume $f\in L^{2}(0,T;\mathcal{H})$, $g\in H^{1}(0,T;\mathcal{Q}^{\ast})$, and $u_{0}\in\mathcal{V}$ with $\mathcal{B}u_{0}=g(0)$. Then, the PDAE (2.4) with right-hand sides $f$, $g$ – independent of $u$ – and initial value $u_{0}$ has a unique solution

[TABLE]

with $u(0)=u^{0}$. The solution depends continuously on the data and satisfies

[TABLE]

Proof.

A sketch of the proof can be found in [Tar06, Lem. 21.1]. For more details we refer to [Zim15, Ch. 3.1.2.2]. ∎

In order to transfer the results of Lemma 2.6 to the semi-linear PDAE (2.4) we need to reinterpret the nonlinearity $f\colon[0,T]\times\mathcal{V}\to\mathcal{H}$ as a function which maps an abstract measurable function $u\colon[0,T]\to\mathcal{V}$ to $f(\,\cdot\,,u(\,\cdot\,))\colon[0,T]\to\mathcal{H}$. For this, we need the classical Carathéodory condition, see [GKT92, Rem. 1], i.e.,

1. i.)

$v\mapsto f(t,v)$ is a continuous function for almost all $t\in[0,T]$, 2. ii.)

$t\mapsto f(t,v)$ is a measurable function for all $v\in\mathcal{V}$.

Furthermore, we need a boundedness condition such that the Nemytskii map induced by $f$ maps $C([0,T];\mathcal{V})$ to $L^{2}(0,T;\mathcal{H})$. We will assume in the following that there exists a function $k\in L^{2}(0,T)$ such that

[TABLE]

for all $v\in\mathcal{V}$ and almost all $t\in[0,T]$. We emphasize that condition (2.7) is sufficient but not necessary for $f$ to induce a Nemytskii map, cf. [GKT92, Th. 1(ii)]. We will use this condition to prove the existence and uniqueness of a global solution to (2.4).

Theorem 2.7**.**

Assume that $\mathcal{A}$ and $\mathcal{B}$ satisfy Assumptions 2.2 and 2.3. Let $g\in H^{1}(0,T;\mathcal{Q}^{\ast})$ and suppose that $f\colon[0,T]\times\mathcal{V}\to\mathcal{H}$ satisfies the Carathéodory conditions as well as the uniform bound (2.7). In addition, for every $v\in\mathcal{V}$ there exists an open ball $B_{r}(v)\subseteq\mathcal{V}$ with radius $r=r(v)>0$ and a constant $L=L(v)\geq 0$, such that for almost every $t\in[0,T]$ it holds that

[TABLE]

for all $v_{1},v_{2}\in B_{r}(v)$. Then, for a consistent initial value $u_{0}\in\mathcal{V}$, i.e., $\mathcal{B}u_{0}=g(0)$, the semi-linear PDAE (2.4) has a unique solution

[TABLE]

with $u(0)=u^{0}$.

Proof.

Without loss of generality, we assume that $\mathcal{A}=\mathcal{A}_{1}$. For this, we redefine $f(t,v)\leftarrow f(t,v)-\mathcal{A}_{2}v$, leading to an update of the involved constants $L\leftarrow L+c_{\mathcal{A}_{2}}$ and $k\leftarrow k+c_{\mathcal{A}_{2}}$ but leaving the radius $r$ unchanged.

To prove the statement we follow the steps of [Paz83, Ch. 6.3]. Let $t^{\prime}\in(0,T]$ be arbitrary but fixed. With (2.7) we notice that the Nemyskii map induced by $f$ maps $C(0,t^{\prime};\mathcal{V})$ into $L^{2}(0,t^{\prime};\mathcal{H})$, cf. [GKT92, Th. 1]. Therefore, the solution map $S_{t^{\prime}}\colon C(0,t^{\prime};\mathcal{V})\to C(0,t^{\prime};\mathcal{V})$, which maps $y\in C(0,t^{\prime};\mathcal{V})$ to the solution of

[TABLE]

with initial value $u_{0}$, is well-defined, cf. Lemma 2.6. To find a solution to (2.4) we have to look for a unique fixed point of $S_{t^{\prime}}$ and show that $t^{\prime}$ can be extended to $T$.

Let $\tilde{u}\in C(0,T;\mathcal{V})$ be the solution of the PDAE (2.4) for $f\equiv 0$ and initial value $u_{0}$. With $r=r(u_{0})$ and $L=L(u_{0})$ we now choose $t_{1}\in(0,T]$ such that

[TABLE]

for all $t\in[0,t_{1}]$. This is well-defined, since $\tilde{u}-u_{0}$ and the integrals in (b) and (d) are continuous functions in $t$, which vanish for $t=0$. We define

[TABLE]

and consider $y_{1},y_{2}\in D$. By (a) we have $\|y_{i}-u_{0}\|_{C([0,t_{1}],\mathcal{V})}\leq r$. Using that $\tilde{u}$ and $S_{t_{1}}y_{i}$ satisfy the constraint (2.9b), we obtain the estimate

[TABLE]

which implies with (b) that $S_{t_{1}}$ maps $D$ into itself. Further, we have

[TABLE]

for all $t\leq t_{1}$, $i=1,2$. Together with the previous estimate and (c), this shows that $S_{t_{1}}$ is a contraction on $D$. Hence, there exists a unique fixed point $u\in D\subset C([0,t_{1}],\mathcal{V})$ of $S_{t_{1}}$ by the Banach fixed point theorem [Zei92, Th. 1.A]. On the other hand, for every fixed point $u^{\star}=S_{t_{1}}u^{\star}$ in $C([0,t_{1}],\mathcal{V})$, we have the estimate

[TABLE]

Using $(a+b+c)^{2}\leq 3\,(a^{2}+b^{2}+c^{2})$ and Gronwall’s inequality it follows that

[TABLE]

for every $t\leq t_{1}$. Because of (d), this shows that $u^{\star}$ is an element of $D$ and thus, $u^{\star}=u$.

By considering problem (2.4) iterativly from $[t_{i-1},T]$, $t_{0}:=0$, to $[t_{i},T]$ with consistent initial value $u_{0}=u(t_{i})$, we can extend $u$ uniquely on an interval $\mathbb{I}$ with $u\in C(\mathbb{I};\mathcal{V})$ and $u=S_{t^{\prime}}u$ for every $t^{\prime}\in\mathbb{I}$. Note that either $\mathbb{I}=[0,T]$ or $\mathbb{I}=[0,T^{\prime})$ with $T^{\prime}\leq T$. The second case is only possible if $\|u(t)\|_{\mathcal{V}}\to\infty$ for $t\to T^{\prime}$, otherwise we can extend $u$ by starting at $T^{\prime}$. But, since the estimate (2.10) also holds for $u=u^{\star}$ and $t<T^{\prime}$, we have in limit that $\|u(T^{\prime})\|_{\mathcal{V}}\leq\|u(T^{\prime})-\tilde{u}(T^{\prime})\|_{\mathcal{V}}+\|\tilde{u}(T^{\prime})\|_{\mathcal{V}}$ is bounded. Therefore, $u=S_{T}u\in C([0,T];\mathcal{V})$. Finally, the stated spaces for $u$ and $\lambda$ follow by Lemma 2.6 with right-hand side $f=f(\,\cdot\,,u(\,\cdot\,)).$ ∎

Remark 2.8*.*

Under the given assumptions on $f$ from Theorem 2.7, one finds a radius $r_{u}>0$ and a Lipschitz constant $L_{u}\in[0,\infty)$, both based on the solution $u$, such that (2.8) holds for all $x,y\in B_{r_{u}}(u(s))$ with $L=L_{u}$ and arbitrary $s\in[0,T].$ With these uniform constants one can show that the mapping of the data $u_{0}\in\mathcal{V}$ and $g\in H^{1}(0,T;\mathcal{Q})$ with $\mathcal{B}u_{0}=g(0)$ to the solution $(u,\lambda)$ is continuous.

Remark 2.9*.*

It is possible to weaken the assumption (2.7) of Theorem 2.7 to $\|f(t,v)\|_{\mathcal{H}}\leq k(t)(1+\|v\|^{p}_{\mathcal{V}})$ for an arbitrary $p>1$. Under this assumption one can show the existence of a unique solution of (2.4), which may only exists locally.

Remark 2.10*.*

The assumptions considered in [Paz83, Ch. 6.3] are stronger then the one in Theorem 2.7. If these additional assumptions are satisfied, then the existence and uniqueness of a solution to (2.4) follows directly by Lemma 2.6, [Paz83, Ch. 6, Th. 3.1 & 3.3], and the fact that every self-adjoint, elliptic operator $\mathcal{A}\in\mathcal{L}(\mathcal{V},\mathcal{V}^{\ast})$ has a unique invertible square root $\mathcal{A}^{\nicefrac{{1}}{{2}}}\in\mathcal{L}(\mathcal{V},\mathcal{H})$ with $\langle\mathcal{A}v_{1},v_{2}\rangle=(\mathcal{A}^{\nicefrac{{1}}{{2}}}v_{1},\mathcal{A}^{\nicefrac{{1}}{{2}}}v_{2})$ for all $v_{1},v_{2}\in\mathcal{V}$. This can be proven by interpreting $\mathcal{A}$ as an (unbounded) operator $\mathbf{A}\colon D(\mathbf{A})\subset\mathcal{H}\to\mathcal{H}$ with domain $D(\mathbf{A}):=\mathcal{A}^{-1}\mathcal{H}\subset\mathcal{V}\hookrightarrow\mathcal{H}$ and the results of [BS87, Ch. 6, Th. 4 & Ch. 10, Th. 1] and [Paz83, Th. 6.8].

2.4. A solution formula for the linear case

In the linear case, the solution of (2.4) can be expressed by the variation-of-constants formula (Duhamel’s principle), cf. [EM13]. In the semi-linear case, we consider the term $f(t,u)$ as a right-hand side which leads to an implicit formula only. This, however, is still of value for the numerical analysis of time integration schemes.

The solution formula is based on the decomposition $u=u_{\text{ker}}+u_{\text{c}}$ with $u_{\text{ker}}\colon[0,T]\to\mathcal{V}_{\ker}$ and $u_{\text{c}}\colon[0,T]\to\mathcal{V}_{\text{c}}$. The latter is fully determined by the constraint (2.4b), namely $u_{\text{c}}(t)=\mathcal{B}^{-}g(t)\in\mathcal{V}_{\text{c}}$. For $u_{\text{ker}}$ we consider the restriction of (2.4a) to the test space $\mathcal{V}_{\ker}$. Since the Lagrange multiplier disappears in this case, we obtain

[TABLE]

Note that the right-hand side is well-defined as functional in $\mathcal{V}_{\ker}^{*}$ using the trivial restriction of $\mathcal{V}^{*}$ to $\mathcal{V}_{\ker}^{*}$. Further, the term $\mathcal{A}u_{\text{c}}$ disappears under test functions in $\mathcal{V}_{\ker}$ due to the definition of $\mathcal{V}_{\text{c}}$. If this orthogonality is not respected within the implementation, then this term needs to be reconsidered.

The solution to (2.11) can be obtained by an application of the variation-of-constants formula. Since the semigroup can only be applied to functions in $\mathcal{H}_{\ker}$, we introduce the operator

[TABLE]

This operator is again based on a simple restriction of test functions and leads to the solution formula

[TABLE]

Assuming a partition of the time interval $[0,T]$ by $0=t_{0}<t_{1}<\dots<t_{N}=T$, we can write the solution formula in the form

[TABLE]

Note that we use here the abbreviation $g_{n}:=g(t_{n})$. In the following two sections we construct exponential integrators for constrained semi-linear systems of the form (2.4). Starting point is a first-order scheme based on the exponential Euler method applied to equation (2.11).

3. The Exponential Euler Scheme

The idea of exponential integrators is to approximate the integral term in (2.12) by an appropriate quadrature rule. Following the construction for PDEs [HO10], we consider in this section the function evaluation at the beginning of the interval. This then leads to the scheme

[TABLE]

As usual, $u_{n}$ denotes the approximation of $u(t_{n})$. Further, we restrict ourselves to a uniform partition of $[0,T]$ with step size $\tau$ for simplicity. Assuming that the resulting approximation satisfies the constraint in every step, we have $u_{n}-\mathcal{B}^{-}g_{n}\in\mathcal{V}_{\ker}\hookrightarrow\mathcal{H}_{\ker}$ such that the semigroup $e^{-\tau\mathcal{A}_{\ker}}$ is applicable. The derived formula (3.1) is beneficial for the numerical analysis but lacks the practical access which we tackle in the following.

3.1. The practical method

Since the evaluation of the $\varphi$-functions with the operator $\mathcal{A}_{\ker}$ is not straightforward, we reformulate the method by a number of saddle point problems. Furthermore, we need evaluations of $\mathcal{B}^{-}$ applied to the right-hand side $g$ (or its time derivative). Also this is replaced by the solution of a saddle point problem.

Consider $x:=\mathcal{B}^{-}g_{n}=\mathcal{B}^{-}g(t_{n})\in\mathcal{V}_{\text{c}}\subseteq\mathcal{V}$. Then, $x$ can be written as the solution of the stationary auxiliary problem

[TABLE]

Note that equation (3.2b) enforces the connection of $x$ to the right-hand side $g$ whereas the first equation of the system guarantees the desired $\mathcal{A}$-orthogonality. The Lagrange multiplier $\nu$ is not of particular interest and simply serves as a dummy variable. The unique solvability of system (3.2) is discussed in the following lemma.

Lemma 3.1**.**

Let the operators $\mathcal{A}$ and $\mathcal{B}$ satisfy Assumptions 2.2 and  2.3. Then, for every $g_{n}\in\mathcal{Q}^{*}$ there exists a unique solution $(x,\nu)\in\mathcal{V}_{\text{c}}\times\mathcal{Q}$ to system (3.2).

Proof.

Under the given assumptions on the operators $\mathcal{A}$ and $\mathcal{B}$ there exists a unique solution $(x,\nu)\in\mathcal{V}\times\mathcal{Q}$ to (3.2), even in the case with an inhomogeneity in the first equation, see [BF91, Ch. II, Prop. 1.3]. It remains to show that $x$ is en element of $\mathcal{V}_{\text{c}}$. For this, note that $x$ satisfies $\langle\mathcal{A}x,w\rangle=0$ for all $w\in\mathcal{V}_{\ker}$, since the $\mathcal{B}^{*}$-term vanishes for these test functions. This, however, is exactly the definition of the complement space $\mathcal{V}_{\text{c}}$. ∎

Being able to compute $\mathcal{B}^{-}g_{n}$, we are now interested in the solution of problems involving the operator $\mathcal{A}_{\ker}$. This will be helpful for the reformulation of the exponential Euler method (3.1). We introduce the auxiliary variable $w_{n}\in\mathcal{V}_{\ker}$ as the solution of

[TABLE]

This is again equivalent to a stationary saddle point problem, namely

[TABLE]

As above, the Lagrange multiplier is only introduced for a proper formulation and not of particular interest in the following. The unique solvability of system (3.3) follows again by Lemma 3.1, since the right-hand side of the first equation is an element of $\mathcal{V}^{*}$. In order to rewrite (3.1), we further note that the recursion formula for $\varphi_{1}$ implies

[TABLE]

for all $h\in\mathcal{H}_{\ker}$. Recall that $\mathcal{A}_{\ker}$ is indeed invertible due to Assumption 2.3. Thus, the exponential Euler scheme can be rewritten as

[TABLE]

Finally, we need a way to compute the action of $\varphi_{0}(-\tau\mathcal{A}_{\ker})$. For this, we consider the corresponding PDAE formulation. The resulting method then reads $u_{n+1}=\mathcal{B}^{-}g_{n+1}+z(t_{n+1})+w_{n}$, where $z$ is the solution of the linear homogeneous PDAE

[TABLE]

with initial condition $z(t_{n})=u_{n}-\mathcal{B}^{-}g_{n}-w_{n}$. Thus, the exponential Euler scheme given in (3.1) can be computed by a number of saddle point problems. We summarize the necessary steps in Algorithm 1.

Remark 3.2*.*

One step of the exponential Euler scheme consists of the solution of four (from the second step on only three) stationary and a single transient saddle point problem, including only one evaluation of the nonlinear function $f$. We emphasize that all these systems are linear such that no Newton iteration is necessary in the solution process. Furthermore, the time-dependent system is homogeneous such that it can be solved without the need of a regularization.

3.2. Convergence analysis

In this section we analyze the convergence order of the exponential Euler method for constrained PDAEs of parabolic type. For the unconstrained case it is well-known that the convergence order is one. Since our approach is based on the unconstrained PDE (2.11) of the dynamical part in $\mathcal{V}_{\ker}$, we expect the same order for the solution of Algorithm 1. For the associated proof we will assume that the approximation $u_{n}$ lies within a strip of radius $r$ around $u$, where $f$ is locally Lipschitz continuous with constant $L>0$. Note that by Remark 2.8 there exists such a uniform radius and local Lipschitz constant. Furthermore, a sufficiently small step size $\tau$ guarantees that $u_{n}$ stays within this strip around $u$, since the solution $z$ of (3.4) and $\mathcal{B}^{-}g$ are continuous.

Theorem 3.3** (Exponential Euler).**

Consider the assumptions of Theorem 2.7 including Assumptions 2.2 and 2.3. Further, let the step size $\tau$ be sufficiently small such that the derived approximation $u_{n}$ lies within a strip along $u$ in which $f$ is locally Lipschitz continuous with a uniform constant $L>0$. For the right-hand side of the constraint we assume $g\in H^{2}(0,T;\mathcal{Q}^{*})$. If the exact solution of (2.4) satisfies $\frac{\mathrm{d}}{\mathrm{d}t}f(\cdot,u(\cdot))\in L^{2}(0,T;\mathcal{H})$, then the approximation $u_{n}$ obtained by the exponential Euler scheme of Algorithm 1 satisfies

[TABLE]

Note that the involved constant only depends on $t_{n}$, $L$, and the operator $\mathcal{A}$.

Proof.

With $w_{n}$ and $z$ from (3.3) and (3.4), respectively, we define $U(t):=z(t)+w_{n}+\mathcal{B}^{-}g(t)$ for $t\in[t_{n},t_{n+1}]$, $n=0,\ldots,N-1$. This function satisfies

[TABLE]

Furthermore, since $\dot{U}(t)=\dot{z}(t)+\mathcal{B}^{-}\dot{g}(t)$, the function $U$ solves the PDAE

[TABLE]

on $[t_{n},t_{n+1}]$, $n=0,\ldots,N-1$ with initial value $U(t_{0})=u_{0}$. To shorten notation we define $\triangle u:=u-U$ and $\triangle\lambda:=\lambda-\Lambda$, which satisfy

[TABLE]

on each interval $[t_{n},t_{n+1}]$ with initial value $\triangle u(t_{0})=0$ if $n=0$ and $\triangle u(t_{n})=u(t_{n})-u_{n}$ otherwise. In the following, we derive estimates of $\triangle u$ on all sub-intervals. Starting with $n=0$, we have by Lemma 2.6 that

[TABLE]

By Gronwall’s lemma and $t=t_{1}=\tau$ we obtain with $c:=2c^{2}_{\mathcal{A}_{2}}{\mu}^{-1}$ the bound

[TABLE]

With the uniform Lipschitz constant $L$ we have for $n\geq 1$ that

[TABLE]

With this, we obtain similarly as in (3.5) and with Young’s inequality,

[TABLE]

Therefore, with $(1+x)\leq e^{x}$, estimate (3.5), and an iterative application of the estimate (3.6) we get

[TABLE]

for all $n=0,\ldots,N-1$. The stated estimate finally follows by the equivalence of $\|\cdot\|_{\mathcal{V}}$ and $\|\cdot\|_{\mathcal{A}}$ on $\mathcal{V}_{\ker}$, see (2.5). ∎

Remark 3.4*.*

The assumption on the step size $\tau$ only depends on the nonlinearity $f$ and not on the operator $\mathcal{A}$. Thus, this condition does not depend on the stiffness of the system and still allows large time steps.

Remark 3.5*.*

In the case of a self-adjoint operator $\mathcal{A}$, i.e., $\mathcal{A}_{2}=0$, the convergence result can also be proven by the restriction to test functions in $\mathcal{V}_{\ker}$ and the application of corresponding results for the unconstrained case, namely [HO10, Th. 2.14]. This requires similar assumptions but with $\frac{\mathrm{d}}{\mathrm{d}t}f(\cdot,u(\cdot))\in L^{\infty}(0,T;\mathcal{H})$.

We like to emphasize that this procedure is also applicable if $\mathcal{A}_{2}\neq 0$ by moving $\mathcal{A}_{2}$ into the nonlinearity $f$. This, however, slightly changes the proposed scheme, since then only $\mathcal{A}_{2}u_{n}$ enters the approximation instead of $\mathcal{A}_{2}u(t)$. In practical applications this would also require an additional effort in order to find the symmetric part of the differential operator $\mathcal{A}$ which is still elliptic on $\mathcal{V}_{\ker}$.

3.3. An alternative approach

A second approach to construct an exponential Euler scheme which is applicable to constrained systems is to formally apply the method to the corresponding singularly perturbed PDE. This approach was also considered in [HLS98] for DAEs of index 1. In the present case, we add a small term $\varepsilon\dot{\lambda}$ into the second equation of (2.4). Thus, we consider the system

[TABLE]

which can be written in operator matrix form as

[TABLE]

For this, an application of the exponential Euler method yields the scheme

[TABLE]

We introduce the auxiliary variables $(\bar{w}_{n},\bar{\nu}_{n})\in\mathcal{V}\times\mathcal{Q}$ as the unique solution to the stationary saddle point problem

[TABLE]

The included parameter $\theta\in[0,1]$ controls the consistency as outlined below. Then, the exponential Euler method can be rewritten as

[TABLE]

which allows an interpretation as the solution of a linear (homogeneous) PDE. Finally, we set $\varepsilon=0$, which leads to the following time integration scheme: Given $\bar{w}_{n}$, solve on $[t_{n},t_{n+1}]$ the linear system

[TABLE]

with initial condition $z(t_{n})=u_{n}-\bar{w}_{n}$. The approximation of $u(t_{n+1})$ is then defined through $u_{n+1}:=z(t_{n+1})+\bar{w}_{n}$.

We emphasize that the initial value of $z$ may be inconsistent. In this case, the initial value needs to be projected to $\mathcal{H}_{\ker}$, cf. Section 2.4. If the previous iterate satisfies $\mathcal{B}u_{n}=g_{n}$, then the choice $\theta=1$ yields $\mathcal{B}z(t_{n})=0$ and thus, consistency. This, however, does not imply $\mathcal{B}u_{n+1}=g_{n+1}$. On the other hand, $\theta=0$ causes an inconsistency for $z$ but guarantees $\mathcal{B}u_{n+1}=g_{n+1}$. We now turn to an exponential integrator of higher order.

4. Exponential Integrators of Second Order

This section is devoted to the construction of an exponential integrator of order two for constrained parabolic systems of the form (2.4). In particular, we aim to transfer the method given in [SWP12, Exp. 11.2.2], described by the Butcher tableau

[TABLE]

to the PDAE case. In the unconstrained case, i.e., for $\dot{v}+\mathcal{A}_{\ker}v=\tilde{f}(t,v)$ in $\mathcal{V}_{\ker}^{*}$, one step of this method is defined through

[TABLE]

Similarly as for the exponential Euler method, we will define a number of auxiliary problems in order to obtain an applicable method for parabolic systems with constraints.

4.1. The practical method

We translate the numerical scheme (4.2) to the constrained case. Let $u_{n}$ denote the given approximation of $u(t_{n})$. Then, the first step is to perform one step of the exponential Euler method, cf. Algorithm 1, leading to $u_{n+1}^{\text{Eul}}$. Second, we compute $w_{n}^{\prime}$ as the solution of the stationary problem

[TABLE]

and $w_{n}^{\prime\prime}$ as the solution of

[TABLE]

Note that, due to the recursion formula (2.1), $w_{n}^{\prime}$ and $w_{n}^{\prime\prime}$ satisfy the identity

[TABLE]

It remains to compute $\varphi_{0}(-\tau\mathcal{A}_{\ker})w_{n}^{\prime\prime}$ and thus, to solve a linear dynamical system with starting value $w_{n}^{\prime\prime}$. More precisely, we consider the homogeneous system (3.4) on the time interval $[t_{n},t_{n+1}]$ with initial value $z(t_{n})=w_{n}^{\prime\prime}$. The solution at time $t_{n+1}$ then defines the new approximation by

[TABLE]

Note that the consistency is already guaranteed by the exponential Euler step which yields  $\mathcal{B}u_{n+1}=\mathcal{B}u^{\text{Eul}}_{n+1}=g_{n+1}$. The resulting exponential integrator is summarized in Algorithm 2.

4.2. Convergence analysis

In this subsection we aim to prove the second-order convergence of Algorithm 2 when applied to parabolic PDAEs of the form (2.4). For this, we examine two cases. First, we consider a nonlinearity with values in $\mathcal{V}$, i.e., we assume $f\colon[0,T]\times\mathcal{V}\to\mathcal{V}$. Further, we assume $\mathcal{A}$ to be self-adjoint, meaning that  $\mathcal{A}_{2}=0$. Note that this may be extended to general $\mathcal{A}$ as mentioned in Remark 3.5. In this case, the convergence analysis is based on the corresponding results for unconstrained systems. Second, we consider the more general case with nonlinearities $f\colon[0,T]\times\mathcal{V}\to\mathcal{H}$. Here, it can be observed that the convergence order drops to $\nicefrac{{3}}{{2}}$. Note, however, that this already happens in the pure PDE case.

Theorem 4.1** (Second-order scheme).**

In the setting of Section 2.2, including Assumptions 2.2 and 2.3, we assume that $\mathcal{A}$ is self-adjoint and that $f\colon[0,T]\times\mathcal{V}\to\mathcal{V}$ is two times Fréchet differentiable in a strip along the exact solution $u$ with uniformly bounded derivatives. Further we assume that the right-hand side $g$ and $u$ are sufficiently smooth, the latter with derivatives in $\mathcal{V}$. Then, the approximation obtained by Algorithm 2 is second-order accurate, i.e.,

[TABLE]

Proof.

We reduce the procedure performed in Algorithm 2 to the unconstrained case. For this, assume that $u_{n}=u_{\text{ker},n}+\mathcal{B}^{-}g_{n}\in\mathcal{V}$ is given with $u_{\text{ker},n}\in\mathcal{V}_{\ker}$ and that $u^{\text{Eul}}_{n+1}$ denotes the outcome of a single Euler step, cf. Algorithm 1. By $u_{\text{ker},n+1}^{\text{Eul}}$ we denote the outcome of a Euler step for the unconstrained system

[TABLE]

with $\tilde{f}$ defined by $\tilde{f}(t,u_{\text{ker}}):=\iota_{0}\,[f(t,u_{\text{ker}}+\mathcal{B}^{-}g(t))-\mathcal{B}^{-}\dot{g}(t)]$ and initial data $u_{\text{ker},n}$. For this, we know that $u^{\text{Eul}}_{n+1}=u_{\text{ker},n+1}^{\text{Eul}}+\mathcal{B}^{-}g_{n+1}$. By the given assumptions, it follows from [HO10, Th. 2.17] that

[TABLE]

defines a second-order approximation of $u_{\text{ker}}(t_{n+1})$. This in turn implies that

[TABLE]

satisfies the error estimate

[TABLE]

It remains to show that $u_{n+1}$ is indeed the outcome of Algorithm 2. Following the construction in Section 4.1, we conclude that

[TABLE]

with $w_{n}^{\prime}$ and $w_{n}^{\prime\prime}$ denoting the solutions of (4.3) and (4.4), respectively. Finally, $\varphi_{0}(-\tau\mathcal{A}_{\ker})w_{n}^{\prime\prime}$ is computed in line 7 of Algorithm 2 which completes the proof. ∎

Up to now we have assumed that $f$ maps to $\mathcal{V}$, leading to the desired second-order convergence. In the following, we reconsider the more general case with $f\colon[0,T]\times\mathcal{V}\to\mathcal{H}$. For PDEs it is well-known that the exponential integrator given by the Butcher tableau (4.1) has a reduced convergence order if we assume $\tfrac{\mathrm{d}^{2}}{\mathrm{d}t^{2}}f(\cdot,u(\cdot))\in L^{\infty}(0,T;\mathcal{H})$, cf. [HO05a, Th. 4.3]. This carries over to the PDAE case.

Theorem 4.2** (Convergence under weaker assumptions on $f$).**

Consider the assumptions from Theorem 2.7 and let the step size $\tau$ be sufficiently small such that the discrete solution $u_{n}$ lies in a strip along $u$, where $f$ is locally Lipschitz continuous with a uniform constant $L>0$. Further assume that $g\in H^{3}(0,T;\mathcal{Q}^{*})$. If the exact solution of (2.4) satisfies $f(\,\cdot\,,u(\,\cdot\,))\in H^{2}(0,T;\mathcal{H})$, then the approximation $u_{n}$ obtained by Algorithm 2 satisfies the error bound

[TABLE]

Note that the involved constant only depends on $t_{n}$, $L$, and the operator $\mathcal{A}$.

Proof.

Let $U^{\text{Eul}}$ be the function constructed in the proof of Theorem 3.3 which satisfies $U^{\text{Eul}}(t_{n})=u_{n}$ and $U^{\text{Eul}}(t_{n+1})=u^{\text{Eul}}_{n+1}$ and set $U(t):=U^{\text{Eul}}(t)+z(t)-w_{n}^{\prime\prime}+\frac{t-t_{n}}{\tau}w_{n}^{\prime}$. This function satisfies

[TABLE]

Note that the estimates (3.5) and (3.6) are still valid if one replaces $u_{n+1}$ by $U^{\text{Eul}}(t_{n+1})$ on the left-hand side of these estimates. As in the proof of Theorem 3.3, we can interpret $U$ as the solution of a PDAE on $[t_{n},t_{n+1}]$. The corresponding right-hand sides are then given by

[TABLE]

for the dynamic equation and $g(t)$ for the constraint. Then, by Young’s inequality, Gronwall’s lemma, and error bounds for the Taylor expansion we get

[TABLE]

with $c=2c^{2}_{\mathcal{A}_{2}}{\mu}^{-1}$. The stated error bound then follows by an iterative application of the previous estimate together with the estimates (3.5), (3.6) and the norm equivalence of $\|\cdot\|_{\mathcal{V}}$ and $\|\cdot\|_{\mathcal{A}_{1}}$. ∎

The actual performance of the proposed scheme is presented in the numerical experiments of Section 5. We close this section with remarks on alternative second-order schemes.

4.3. A class of second-order schemes

The analyzed scheme (4.1) is a special case of a one-parameter family of exponential Runge-Kutta methods described by the tableau

[TABLE]

with positive parameter $c_{2}>0$, cf. [HO10]. Therein, $\varphi_{1}$ stands for $\varphi_{1}(-\tau\mathcal{A}_{\ker})$, whereas $\varphi_{1,2}$ is defined by $\varphi_{1}(-c_{2}\tau\mathcal{A}_{\ker})$. Obviously, we regain (4.1) for $c_{2}=1$.

For $c_{2}\neq 1$, the resulting scheme for constrained systems calls for two additional saddle point problems in order to compute $\mathcal{B}^{-}g(t_{n}+c_{2}\tau)$ and $\mathcal{B}^{-}\dot{g}(t_{n}+c_{2}\tau)$. This then leads to an exponential integrator summarized in Algorithm 3 with the abbreviations

[TABLE]

We emphasize that all convergence results of Theorems 4.1 and 4.2 transfer to this family of second-order integrators. In a similar manner, Runge-Kutta schemes of higher order may be translated to the here considered constrained case.

5. Numerical Examples

In this final section we illustrate the performance of the introduced time integration schemes for two numerical examples. The first example is a heat equation with nonlinear dynamic boundary conditions. In the second experiment, we consider the case of a non-symmetric differential operator for which the theory is not applicable.

Since exponential integrators for PDAEs are based on the exact solution of homogeneous systems of the form (3.4), we first discuss the efficient solution of such systems.

5.1. Efficient solution of homogeneous DAEs with saddle-point structure

This subsection is devoted to the approximation of $z(t)$, which is needed in line 5 of Algorithm 1 and in line 7 of Algorithm 2. Given a spatial discretization, e.g., by a finite element method, the PDAE (3.4) turns into a DAE of index $2$, namely

[TABLE]

with consistent initial value $x(0)=x_{0}$, $Bx_{0}=0$. The matrices satisfy $M,A\in\mathbb{R}^{n\times n}$ and $B\in\mathbb{R}^{m\times n}$ with $m\leq n$. Here, the mass matrix $M$ is symmetric, positive definite and $B$ has full rank. The goal is to find an efficient method to calculate the solution $x$ at a specific time point $t\in[0,T]$.

Let us first recall the corresponding ODE case. There exist various methods to approximate the solution $x(t)=e^{-At}x_{0}$ of the linear ODE $\dot{x}(t)+Ax(t)=0$ with initial condition $x(0)=x_{0}$, $A\in\mathbb{R}^{n\times n}$, for an overview see [MVL78]. This includes Krylov subspace methods to approximate the action of the matrix exponential $e^{-At}$ to a vector, see [Saa92, HL97, EE06], but also methods based on an interpolation of $e^{-At}x_{0}$ by Newton polynomials [CO09]. The first approach is based on the fact that the solution $e^{-At}x_{0}=\sum_{k=0}^{\infty}\frac{1}{k!}(-At)^{k}x_{0}$ is an element of the Krylov subspace

[TABLE]

Now, we approximate $e^{-At}x_{0}$ by an element of $\mathcal{K}_{r}$ with $r$ relatively small compared to $n$. For this, we generate an orthogonal basis of $\mathcal{K}_{r}$ using the Arnoldi algorithm with $v_{1}={x_{0}}/{\|x_{0}\|}$ as initial vector. This yields $-V_{r}^{T}AV_{r}=H_{r}$ with an isometric matrix $V_{r}\in\mathbb{R}^{n\times r}$ and an upper Hessenberg matrix $H_{r}\in\mathbb{R}^{r\times r}$. Since the columns of $V_{r}$ are orthonormal and span $\mathcal{K}_{r}$, $H_{r}$ is the orthogonal projection of $-A$ onto $\mathcal{K}_{r}$. Therefore, it is reasonable to use the approximation

[TABLE]

with unit basis vector $e_{1}\in\mathbb{R}^{r}$, cf. [HL97]. We like to emphasize that the Arnoldi algorithm does not use the explicit representation of $A$ but only its action onto a vector.

We return to the DAE case (5.1). By [EM13, Th. 2.2] there exists a matrix $X\in\mathbb{R}^{n\times n}$ such that the solution $x$ of (5.1) with arbitrary consistent initial value $x_{0}$ is given by $x(t)=e^{Xt}x_{0}$. Furthermore, there exists a matrix-valued function $\Lambda\in C^{\infty}([0,\infty);\mathbb{R}^{m\times n})$ with $\lambda(t)=\Lambda(t)x_{0}$. To calculate the action of $X$ we note that by (5.1b) also $B\dot{x}=0$ holds. We define $y:=Xx_{0}$ and $\mu:=\Lambda(0)x_{0}$. Then with equation (5.1a), $B\dot{x}=0$, and $t\to 0^{+}$ we get

[TABLE]

Since the solution of (5.2) is unique, its solution $y$ describes the action of $X$ applied to $x_{0}$. As a result, we can approximate the solution of the DAE (5.1) in an efficient manner by using $x(t)=e^{Xt}x_{0}$, the saddle point problem (5.2), and Krylov subspace methods. For the numerical experiments we have adapted the code provided in [NW12].

Remark 5.1*.*

Given an approximation $x_{t}\approx x(t)$, the solution $\mu$ of (5.2) with right-hand side $-Ax_{t}$ provides an approximation of the Lagrange multiplier $\lambda(t)$.

Remark 5.2*.*

Since the saddle point problem (5.2) has to be solved several times in every time step, the numerical solution $\tilde{x}$ of (5.1) may not satisfy the constraint (5.1b) due to round-off errors. To prevent a drift-off, one can project $\tilde{x}$ onto the kernel of $B$ – by solving an additional saddle point problem.

5.2. Nonlinear dynamic boundary conditions

In this first example we revisit Example 2.5 and consider the linear heat equation with nonlinear dynamic boundary conditions, cf. [SW10]. More precisely, we consider the system

[TABLE]

with $\alpha=1$, $\kappa=0.02$, and the nonlinearity $f_{\Gamma}(t,u)(x)=3\cos(2\pi t)-\sin(2\pi x)-u^{3}(x)$. As initial condition we set $u(0)=u_{0}=\sin(\pi x)\cos(5\pi y/2)$. Following [Alt19], we can write this in form of a PDAE, namely as

[TABLE]

with spaces $\mathcal{V}=H^{1}_{\Gamma_{D}}(\Omega)\times H^{1/2}(\Gamma_{\text{dyn}})$, $\mathcal{H}=L^{2}(\Omega)\times L^{2}(\Gamma_{\text{dyn}})$, $\mathcal{Q}=[H^{1/2}(\Gamma_{\text{dyn}})]^{*}$ and constraint operator $\mathcal{B}(u,p)=u|_{\Gamma_{\text{dyn}}}-p$. Here, $p$ denotes a dummy variable modeling the dynamics on the boundary $\Gamma_{\text{dyn}}$. The constraint (5.4b) couples the two variables $u$ and $p$. This example fits into the framework of this paper with $g\equiv 0$. Further, the nonlinearity satisfies the assumptions of the convergence results in Theorems 3.3 and 4.2 due to well-known Sobolev embeddings, see [Rou05, p. 17f].

For the simulation we consider a spatial discretization by bilinear finite elements on a uniform mesh with mesh size $h=1/128$. The initial value of $p$ is chosen in a consistent manner, i.e., by $u_{0}|_{\Gamma_{\text{dyn}}}$. An illustration of the dynamics is given in Figure 5.1. The convergence results of the exponential Euler scheme of Section 3 and the exponential integrator introduced in Section 4 are displayed in Figure 5.2 and show first and second-order convergence, respectively. Note that the theory only ensures an order $\nicefrac{{3}}{{2}}$ for the second integrator. However, the spatial discretization acts as a regularization in space, i.e., we expect order $\nicefrac{{3}}{{2}}$ only for $h\to 0$ and non-smooth initial data.

Finally, we note that the computations remain stable for very coarse step sizes $\tau$, since we do not rely on a CFL condition here.

5.3. A non-symmetric example

In this final example we consider a case for which Assumption 2.3 is not satisfied. More precisely, we consider the coupled system

[TABLE]

with initial value

[TABLE]

and the constraint $u(t,1)-v(t,1)=g(t)=e^{2t}-1$. At the other boundary point $x=0$ we prescribe homogeneous Dirichlet boundary conditions. In this example, the operator $\mathcal{A}$ has the form $-[\partial_{xx},\,\partial_{xx};\,-\operatorname{id},\,\partial_{xx}]$. Thus, the non-symmetric part $\mathcal{A}_{2}$ includes a second-order differential operator which contradicts Assumption 2.3. As a consequence, non of the convergence results in this paper apply.

The numerical results are shown in Figure 5.3, using a finite element discretization with varying mesh sizes $h$. One can observe that the exponential Euler scheme still converges with order $1$, whereas the second-order scheme introduced in Section 4 clearly converges with a reduced rate. Moreover, the rate decreases as the mesh size $h$ gets smaller. By linear regression one can approximate the convergence rate as a value between $1.40$ (coarsest mesh, $h=1/32$) and $1.34$ (finest mesh, $h=1/256$). Thus, the convergence rate is strictly below $\nicefrac{{3}}{{2}}$.

A deeper analysis with fractional powers of $\mathcal{A}$ may predict the exact convergence rate, cf. [HO05a, Th. 4.2 & Th. 4.3]. However, this is a task for future work.

6. Conclusion

In this paper, we have introduced a novel class of time integration schemes for semi-linear parabolic equations restricted by a linear constraint. For this, we have combined exponential integrators for the dynamical part of the system with (stationary) saddle point problems for the ’algebraic part’ of the solution. As a result, we obtain exponential integrators for constrained systems of parabolic type for which we have proven convergence of first and second order, respectively. Numerical experiments illustrate the obtained results.