A high-order exponential integrator for nonlinear parabolic equations with nonsmooth initial data

Let $A$ be the generator of a bounded analytic semigroup on a Banach space $X$, with domain $D(A)\subset X$, and consider the abstract semilinear initial-value problem

where $u_{0}\in X$ and $f:[0,\infty)\times X\rightarrow X$ is a smooth (locally Lipschitz continuous) function. A function $u\in C([0,T];X)$ is called a mild solution of (1) if it satisfies the integral equation

where $e^{tA}$ denotes the semigroup generated by the operator $A$.

In the linear case $f(t,u)\equiv f(t)$, time discretization of (1) by a $k^{\rm th}$-order Runge–Kutta method satisfies the following error estimate:

where $\tau$ denotes the stepsize of time discretization; see Lubich-Ostermann-1996 ; Thomee2006 . In particular, for a nonsmooth initial value $u_{0}\in X$, the methods have $k^{\rm th}$-order accuracy when $t_{n}$ is not close to zero. This result also holds for implicit backward difference formulae (BDF), exponential integrators Hochbruck-Ostermann-2010 ; Koskela-Ostermann-2013 , and fractional-order evolution equations (Jin-Li-Zhou-2017, , Remark 2.6).

However, such high-order convergence as (3) does not hold when the source function $f(t,u)$ is nonlinear with respect to $u$. A counter example constructed in Crouzeix-Thomee-1987 shows that a $k^{\rm th}$-order Runge–Kutta method normally has only first-order convergence for a general nonsmooth initial data $u_{0}\in X$, i.e.,

Similarly, semi-implicit Runge–Kutta methods also suffer from this barrier of convergence rate Ostermann-Thalhammer-2000 . For nonlinear problems with nonsmooth initial data, existing error estimates for exponential integrators also yield only first-order convergence (see Hochbruck-Ostermann-2005 ; Mukam-Tambue-2018 ; Ostermann-Thalhammer-2000 )

No method has been proved to have high-order convergence for semilinear parabolic equations with general nonsmooth initial data $u_{0}\in X$.

Of course, if the initial value is sufficiently smooth and satisfies certain compatibility conditions, e.g., $u_{0}\in D(A^{k})$, then $O(\tau^{k})$ convergence can be achieved uniformly for $t_{n}\in[0,T]$ for the nonlinear problem (1). This has been proved for most time discretization methods, including Runge–Kutta methods Crouzeix-Thomee-1987 , implicit A$(\alpha)$-stable multistep methods Lubich-1990 , implicit–explicit BDF methods Akrivis-2015 ; Akrivis-Lubich-2015 , splitting methods Einkemmer-Ostermann-2015 ; Einkemmer-Ostermann-2016 ; Hansen-Kramer-Ostermann-2012 and several types of exponential integrators Calvo-Palencia-2006 ; Hochbruck-Ostermann-2005 ; Hochbruck-Ostermann-2010 ; Ostermann-Thalhammer-Wright-2006 . Extension to quasi-linear parabolic problems has also been done; see Gonzalez-Thalhammer-2007 ; Gonzalez-Thalhammer-2015 ; Gonzalez-Thalhammer-2016 ; Hochbruck-Ostermann-2011 ; Hochbruck-Ostermann-Schweitzer-2009 ; Lubich-Ostermann-1995 . The error estimates presented in these articles do not apply to nonsmooth initial data.

Due to the presence of the factor $t_{n}^{-1}$, the first-order convergence of Runge–Kutta methods cannot be improved by using variable stepsizes. However, compared with other time discretization methods, exponential integrators were proved to have an error bound of $O(\tau)$ uniformly for $t_{n}\in[0,T]$ for nonsmooth initial data $u_{0}\in X$, without the factor $t_{n}^{-1}$ appearing in the error estimates for other methods (see Hochbruck-Ostermann-2005 ; Mukam-Tambue-2018 ; Ostermann-Thalhammer-2000 ). This uniform convergence motivates us to consider the possibility of constructing high-order exponential integrators with variable stepsizes.

In this paper, we propose a variable stepsize exponential $k$-step integrator for (1) with general nonsmooth initial data $u_{0}\in X$, by choosing

where $\tau_{n}=t_{n}-t_{n-1}$ denotes the $n^{\rm th}$ stepsize in the partition $0=t_{0}<t_{1}<\dots<t_{N}=T$, and $\tau$ the maximal stepsize. For the convenience of implementation, we also integrate in the numerical method (and the error analysis) an algorithm for approximating the exponential integrator by using the contour integral techniques developed in Lopez-Fernandez-2010 ; LFPS-2006 ; Schadle-LF-Lubich-2006 .

The proposed variable stepsize method, with contour integral approximation of the exponential integrator, can achieve $k^{\rm th}$-order accuracy in approximating a mild solution of (1), i.e.,

This is the first high-order convergence result in approximating semilinear parabolic equations with nonsmooth initial data (without any regularity in addition to $u_{0}\in X$). In view of the result (4) for the Runge–Kutta methods, the convergence result (7) shows the superiority of the variable stepsize exponential integrator for problems with nonsmooth initial data. The approximation of exponential integrator would require $O(\ln(\tau^{-1}))$ parallel solutions of linear equations, and there are $N=O(\tau^{-1})$ time levels by using the stepsize in (6). Therefore, the total computational cost is $O(\tau^{-1}\ln(\tau^{-1}))$ for an accuracy of $O(\tau^{k})$.

For rigorous analysis without extra regularity assumptions on the solution, we assume that the nonlinear source function satisfies the following estimates:

where $\|\cdot\|$ denotes the norm of $X$, $C_{f,u,v}$ is a constant depending on $f$, $\|u\|$ and $\|v\|$; similarly, $C_{f,u,\ell}$ is a constant depending on $f$, $\|u(t)\|$, $\ell$, and the summation above extends over all possible positive integers $m_{1},\dots,m_{j}$ satisfying $m_{1}+\dots+m_{j}\leq\ell$ for a given $j$.

Assumptions (8)–(9) are naturally satisfied by a general smooth function $f:{\mathbb{R}}\rightarrow{\mathbb{R}}$ in a semilinear parabolic partial differential equation (PDE)

In this case, the Dirichlet Laplacian $\Delta$ generates a bounded analytic semigroup on $X=C_{0}(\overline{\Omega})$, the space of continuous functions on $\overline{\Omega}$ which equal zero on the boundary $\partial\Omega$; see Ouhabaz1995 . Furthermore, the smooth function $f$ naturally extends to a function of $u\in C_{0}(\overline{\Omega})$, satisfying

and

Obviously, all these time derivatives of $f(u(x,t))$ satisfy (9). Hence, the semilinear parabolic PDE (10) with a general smooth function $f:{\mathbb{R}}\rightarrow{\mathbb{R}}$ is an example of the abstract problem (1) satisfying assumptions (8)–(9).

Assumption (8) is the same as the local Lipschitz continuity assumption used in Calvo-Palencia-2006 and Hochbruck-Ostermann-2005 . In Calvo-Palencia-2006 , authors proved high-order convergence with an addition assumption that the solution is in $C^{k}([0,T],X)$, which is not satisfied when the initial data is nonsmooth, i.e., $u_{0}\in X$ instead of $D(A^{k})$. In Hochbruck-Ostermann-2005 , authors proved high-order convergence of exponential integrators with an additional assumption that $\partial_{t}^{k}f(u(t))$ is uniformly bounded for $t\in[0,T]$, which is also not satisfied when the initial data is nonsmooth. These additional assumptions in Calvo-Palencia-2006 ; Hochbruck-Ostermann-2005 are replaced by (9) in this paper, which is used to prove the weighted estimates

which allow the solution to be nonsmooth at $t=0$. These weighted estimates are used to prove high-order convergence of the exponential integrator in this paper.

Both the regularity analysis and the error analysis in this paper can be similarly extended to semilinear parabolic equations with smoothly varying time-dependent operators. However, the extension to quasilinear parabolic equations with nonsmooth initial data is still not obvious.

We denote by $\widehat{g}(z):=\int_{0}^{\infty}e^{-zt}g(t){\mathrm{d}}t$ the Laplace transform of a given function $g$. Then we let $g(t):=f(t,u(t))$ and take the Laplace transform of (2) in time. This yields

Since $A$ generates a bounded analytic semigroup on $X$, there exists an angle $\phi\in(0,\frac{\pi}{2})$ such that the operator $(z-A)^{-1}$ is analytic with respect to $z$ in the sector

In order to use the established contour integral techniques of Lopez-Fernandez-2010 ; LFPS-2006 , we take inverse Laplace transform of (11) along the contour

where $\alpha=\frac{\pi}{4}-\frac{\phi}{2}$ and $\lambda$ is to be determined. Then we have

Similarly, by considering $t_{n-1}$ as the initial time, the solution at $t=t_{n}$ can be written as

where $g_{n}(s)=f(t_{n-1}+s,u(t_{n-1}+s))$.

In (Lopez-Fernandez-2010, , Theorem 1) the authors proved that, by choosing the parameter

with

there are quadrature nodes and weights on the contour $\Gamma_{\lambda_{n}}$,

such that (12) can be approximated by a quadrature

with an error of $O(e^{-K/C})$.

Therefore, if $u^{(\tau)}=(u_{n})_{n=0}^{N}$ denotes the numerical approximation of $(u(t_{n}))_{n=1}^{N}$, then we approximate the source function $f(t,u(t))$ by an extrapolation polynomial of degree $k-1$:

where $L_{j}(t)$ is the unique polynomial of degree $k-1$ such that

For $n\geq k+1$ and given numerical solutions $u_{n-j}$, $j=1,\dots,k$, we denote

and compute

The numerical solution at the starting $k$ steps can be computed by using the exponential Euler method

Since the stepsize choice in (6) implies $\tau_{n}=O(\tau^{\frac{1}{1-\beta}})=O(\tau^{k})$ for the starting $k$ steps, the exponential Euler scheme (2) can keep the errors of numerical solutions within $O(\tau^{k})$ at the starting $k$ steps.

The main result of this paper is the following theorem.

The proof consists of two parts. In section 3.1, we prove the regularity of the solution $u$. By using the regularity result, we estimate errors of numerical solutions in section 3.2.