Strang splitting method for semilinear parabolic problems with inhomogeneous boundary conditions: a correction based on the flow of the nonlinearity

TL;DR

This paper introduces a new correction method for the Strang splitting technique applied to semilinear parabolic problems with inhomogeneous boundary conditions, reducing computational effort and improving performance especially for stiff nonlinearities.

Contribution

A novel correction based on the flow of the nonlinearity that requires only one evaluation per step, simplifying implementation and enhancing efficiency for stiff problems.

Findings

The new correction reduces computational effort compared to previous methods.

Numerical experiments show improved performance with stiff nonlinearities.

The method achieves convergence for smooth nonlinearities.

Abstract

The Strang splitting method, formally of order two, can suffer from order reduction when applied to semilinear parabolic problems with inhomogeneous boundary conditions. The recent work [L .Einkemmer and A. Ostermann. Overcoming order reduction in diffusion-reaction splitting. Part 1. Dirichlet boundary conditions. SIAM J. Sci. Comput., 37, 2015. Part 2: Oblique boundary conditions, SIAM J. Sci. Comput., 38, 2016] introduces a modification of the method to avoid the reduction of order based on the nonlinearity. In this paper we introduce a new correction constructed directly from the flow of the nonlinearity and which requires no evaluation of the source term or its derivatives. The goal is twofold. One, this new modification requires only one evaluation of the diffusion flow and one evaluation of the source term flow at each step of the algorithm and it reduces the computational effort to construct the correction. Second, numerical experiments suggest it is well suited in the case where the nonlinearity is stiff. We provide a convergence analysis of the method for a smooth nonlinearity and perform numerical experiments to illustrate the performances of the new approach.

Tables

Table 1. Table 1 : Number of evaluations of the diffusion flows ( 1.5 ,  1.8 ) and of the source term flows ( 1.4 ,  1.6 ) needed by the considered splitting methods.

Splitting methods	Diffusion flows	Source term flows	Total number of
	for one step	for one step	flows for $n$ steps
Classical Strang (1.2)	1	1	$2n+1$
StrangM3 (1.7)	2	1	$3n$
StrangM5a (3.2, 3.3)	1	2	$3n$
StrangM5b (3.2, 3.4)	1	1	$2n+1$

Equations

$${\partial }_{t}u=Du+f(u)\mathrm{in}\Omega ,Bu=b\mathrm{on}\partial \Omega ,u(0)=u_0,$$

$$u_{n+1}={\varphi }_{\frac{\tau }{2}}^D\circ {\varphi }_{\tau }^f\circ {\varphi }_{\frac{\tau }{2}}^D(u_n)$$

$$u_{n+1}={\varphi }_{\frac{\tau }{2}}^f\circ {\varphi }_{\tau }^D\circ {\varphi }_{\frac{\tau }{2}}^f(u_n),$$

$${\partial }_{t}u=f(u),u(0)=u_0,$$

$${\partial }_{t}u=Du\mathrm{in}\Omega ,Bu=b\mathrm{on}\partial \Omega ,u(0)=u_0.$$

$${\partial }_{t}u=f(u)-q_n.$$

$$u_{n+1}={\varphi }_{\frac{\tau }{2}}^{D+q_n}\circ {\varphi }_{\tau }^{f-q_n}\circ {\varphi }_{\frac{\tau }{2}}^{D+q_n}(u_n),$$

$${\partial }_{t}u=Du+q_n\mathrm{in}\Omega ,Bu=b\mathrm{on}\partial \Omega .$$

$$B{q}_n=Bf(u_n)+\mathcal{O}(\tau )\mathrm{on}\partial \Omega .$$

$${\mathcal{S}}_{\tau }(u_n)=({\varphi }_{\frac{\tau }{2}}^f\circ {\varphi }_{\frac{\tau }{2}}^{-q_n}\circ {\varphi }_{\tau }^{D+q_n}\circ {\varphi }_{\frac{\tau }{2}}^{-q_n}\circ {\varphi }_{\frac{\tau }{2}}^f)(u_n)$$

$${\partial }_{t}u=Du+f(u)\mathrm{in}\Omega ,Bu=b\mathrm{on}\partial \Omega ,u(0)=u_0.$$

$$D=\sum _{i,j=1}^n{a}_{ij}(x){\partial }_{ij}+\sum _{i=1}^n{a}_i(x){\partial }_i+a(x)I,$$

$$\forall x\in \overline{\Omega },\forall \xi \in {\mathbb{R}}^n,\sum _{i,j=1}^n{a}_{ij}(x){\xi }_i{\xi }_j\ge \lambda |\xi {|}^2,$$

$$B=\sum _{i=1}^n{\beta }_i(x){\partial }_i+\alpha (x)I,$$

$$\underset{x\in \partial \Omega }{\inf }\left|\sum _{i=1}^n{\beta }_i(x){\nu }_i(x)\right|>0,$$

$${\partial }_t\tilde{u}=D\tilde{u}+f(\tilde{u}+z)+Dz-{\partial }_{t}z\mathrm{in}\Omega ,B\tilde{u}=0\mathrm{on}\partial \Omega ,\tilde{u}(0)=u_0-z(0).$$

$$Av=Dv\forall v\in D(A)=\{u\in W^{2,p}(\Omega ):Bu=0\mathrm{in}\partial \Omega \}.$$

$$\Vert (-A{)}^{\alpha }{\mathrm{e}}^{tA}\Vert \le \frac{C}{t^{\alpha }}\alpha \ge 0t>0.$$

$${\phi }_1(\tau A)={\int }_0^{\tau }{\mathrm{e}}^{(\tau -s)A}\frac{1}{\tau }\mathrm{ds}.$$

$$W^{1,p}(\Omega )\subset D((-A{)}^{\alpha }).$$

$$u(t_{n+1})$$

$$+{\int }_0^{\tau }{\mathrm{e}}^{(\tau -s)A}(f(u(t_n+s))+D{z}_n(s)-{\partial }_t{z}_n(s))\mathrm{ds}.$$

$${\partial }_{t}u=-q_{n}u(0)=w_n,$$

$$u_{n+1}={\mathcal{S}}_n(u_n)={\varphi }_{\frac{\tau }{2}}^f\circ {\varphi }_{\frac{\tau }{2}}^{-q_n}\circ {\varphi }_{\tau }^{D+q_n}\circ {\varphi }_{\frac{\tau }{2}}^{-q_n}\circ {\varphi }_{\frac{\tau }{2}}^f(u_n),$$

$$B{q}_n$$

$$u_n={\varphi }_{\frac{\tau }{2}}^f\circ {\varphi }_{\tau }^D\circ ({\varphi }_{\tau }^f\circ {\varphi }_{\tau }^D{)}^{n-1}\circ {\varphi }_{\frac{\tau }{2}}^f(u_0),$$

$$\Vert u_n-u(t_n)\Vert \le C{\tau }^{2}0\le n\tau \le T,$$

$$\Vert u_n-u(t_n)\Vert \le C\tau 0\le n\tau \le T,$$

$$f(u(t_n+s))-q_n={\varphi }_n+\mathcal{O}(\tau ),$$

$$\Vert u_n-u(t_n)\Vert \le C{\tau }^2(1+|\log \tau |)0\le n\tau \le T,$$

Full text

Key words. Strang splitting, diffusion-reaction equation, nonhomogeneous boundary conditions, order reduction, stiff nonlinearity

AMS subject classifications. 65M12, 65L04

1 Introduction

In this paper, we consider a parabolic differential equation of the form

[TABLE]

where $D$ is a linear diffusion operator and $f$ is a nonlinearity. A natural method for approximating (1.1) are splitting methods. The idea is to divide the main equation (1.1) into two auxiliary subproblems (1.4) and (1.5) so one can use specific numerical methods for both subproblems to enhance the global efficiency of the computation of (1.1). Let $N\in\mathbb{N}$ and let $\tau=\frac{T}{N}$ be the time step. Then, one step of the classical Strang splitting is either

[TABLE]

or alternatively

[TABLE]

where $\phi^{f}_{t}(u_{0})$ is the flow after time $t$ of

[TABLE]

and $\phi_{t}^{D}(u_{0})$ is the flow after time $t$ of

[TABLE]

The Strang splitting, when applied to ODE with a sufficiently smooth solution, is a method of order two. However, when the Strang splitting is applied to solve the problem (1.1), a reduction of order can be observed in the case of non homogeneous boundary conditions as shown in [2, 3]. The reason is that $Bu$ is not left invariant through the flow $\phi_{t}^{f}$ and therefore leaves the domain of $D$, which creates a discontinuity at $t=0$ in the flow $\phi_{t}^{D}$. In this case the Strang splitting has in general a fractional order of convergence between one and two [3, section 4.3]. In [2, 3], a modification of the Strang splitting is given to recover the order two. The main idea in [3] is to find a function $q_{n}$ such that $Bu$ is now left invariant by $\phi_{t}^{f-q_{n}}$, the exact flow of

[TABLE]

One step of the modified splitting in [3] is then

[TABLE]

where $\phi_{t}^{D+q_{n}}$ is the exact flow of

[TABLE]

Numerically, one can choose any smooth function $q_{n}$ such that

[TABLE]

Several options to construct $q_{n}$ are presented in [1]. One challenge is then to find a correction $q_{n}$ that is both cheap to compute and minimizes the constant of error.

In this paper, we give a new modification of the classical Strang splitting that removes the order reduction and requires only one evaluation of the diffusion flow and one evaluation of the source term flow at each step of the algorithm and allows a cheaper and easier to implement construction of $q_{n}$. As illustrated in the experiments, this new construction performs better for the case of a stiff reaction. The idea is to leave $Bu$ unpreserved at the boundary through the flow $\phi^{f}_{\tau}$ and then apply a correction $q_{n}$ afterwards that brings back the solution to the domain of $D$. This new splitting is then

[TABLE]

and the correction $q_{n}$ is constructed such that $Bq_{n}=\frac{2}{\tau}(B\phi^{f}_{\tau/2}(u_{n})-Bu_{n})$ on $\partial\Omega$. The correction $q_{n}$ is now constructed from the output of the flow $\phi^{f}_{\tau/2}(u_{n})$ and not directly from the nonlinearity $f$. Note that the computation of $\phi^{-q_{n}}_{\frac{\tau}{2}}$ and $q_{n}$ requires no evaluation of $f$, which is particularly useful when $f$ is costly. More importantly, in many situations, the flow $\phi^{f}_{\tau/2}(u_{n})$ is smoother than the nonlinearity $f$ itself, which can avoid the possible instability due to the eventual stiffness of the reaction.

In section 2, we give the appropriate framework for the convergence analysis of this modified splitting. In section 3, we describe the new modification we consider in this paper. In section 4, we prove that the method is of global order two under the hypotheses made in section 2. In section 5, we present some numerical experiments to illustrate the performance of the new approach.

2 Analytical framework

In this section, we describe the appropriate analytical framework that we consider in this paper. We choose the framework described in [8, Chapter 3]. The notation is similar to that used in [8]. Let $\Omega\subset\mathbb{R}^{d}$ be a bounded connected open set with a $C^{2}$ boundary $\partial\Omega$. We consider the following semilinear parabolic problem on $\Omega\times[0,T]$, $T\geq 0$.

[TABLE]

The differential operator $D$ is defined by

[TABLE]

where the matrix $(a_{ij}(x))\in\mathbb{R}^{d\times d}$ is assumed symmetric and there exists $\lambda>0$ such that

[TABLE]

and $a_{ij},\ a_{i},\ a$ are assumed continuous, $a_{i,j},\ a_{i},\ a\in C(\Omega,\mathbb{R})$ . Let $B$ be the linear operator

[TABLE]

where we assume the uniform non tangentiality condition

[TABLE]

where $\nu(x)$ is the exterior normal unit vector at $x\in\partial\Omega$. We assume that the functions $\beta_{i}$ and $\alpha$ are continuously differentiable, $\beta_{i},\alpha\in C^{1}(\partial\Omega,\mathbb{R})$ and $b$ is continuously differentiable, $b\in C^{1}([0,T],W^{2,p}(\partial\Omega))$. We follow next the construction made in [3] to take benefit of homogeneous boundary conditions. Let $z\in C^{1}([0,T],W^{2,p}(\Omega))$ satisfying $Bz=b$ on $\partial\Omega$. We define $\tilde{u}=u-z$ and (1.1) becomes

[TABLE]

We define the linear operator $A$ as

[TABLE]

Under those conditions $-A$ is a sectorial operator and therefore $A$ is the generator of an analytic semigroups $e^{tA}$ (see [8], Chapter 3). In particular, the operator $A$ satisfies the parabolic smoothing property that we use intensively throughout this paper,

[TABLE]

We denote

[TABLE]

We observe that $\tau A\varphi_{1}(\tau A)=\mathcal{O}(1)$ is a bounded operator. We recall the following theorem, a direct consequence of [5, Theorem 8.1′], valid for $1<p<\infty$ (see also [4, Theorem 1 and 2] in English for the case $p=2$), which states that there exists $\alpha>0$ such that $D((-A)^{\alpha})$ becomes free of the boundary conditions.

Theorem 2.1**.**

Let the differential operator $A$ be define as in section 2. Then, there exists $\alpha>0$ such that

[TABLE]

We ask $f$ to satisfy the following. Let $U\subset W^{2,p}(\Omega)$ be a neighborhood of the exact solution $u$. Then we require the nonlinearity $f$ to be twice continuously differentiable in $U$ with values in $W^{2,p}(\Omega)$, $f\in C^{2}(U,W^{2,p}(\Omega))$. We refer to the discussion in [3, Section 4] for possibly relaxing the hypotheses made on $f$. We assume that the solution $u$ of (1.1) is twice continuously differentiable, $u\in C^{2}([0,T],W^{2,p}(\Omega))$. The exact solution of (1.1) can be expressed using the variation of constant formula,

[TABLE]

3 Description of the method

In this paper, we describe a new modification for the Strang splitting that we call the five parts modified Strang splitting. The idea of this new modification is to compose the flow of the nonlinearity $w_{n}=\phi^{f}_{\frac{\tau}{2}}(u_{n})$ with a projection $\phi^{-q_{n}}_{\frac{\tau}{2}}(w_{n})=w_{n}-\frac{\tau}{2}q_{n}$, the exact flow of

[TABLE]

where $q_{n}$ is independent of time, in the spirit of projection methods used in the context of geometric numerical integration (see [6, Chapter IV.4]). The splitting algorithm that we propose and analyze in this paper is given by

[TABLE]

where $\phi^{D+q_{n}}_{\tau}$ is the flow of (1.8) and $(q_{n})_{n\in\{0,\ldots,N\}}$ is a sequence of correctors satisfying one of the two following conditions on the boundary $\partial\Omega$,

[TABLE]

or alternatively

[TABLE]

(see Remark 3.2 below). In the interior of $\Omega$, we require $q_{n}$ to be in $W^{2,p}(\Omega)$. A possibility, to construct $q_{n}$ in $\Omega$, is to choose $q_{n}$ to be harmonic if this is possible or to use a smoothing iterative algorithm. For more details on how to construct the correction $q_{n}$ on the interior of the domain, see [1]. We also assume $(q_{n})_{n\in\{0,\ldots,N\}}$ uniformly bounded, that is there exists a constant $C$ independent of $n,\tau$ and $N$, such that $\|q_{n}\|_{L^{p}(\Omega)}\leq C$. We observe that $\frac{2}{\tau}(B\phi^{f}_{\tau/2}(u_{n})-Bu_{n})$ is a finite difference approximation of $\partial_{t}\phi^{f}_{\tau/2}(u_{n})=f(u_{n})$, and hence this new condition is close to (1.9).

Remark 3.1**.**

In contrast to the correction of [3], the correction $q_{n}$ for the five parts modified Strang splitting (3.2) is constructed directly from the flow of the nonlinearity $\phi^{f}_{t}$ and not from the nonlinearity $f$ itself. The modified splitting of [3] has a good behavior when the nonlinearity $f$ is not stiff and cheap to compute as analyzed and illustrated numerically in [3]. However, in the case of a stiff nonlinearity $f$, the modification for the splitting in [3] can lead to instability in contrast to (3.2) as shown in the experiments (see section 5). Furthermore, if the nonlinearity $f$ is very costly to compute, the correction in [3] requires an additional cost that can be substantial. In comparison, the construction of the correction $q_{n}$ for (3.2) requires no evaluation of $f$ or its derivatives. We also observe that, in the extreme case where the diffusion D is zero, the five parts modified Strang splitting (3.2) becomes exact analogously to the classical Strang splitting methods (1.2) and (1.3). This later property does not hold for the modified splitting in [3]. (Note that the flows $\phi^{q}_{\frac{\tau}{2}}$ and $\phi^{f-q}_{\tau}$ do not commute in general).

Remark 3.2**.**

When implementing the classical Strang splitting, it is often computationally advantageous to compose the flows $\phi^{f}_{\frac{\tau}{2}}$ that appear in the splitting, that is $\phi^{f}_{\frac{\tau}{2}}\circ\phi^{f}_{\frac{\tau}{2}}=\phi^{f}_{\tau}$. The numerical approximation $u_{n}$ of $u(t_{n})$ is then

[TABLE]

which makes the classical Strang splitting have the same cost as the Lie Trotter splitting with only one evaluation of $\phi^{f}_{\tau}$ per time step. If we use the correction (3.3), we need then to compute $Bu_{k}$ and this idea does not apply since the algorithm requires $Bu_{k}$. However, if we use the correction (3.4) instead, we can implement the five parts modified Strang splitting (3.2) as explained above for the classical Strang splitting. Note that this is an advantageous implementation that cannot be used with the method presented in [3].

4 Convergence analysis

We prove in this section that, using the framework and assumptions described in section 2, the five parts modified Strang splitting method (3.2) is of global order of convergence two and thus avoids order reduction phenomena.

Theorem 4.1**.**

Under the assumption of section 2 the five parts modified Strang splitting (3.2) satisfies the bound

[TABLE]

for all $\tau$ small enough, and where the constant C depends on T but is independent on $\tau$ and n.

We start by showing the following proposition, which states that the five parts splitting is at least first order convergent.

Proposition 4.2**.**

Under the assumptions of section 2 the five parts modified Strang splitting (3.2) satisfies

[TABLE]

for all $\tau$ small enough, and where the constant C depends on T but is independent on $\tau$ and n.

We need this result to justify the condition (3.3) and (3.4) for the construction of $q_{n}$, that is we need to show that $q_{n}$ satisfies

[TABLE]

with $\phi_{n}\in D(A)$ and $A\phi_{n}=\mathcal{O}(1)$.

The proof of Theorem 4.1 relies on Theorem 2.1 (from [5]). The proof of Theorem 2.1 uses sophisticated tools from interpolation theory. Since all the arguments in our proofs do not require any knowledges of interpolation theory, we decide to present first the proof without using Theorem 2.1 and obtain Proposition 4.3 below. We then explain how we use Theorem 2.1.

Proposition 4.3**.**

Under the assumptions of section 2 the five parts modified Strang splitting (3.2) satisfies

[TABLE]

for all $\tau$ small enough, and where the constant C depends on T but is independent on $\tau$ and n.

For all the convergence analysis, for a function $\phi$ in $L^{p}(\Omega)$, we shall use the following notation for $k=0,1,\ldots$:

[TABLE]

where $C$ is independent of $\tau$ and where $\tau$ is assumed small enough.

4.1 Quadrature error analysis

The main idea of the convergence analysis is to approximate the integrals of the form $\int_{0}^{\tau}\mathrm{e}^{(\tau-s)A}\widehat{\psi}(s)ds$ with quadrature formulas, using the Peano integral representation of the error. This idea is not new in the literature and is used, for example, in [7, 3, 2]. Indeed one cannot apply quadrature formulas naively to such integrals since, for a smooth function $s\mapsto\hat{\psi}(s)$, the integrand $\mathrm{e}^{(\tau-s)A}\widehat{\psi}(s)$ is not smooth enough in general with respect to $s$. We recall that to state that the (local) error of a $k$th order quadrature formula is $\mathcal{O}(\tau^{k+1})$, the integrand needs to be $k$ times continuously differentiable with respect to $s$ but this assumption does not hold in the proofs in general. Hence, we shall prove refined estimates for the order $k=2,3$ (see Lemmas 4.4, 4.5, and 4.6). We show in Lemma 4.4 below, with the help of the parabolic smoothing property, that if $\widehat{\psi}$ is close to the domain of $A$, that is if

[TABLE]

is satisfied, then first and second order quadrature formulas regain partially their accuracy. In [3], the authors prove this statement is true for the left rectangle quadrature formula and the midpoint rule. Since we need such results for various quadrature formulas, we prove instead it is true for a general quadrature formula since it adds no difficulties to the proof. The first of the two lemmas that follow deals with quadrature formulas of order one. The second lemma deals with quadrature formulas of order two.

Lemma 4.4**.**

Let $s\mapsto\widehat{\psi}(s)$ be a continuously differentiable function with values in $L^{p}(\Omega)$, $1<p<\infty$, and let $\psi(s)=\mathrm{e}^{(\tau-s)A}\widehat{\psi}(s)$. Let $Q(\psi)=\tau\sum_{k=1}^{m}b_{k}\psi(\tau c_{k})$ be a quadrature formula that approaches the integral $\int_{0}^{\tau}\psi(s)ds$ such that $\sum_{k=1}^{m}b_{k}=1$.

Then the quadrature error $E$ satisfies

[TABLE]

If additionally $\widehat{\psi}(s)=\phi_{0}+\mathcal{O}(\tau)$ with $\phi_{0}\in D(A)$ and $A\phi_{0}=\mathcal{O}(1)$, then

[TABLE]

Proof.

Since $\psi$ is uniformly bounded on $[0,\tau]$, the first result follows. Let us assume that the condition (4.2) is satisfied.

Let us compute the first derivative of $\psi$,

[TABLE]

Let $Q_{l}(\psi)=\tau\psi(0)$ be the left rectangle quadrature formula. We prove that every first order quadrature formula $Q$ satisfies

[TABLE]

First, we observe that $\left\|\tau A\varphi_{1}(\tau A)\widehat{\psi}(0)\right\|=\mathcal{O}(\tau)$. Indeed

[TABLE]

We extend $\psi(c\tau)$ around [math],

[TABLE]

which proves (4.3) since

[TABLE]

Therefore, we only need to show that

[TABLE]

We write the quadrature error as follows using the Peano kernel representation of the error for a first order quadrature formula,

[TABLE]

which gives for $Q_{l}$,

[TABLE]

We need to bound the integral $\int_{0}^{1}(1-s)\psi^{\prime}(\tau s)ds$. We first bound $\psi^{\prime}(s)$. For that, we need to bound $\|-A\mathrm{e}^{\tau(1-s)A}\widehat{\psi}(\tau s)\|$. We get

[TABLE]

Therefore

[TABLE]

We can now compute the error of the quadrature formula $Q_{l}$ and show (4.4). This follows from the inequality

[TABLE]

Which gives us, with (4.3), the desired result for any first order quadrature formula,

[TABLE]

which concludes the proof of Lemma 4.4. $\square$

Lemma 4.5**.**

Let $Q$ and $\psi$ be as in Lemma 4.4. We assume that $\widehat{\psi}$ is twice continuously differentiable and that $Q$ is a second order quadrature formula ($\sum_{k=1}^{m}c_{k}b_{k}=\frac{1}{2}$). Then

[TABLE]

If, in addition, $\widehat{\psi}(s)=\phi_{0}+\mathcal{O}(\tau)$ with $\phi_{0}\in D(A)$ and $A\phi_{0}=\mathcal{O}(1)$, then

[TABLE]

Using the notation (4.1), we remark that the term $A\mathcal{O}(\tau^{k})$, where $k=2$ in (4.5) and $k=3$ in (4.6), is a function of the form $A\phi$, where $\|\phi\|_{L^{p}(\Omega)}\leq C\tau^{k}$. In particular $\phi\in L^{p}(\Omega)$ is not in the domain of $A$ in general, and hence $A\phi\notin L^{p}(\Omega)$ in general and $A\phi\in W^{-2,p}(\Omega)$.

Proof.

The second derivative of $\psi$ is

[TABLE]

If $Q$ is a second order quadrature formula we write the quadrature error as follows using the Peano kernel representation of the error for a second order quadrature formula:

[TABLE]

It remains to estimate

[TABLE]

and

[TABLE]

We first bound $\|A\mathrm{e}^{\tau(1-s)A}\widehat{\psi}(\tau s)\|$. We get

[TABLE]

Then

[TABLE]

and

[TABLE]

This gives the following estimation for the integrals $\int_{0}^{1}\frac{(1-s)^{2}}{2}P_{i}(\tau s)ds$:

[TABLE]

[TABLE]

We show that $\|\sum_{i=1}^{m}b_{i}\int_{0}^{c_{i}}(c_{i}-s)P_{1}(\tau s)ds\|\leq\frac{C}{\tau}$.

[TABLE]

If $c_{i}=1$, we have

[TABLE]

If $c_{i}\neq 1$, then

[TABLE]

Therefore

[TABLE]

For the integral of $P_{2}$, we obtain

[TABLE]

which gives the desired bound for the error,

[TABLE]

If condition (4.2) is satisfied we can obtain a better bound for $\|A\mathrm{e}^{\tau(1-s)A}\widehat{\psi}(\tau s)\|$ and thus also for $\|P_{1}(\tau s)\|$,

[TABLE]

We obtain the following estimation for $P_{1}(\tau s)$,

[TABLE]

which gives us the estimation $\|\sum_{i=1}^{m}b_{i}\int_{0}^{c_{i}}(c_{i}-s)P_{1}(\tau s)ds\|\leq C$. Finally, we have the desired error bound

[TABLE]

which concludes the proof of Lemma 4.5. $\square$

Using Theorem 2.1, we can improve Lemma 4.5 as follows.

Lemma 4.6**.**

Under the hypotheses of Lemma 4.5, there exists $\alpha>0$ such that

[TABLE]

If additionally condition (4.2) is satisfied, then

[TABLE]

Proof.

We use Theorem 2.1, which states that for sufficiently small $\alpha>0$, $W^{1,p}(\Omega)$ is included in the domain of $(-A)^{\alpha}$, which does not involve any condition on the boundary. One then obtains

[TABLE]

and

[TABLE]

which gives

[TABLE]

instead of (4.7). Similarly, if condition (4.2) is satisfied, one obtains

[TABLE]

instead of (4.8), and this concludes the proof. $\square$

4.2 Order one error estimate for the five parts Strang splitting

In this section, we prove that the five parts modified splitting (3.2) is of global order one because this is needed in the proof of the global order of the method. We start to give two estimations for the local error. To perform our convergence analysis, we need an exact formula for (3.2). We expand each flow that appears in the Strang splitting,

[TABLE]

We obtain the following exact formula for the numerical flow:

[TABLE]

We define the local error at time $t_{n+1}$, $\delta_{n+1}$, as follows:

[TABLE]

Using the formula (2.4) of the exact solution and formula (4.2) of the numerical solution we obtain

[TABLE]

Since all the integrands are uniformly bounded on $[0,\tau]$, we obtain the following result, which states that the five parts modified Strang splitting (3.2) is locally of first order.

Lemma 4.7**.**

Under the assumption of section 2, the five parts modified Strang splitting (3.2) satisfies the following local error estimate:

[TABLE]

We prove the next local error estimate we use in the theorem for the global error.

Lemma 4.8**.**

Under the assumption of section 2, the five parts modified Strang splitting (3.2) satisfies the following local error estimate,

[TABLE]

Proof.

In formula (4.2) of the local error, we use the trapezoidal quadrature formula to approximate the integrals. By Lemma 4.5, the quadrature error made to approximate $\int_{0}^{\tau}\mathrm{e}^{(\tau-s)A}(q_{n}-f(u(t_{n}+s)))ds$ is equal to $A\mathcal{O}(\tau^{2})+\mathcal{O}(\tau^{2})$. We get

[TABLE]

Since $\phi^{f}_{\tau/2}(u(t_{n}))=u(t_{n})+\mathcal{O}(\tau)$, $f(u(t_{n}+\tau))=f(u(t_{n}))+\mathcal{O}(\tau)$ and $\tilde{v}_{n}=u(t_{n})+\mathcal{O}(\tau)$, and expending $\mathrm{e}^{\tau A}=Id+\tau A\varphi_{1}(\tau A)$, we obtain

[TABLE]

which concludes the proof. $\square$

Using Theorem 2.1 and Lemma 4.6, we can improve Lemma 4.8 as follows.

Lemma 4.9**.**

Under the assumption of section 2, the five parts modified Strang splitting (3.2) satisfies the following. There exists $\alpha>0$ such that

[TABLE]

Proof.

Using Lemma 4.6, we can obtain that the quadrature error made to approximate $\int_{0}^{\tau}\mathrm{e}^{(\tau-s)A}(q_{n}-f(u(t_{n}+s)))ds$ is equal to $(-A)^{1-{\alpha}}\mathcal{O}(\tau^{2})+\mathcal{O}(\tau^{2})$. We then use Theorem 2.1 to bound $\frac{\tau^{2}}{2}A\varphi_{1}(\tau A)\left(f(u(t_{n}))-q_{n}\right)+\frac{\tau^{2}}{2}A\varphi_{1}(\tau A)(q_{n}-f(u(t_{n}))$. We obtain

[TABLE]

This gives us the desired result. $\square$

Using the previous results for the local error, we can now prove the following order estimation for the global error.

Proposition 4.10**.**

Under the assumptions of section 2, the five parts modified Strang splitting (3.2) satisfies

[TABLE]

The constant C depends on T but is independent on $\tau$ and n.

Proof.

The global error is defined as $e_{n}=u_{n}-u(t_{n})$.

[TABLE]

Using the exact formula (4.2) for $\mathcal{S}_{\tau}u_{n}$ and $\mathcal{S}_{\tau}u(t_{n})$, we obtain for $e_{n+1}$,

[TABLE]

with

[TABLE]

Let us bound $E(u(t_{n}),u_{n})$. We use the Lipschitz continuity of $f$ and $\phi_{s}^{f}$. For the first integral in (4.2), we have

[TABLE]

For the second integral that appears in (4.2), we observe that

[TABLE]

We obtain

[TABLE]

Writing the exact formula for $\phi_{\tau}^{D+q_{n}}\circ\phi_{\frac{\tau}{2}}^{-q_{n}}\circ\phi_{\frac{\tau}{2}}^{f}$, we have

[TABLE]

Therefore

[TABLE]

The global error $e_{n}$ satisfies the following recursive formula:

[TABLE]

This gives us, thanks to to previous estimation for $\|E(u_{k},u(t_{k}))\|$ and $\|\delta_{k}\|$ and since $\|e_{0}\|=0$,

[TABLE]

Since, by Lemma 4.7, $\|\delta_{n}\|=\mathcal{O}(\tau)$ and using the parabolic smoothing property, we get

[TABLE]

We rearrange the second sum and observe that $nC\tau^{2}=C\tau$, which gives

[TABLE]

The second sum can be bounded as,

[TABLE]

Using the discrete Gronwall’s lemma we obtain the desired result:

[TABLE]

which concludes the proof. $\square$

Using Theorem 2.1 and Proposition 4.10 we are now in position to prove Proposition 4.2, which provides a first order estimation for the global error.

Proof of Proposition 4.2.

We use Lemma 4.6 to remove the term $\log(\tau)$ in the global error estimate. Indeed, we obtain,

[TABLE]

We then estimate

[TABLE]

which concludes the proof. $\square$

Remark 4.11**.**

In Lemma 4.12, to show that a function $q_{n}$ satisfying the boundary condition (3.3) or (3.4) satisfies the condition (4.2), we need to use Proposition 4.2, which we prove using Theorem 2.1. To prove $\|u_{n}-u(t_{n})\|\leq C\tau^{2}(1+|\log\tau|)$ without using Theorem 2.1, we need a weaker condition that $q_{n}$ must satisfy, for example, $f(u(t_{n}+s))-q_{n}=\phi_{n}+\mathcal{O}(\tau(1+|\log(\tau)|))$ with $\phi_{n}\in D(A)$ and $A\phi_{n}=\mathcal{O}(1)$, instead of (4.2). We can then prove, with the help of the bound, $\|u_{n}-u(t_{n})\|\leq C\tau(1+|\log\tau|)$, that a function satisfying (3.3) satisfies this new condition. We decide not to follow this approach as we think it simplifies our arguments to only have condition (4.2) throughout the paper.

4.3 Analysis of the corrector function

We show that the conditions (3.3) and (3.4) for $q_{n}$ are properly chosen, that is $\widehat{\psi}(s)=q_{n}-f(u(t_{n}+s))$ satisfies the hypothesis (4.2) when one of those conditions is satisfied. For that purpose, we use Proposition 4.2, which states that $u(t_{n})-u_{n}=\mathcal{O}(\tau)$. We stress that no condition on $q_{n}$ is required in the proof of Proposition 4.2. We first consider the boundary condition (3.3) for $q_{n}$.

Lemma 4.12**.**

Let $q_{n}$ be chosen such that (3.3) is satisfied,

[TABLE]

Let $\widehat{\psi}(s)=q_{n}-f(u(t_{n}+s))$. Then $\widehat{\psi}(s)=\phi_{n}+\mathcal{O}(\tau)$ with $\phi_{n}\in D(A)$ and $A\phi_{n}=\mathcal{O}(1)$.

Proof.

We observe that

[TABLE]

We obtain, with Proposition 4.2, that

[TABLE]

We define $\phi_{n}$ as follows,

[TABLE]

Since $B(q_{n}-(\frac{2}{\tau}\phi^{f}_{\frac{\tau}{2}}(u_{n})-\frac{2}{\tau}u_{n}))=0$, $\phi_{n}$ is in $D(A)$. Furthermore, a simple calculation yields

[TABLE]

Therefore $A\phi_{n}=\mathcal{O}(1)$. $\square$

We now consider the boundary condition (3.4) for the corrector functions.

Lemma 4.13**.**

Let $q_{n}$ be chosen such that (3.4) is satisfied,

[TABLE]

Let $\widehat{\psi}(s)=q_{n}-f(u(t_{n}+s))$. Then $\widehat{\psi}(s)=\phi_{n}+\mathcal{O}(\tau)$, with $\phi_{n}\in D(A)$ and $A\phi_{n}=\mathcal{O}(1)$.

Proof.

The proof is conducted by induction. Since $Bu_{0}=b_{0}$, we know by Lemma 4.12 that the result is true for $n=0$.

We assume that the statement is true for $n=0,\ldots,k-1$. Let us show that it is true for $n=k$. We write the exact formula for $\phi^{f}_{\frac{\tau}{2}}(u_{k})$ in function of $v_{k-1}$ and $q_{k}$:

[TABLE]

We then apply the midpoint quadrature formula to the integral and obtain an error of size $\mathcal{O}(\tau^{2})$ since $f$ is twice continuously differentiable:

[TABLE]

We observe that $\phi^{f}_{\frac{\tau}{2}}(v_{k-1}-\frac{\tau}{2}q_{k-1})=u_{k}$. Since by Proposition 4.2, $u_{k}=f(u(t_{k}))+\mathcal{O}(\tau)$ and since $f(u(t_{k}))=f(u(t_{k-1}))+\mathcal{O}(\tau)$, we obtain

[TABLE]

Since $Bv_{k-1}=b_{k}$, and since by hypothesis $f(u(t_{k-1}))-q_{k-1}=\phi_{k-1}+\mathcal{O}(\tau)$ with $\phi_{k-1}\in D(A)$ and $A\phi_{k-1}=\mathcal{O}(1)$, we obtain

[TABLE]

We can therefore decompose $q_{k}$ as $q_{k}=\tilde{q}_{k}+r_{k}$, where $B\tilde{q}_{k}=Bf(u(t_{k}))$ and $r_{k}=\mathcal{O}(\tau)$ and choose $\phi_{k}=\tilde{q}_{k}-f(u(t_{k}))$. Then $B\phi_{k}=0$ and a simple calculation yields $\|A\phi_{k}\|\leq C$. This concludes the proof. $\square$

4.4 Order two error estimate for the five parts modified Strang splitting

The following lemma is an estimation of the local error for the modified Stang splitting, which states that the five parts modified Strang splitting (3.2) is locally a method of second order.

Lemma 4.14**.**

Under the assumption of section 2, the five parts modified Strang splitting (3.2) satisfies

[TABLE]

Proof.

In the exact formula of $\delta_{n+1}$ (4.2), we use the left rectangle quadrature formula for the first and third integral and a right quadrature formula for the second integral. By Lemma 4.4, the quadrature error is $\mathcal{O}(\tau^{2})$. We get

[TABLE]

Expending $\mathrm{e}^{\tau A}$ and $u(t_{n})$ around $\tau$, that is, $\mathrm{e}^{\tau A}\frac{\tau}{2}(f(u(t_{n}))-q_{n})=\frac{\tau}{2}(f(u(t_{n+1}))-q_{n})+\mathcal{O}(\tau^{2})$, we obtain the following result:

[TABLE]

We use the exact formula for $\tilde{v}_{n}$ and $u(t_{n+1})$ and the Lipschitz continuity of $f$:

[TABLE]

Since the integrands are uniformly bounded on $[0,\tau]$, we get $\delta_{n+1}=\mathcal{O}(\tau^{2})$, which concludes the proof. $\square$

The following lemma gives the second local error estimates that we need in the proof for the global convergence of the method.

Lemma 4.15**.**

Under the assumption of section 2, the five parts modified Strang splitting (3.2) satisfies

[TABLE]

Proof.

We start as in the proof of Lemma 4.8, by using trapezoidal quadrature formulas to approximate the integrals in formula (4.2) of the local error. By Lemma 4.5, the quadrature error made to approximate $\int_{0}^{\tau}\mathrm{e}^{(\tau-s)A}(q_{n}-f(u(t_{n}+s)))ds$ is equal to $A\mathcal{O}(\tau^{3})+\mathcal{O}(\tau^{3})$. We get

[TABLE]

By Lemma 4.14, we have the following equality:

[TABLE]

We obtain

[TABLE]

We observe that $\frac{\tau}{2}\mathrm{e}^{\tau A}(f(u(t_{n}))-q_{n})+\frac{\tau}{2}(f(u(t_{n}+\tau))-q_{n})$ is the trapezoidal quadrature formula for $\int_{0}^{\tau}\mathrm{e}^{(\tau-s)A}(f(u(t_{n}+s)))-q_{n})ds$ and that $\mathrm{e}^{\tau A}\frac{\tau}{2}(f(\phi^{f}_{\tau/2}(u(t_{n})))-q_{n})+\frac{\tau}{2}(f(\tilde{v}_{n})-q_{n})$ is the trapezoidal quadrature formula for $\int_{0}^{\tau}\mathrm{e}^{(\tau-s)A}(f(\phi_{(\tau-s)/2}^{f}\circ\phi^{-q_{n}}_{s/2}\circ\phi_{s}^{D+q_{n}}\circ\phi^{-q_{n}}_{s/2}\circ\phi_{s/2}^{f}(u(t_{n})))-q_{n})ds$. Since $\phi_{(\tau-s)/2}^{f}\circ\phi^{-q_{n}}_{s/2}\circ\phi_{s}^{D+q_{n}}\circ\phi^{-q_{n}}_{s/2}\circ\phi_{s/2}^{f}(u(t_{n}))=u(t_{n})+\mathcal{O}(\tau)$, the second integrand satisfies the hypotheses of Lemma 4.5, by Lemmas 4.12 and 4.13. Applying Lemma 4.5, we therefore have a quadrature error of the form $A\mathcal{O}(\tau^{3})+\mathcal{O}(\tau^{3})$,

[TABLE]

Applying the midpoint quadrature method to both integrals, we obtain

[TABLE]

Using Lemma 4.14, the Lipschitz continuity of $f$, and the boundedness of $\mathrm{e}^{\frac{\tau}{2}A}$ we have the desired result,

[TABLE]

which concludes the proof. $\square$

We now improve Lemma 4.15 as follows using Lemma 4.6 instead of Lemma 4.5 for the quadrature error.

Lemma 4.16**.**

Under the assumption of section 2, the five parts modified Strang splitting (3.2) satisfies the following. There exists $\alpha>0$ such that

[TABLE]

Using the previous results for the local error, we can prove Proposition 4.3. It is the main global error estimate that we obtain without using Theorem 2.1.

Proof of Proposition 4.3.

The proof is similar to the proof of the Proposition 4.10, except for the bounds of the local errors $\delta_{n}$, since if (3.3) is satisfied, $\delta_{k}=A\mathcal{O}(\tau^{3})+\mathcal{O}(\tau^{3})$ and $\delta_{n}=\mathcal{O}(\tau^{2})$. Hence

[TABLE]

which concludes the proof. $\square$

Proof of Theorem 4.1.

We follow the proof of Proposition 4.2, using Lemma 4.16 to remove the $\log(\tau)$ term in the global error estimate of Proposition 4.3. $\square$

5 Numerical experiments

In this section we perform several numerical experiments in dimension $d=1,2$ to illustrate the performance of the five parts modified Strang splitting (3.2) when applied to diffusion problems with various nonlinearities. The norm we use to compute the numerical error is

[TABLE]

where we use the trapezoidal rule to approximate the $L^{2}(\Omega)$ norm. Considering the domain $(0,1)^{d}$ equipped with an uniform spatial grid with size $\Delta x$, we use the usual second order finite difference approximation of the Laplacian. We compute the diffusion flows in time exactly using the Matlab package described in [9] for approaching matrix exponentials and the $\varphi_{1}$ matrix function (2.3).

In the numerical experiments that follow, we denote the classical Strang splitting (1.2) by Strang, the modified Strang splitting (1.7) constructed in [2, 3] by StrangM3, the five parts modified Strang splitting (3.2) with correction (1.2) by StrangM5a and the five parts modified Strang splitting (3.2) with correction (1.3) by StrangM5b. In Table 1, we recall the number of evaluations of the diffusion flows (1.5, 1.8) and the number of evaluations of the source term flows (1.4, 1.6) required at each step of the algorithm for the four considered numerical methods. We also collect the total number of flows after $n$ steps, where we consider that the correction steps (3.1) are negligible. We recall that for Strang and StrangM5b, only one evaluation of the diffusion flow and one evaluation of the source term flow is needed after the first step, due to the semigroup property of the exact flows (see Remark 3.2).

A quadratic nonlinearity (see Figure 1)

First, we apply the splitting methods Strang (1.2), StrangM3 (1.7), StrangM5a (3.2, 3.3) and StrangM5b (3.2, 3.4) to a problem given in [3, Example 5.2]. We then change the nonlinearity to see how the splitting methods behave. The non linearities we consider are $f(u)=u^{2}$ and $f(u)=5u^{2}$. The case $f(u)=u^{2}$ is the one presented in [3, Example 5.2]. We perform the experiment with mixed boundary conditions, $u(0)=1$, $\partial_{n}u(1)=1$. We choose a smooth initial condition that satisfies the prescribed boundary conditions. We obtain the following equation with $m=1$ and $m=5$.

[TABLE]

The correction we use for StrangM3, given in [3] for $m=1$, is $q_{n}=m+2mxu_{n}(1)$. We use $500$ spatial points to discretize the interior of $\Omega$. We compute the solution at final time $T=0.1$. The chosen time steps are $\tau=0.02\cdot 2^{-k}$, $k=0,\ldots,6$. The reference solution is computed with the classical fourth order explicit Runge-Kutta method and a very small time step $\tau=0.02\cdot 2^{-14}\approx 10^{-6}$. In the splitting algorithms, we use the analytic formulas for computing the flows $\phi^{f}_{\frac{\tau}{2}}(u_{n})$ and $\phi^{f-q_{n}}_{\tau}(u_{n})$.

We observe in Figure 1 that the splittings StrangM3, StrangM5a and StrangM5b are all of order two compared to the classical splitting Strang. The methods StrangM5a and StrangM5b have a slightly worse constant of error compared to the splitting StrangM3 of [3] for $f(u)=u^{2}$. Recall that StrangM3 costs two diffusion flows per time step while Strang, StrangM5a and StrangM5b cost only one. For $f(u)=5u^{2}$, StrangM5a and StrangM5b become a lot more accurate.

A meteorology model with an integral source term (see Figure 2)

We apply the splitting methods Strang (1.2), StrangM3 (1.7), StrangM5a (3.2, 3.3) andStrangM5b (3.2, 3.4) to a problem presented in [10, Equation 3.1], where we replace the left Dirichlet boundary condition $u(0,t)=2(2-\sqrt{t})$ by $u(0,t)=2(2-t)$ to have a time continuously differentiable boundary condition. The considered differential equation is the following

[TABLE]

for $\Omega=[0,1]$. We choose $500$ points to discretize the interior of $[0,1]$. We then apply the splitting methods with different time steps $\tau=2\cdot 10^{-2}\cdot 2^{-k}$, $k=0,\ldots,6$. The reference solution is computed with the classical fourth order explicit Runge-Kutta method and a very small time step $\tau=0.02\cdot 2^{-14}\approx 10^{-6}$. To solve the integral, we use the trapezoidal quadrature formula with the $502$ nodes given by the space discretization. To solve $\partial_{t}u=f(u)$ and $\partial_{t}u=f(u)-q_{n}$, we apply five steps of the classical order four explicit Runge-Kutta method with a time step $\frac{\tau}{10}$. Note that we compute $\partial_{t}u=f(u)$ and $\partial_{t}u=f(u)-q_{n}$ with more precision than is needed in practice to emphasize on the error of the splitting itself. Instead, one can simply use one step of a second order Runge-Kutta method. We compute the solution at final time $T=0.1$.

We observe in the left picture of Figure 2 that the modified splitting StrangM3 given in [3] has a slightly better constant of error than StrangM5a and StrangM5b. In the right picture of Figure 2 however, we recover that the splitting StrangM5b requires less evaluations of the flows to obtain a given precision (see Table 1 and Remark 3.2). Since both flows have similar computational costs, the splitting StrangM5b is slightly more accurate for the same computational cost. Moreover the construction of the correction $q_{n}$ for the splitting StrangM5b requires less computational cost and is easier to implement than the correction $q_{n}$ of StrangM3 given by (1.9). Indeed, one needs then to evaluate $f$ on the boundary at each step of the algorithm.

Case of a stiff nonlinearity (see Figure 3)

In the next experiments, we compare Strang (1.2), StrangM3 (1.7), StrangM5a (3.2, 3.3) and StrangM5b (3.2, 3.4) when applied to a stiff problem in a two dimensional domain $\Omega=[0,1]\times[0,1]$. We choose the nonlinearity $f(u)=(1-M\sin(\pi x)\sin(\pi y))u^{2}$ where $M=1$ in the nonstiff case and $M=100$ in the stiff case. We impose Dirichlet boundary conditions on the left boundary and Neumann boundary conditions on the bottom, top and right boundaries. The boundary conditions are chosen to be consistent with the initial condition $u_{0}=\frac{\mathrm{e}^{x}+\mathrm{e}^{y}}{2}$. We obtain the following problem

[TABLE]

where $x,y\in\Omega$, $t\in[0,T]$ and with $M=1$ or $M=100$. The correction functions $q_{n}$ are constructed with a multigrid iterative smoothing algorithm, analogously to [3, Algorithm 1]. We recall that the method StrangM3 in [3] was proposed for nonstiff nonlinearities. Indeed note that $q_{n}$ is of size $\mathcal{O}(M)$ and becomes large for the stiff case $M=100$. For instance in dimension one, on the domain $[0,1]$ with $f(u)=(1-M\sin(\pi x))u^{2}$ and boundary conditions $u(0,t)=1$, $\partial_{n}u(1,t)=1$, one obtains $q_{n}=1+(M\pi u_{n}(1)^{2}+2u_{n}(1))x$. We use a mesh with size $\Delta x=\frac{1}{128}$ to discretize the interior of $\Omega$. We compute the solution at final time $T=0.1$. The chosen time steps are $\tau=0.1\cdot 2^{-k}$, $k=0,\ldots,8$. The reference solution is computed with the classical fourth order explicit Runge-Kutta method and a really small time step $\tau=0.1\cdot 2^{-14}\approx 10^{-6}$. In the splitting algorithms, we use the analytic formulas of $\phi^{f}_{\frac{\tau}{2}}(u_{n})$ and $\phi^{f-q_{n}}_{\tau}(u_{n})$ .

We observe in the top pictures of Figure 3 that StrangM5a and StrangM5b are slightly more accurate than the modified Strang splitting StrangM3 for $M=1$ and drastically more accurate for the stiff case $M=100$. Indeed, as explained above, in this later stiff case, the accuracy of the splitting StrangM3 deteriorates due to the large correction involved. Moreover, in the bottom pictures of Figure 3, we observe that for the same number of evaluations of the diffusion flows, StrangM5b is more accurate than all the other methods also for $M=1$. Indeed, here, the diffusion flows dominates the total cost and we recall that StrangM5b only requires one evaluation of the diffusion flow and one evaluation of the source term flow at each step (see Table 1 and Remark 3.2), analogously to the classical Strang splitting Strang.

Acknowledgements.

The authors would like to thank Christophe Besse and Lukas Einkemmer for helpful discussions. This work was partially supported by the Swiss National Science Foundation, grants No. 200021_162404, No. 200020_178752, No. 200020_144313/1, No. 200020_184614, No. 200021_162404 and No. 200020_178752.