Incomplete iterative solution of subdiffusion

Abstract

In this work, we develop an efficient incomplete iterative scheme for the numerical solution of the subdiffusion model involving a Caputo derivative of order $\alpha \in (0,1)$ in time. It is based on piecewise linear Galerkin finite element method in space and backward Euler convolution quadrature in time and solves one linear algebraic system inexactly by an iterative algorithm at each time step. We present theoretical results for both smooth and nonsmooth solutions, using novel weighted estimates of the time-stepping scheme. The analysis indicates that with the number of iterations at each time level chosen properly, the error estimates are nearly identical with that for the exact linear solver, and the theoretical findings provide guidelines on the choice. Illustrative numerical results are presented to complement the theoretical analysis.

Appendices

Appendix A: Basic estimates

Lemma 11

For $\beta ,\gamma \ge 0$, there holds

$$\begin{aligned} \sum _{i=1}^n(n+1-i)^{-\beta }i^{-\gamma } \le \left\{ \begin{array}{ll} cn^{\max (1-\gamma ,0)-\beta }, &{}\quad 0\le \beta <1,\gamma \ne 1,\\ cn^{-\beta }\ln (1+n), &{}\quad 0\le \beta \le 1,\gamma =1,\\ c{n^{-\min (\beta ,\gamma )}}, &{}\quad \beta>1, \gamma >1. \end{array}\right. \end{aligned}$$

Proof

We denote by $[\cdot ]$ the integral part of a real number. Then

$$\begin{aligned} \sum _{i=1}^n(n+1-i)^{-\beta }i^{-\gamma } = \sum _{i=1}^{[{n/2}]}(n+1-i)^{-\beta }i^{-\gamma } + \sum _{i=[{n/2}]+1}^n(n+1-i)^{-\beta }i^{-\gamma }:=\mathrm{I}+\mathrm{II}. \end{aligned}$$

Then, by the trivial inequalities: for $1\le i\le [{n/2}]$, there holds $(n+1-i)^{-\beta }\le cn^{-\beta }$ and for $[{n/2}]+1\le i\le n$, there holds $i^{-\gamma }\le cn^{-\gamma }$, we deduce

$$\begin{aligned} \mathrm{I} \le cn^{-\beta } \sum _{i=1}^{[{n/2}]}i^{-\gamma }\quad \text{ and }\quad \mathrm{II}\le cn^{-\gamma } \sum _{i=[{n/2}]+1}^n(n+1-i)^{-\beta }. \end{aligned}$$

Simple computation gives $\sum _{i=1}^ji^{-\gamma }\le cj^{\max (1-\gamma ,0)}$ if $\gamma \ne 1$ and $\sum _{i=1}^ji^{-1}\le c\ln (j+1)$. Combining these estimates yields the desired assertion.$\square$

Next we give an upper bound on the CQ weights $b_j^{(\alpha )}$.

Lemma 12

For the weights $b_j^{(\alpha )}$, $|b_j^{(\alpha )}| \le e^{2\alpha }(j+1)^{-\alpha -1}$.

Proof

The weight $b_j^{(\alpha )}$ is given by $b_0^{(\alpha )}=1$ and $b_j^{(\alpha )}= -\varPi _{\ell =1}^j(1-\frac{1+\alpha }{\ell })$ for any $j\ge 1$. Note the elementary inequality $\ln (1-x)\le -x$ for any $x\in (0,1)$, and the estimate $\sum _{\ell =1}^j\ell ^{-1}\ge \int _1^{j+1}s^{-1}\mathrm{d}s = \ln (j+1).$ Since $\ln \alpha = \ln (1-(1-\alpha ))\le \alpha -1$, for any $j\ge 1$,

$$\begin{aligned} \ln |b_j^{(\alpha )}|&= \ln \alpha +\sum _{\ell =2}^j \ln \left( 1-\frac{1+\alpha }{\ell }\right) \le \ln \alpha - \sum _{\ell =2}^j\frac{1+\alpha }{\ell } \\&= \ln \alpha + (1+\alpha ) - \sum _{\ell =1}^j\frac{1+\alpha }{\ell } \le 2\alpha -(1+\alpha ) \ln (j+1). \end{aligned}$$

This completes the proof of the lemma. $\square$

Appendix B: Proof of Lemmas 4 and 5

In this part, we provide the proof of Lemmas 4 and 5. The proof of Corollary 1 is identical with that for Lemma 5 and thus it is omitted. The proof relies on the discrete Laplace transform, and the following two well-known estimates

$$\begin{aligned} \quad c_1 |z|&\le |\delta _\tau (e^{-z\tau })| \le c_2|z| \quad \forall \; z\in \varGamma _{\theta ,\delta }^\tau , \end{aligned}$$

(B.1)

$$\begin{aligned} |\delta _\tau (e^{-z\tau })|&\le |z| \sum _{k=1}^\infty \frac{|z\tau |^{k-1}}{k!} \le |z|e^{|z|\tau }, \quad \forall \; z\in \varSigma _\theta , \end{aligned}$$

(B.2)

and the resolvent estimate: for any $\theta \in (\pi /2,\pi )$,

$$\begin{aligned} \Vert (z+A_h)^{-1}\Vert \le c|z|^{-1},\quad \forall \; z\in \varSigma _\theta . \end{aligned}$$

(B.3)

Now we can give the proof of Lemma 4.

Proof of Lemma 4

By Laplace transform, $w_h(t_n)={\bar{\partial }}_\tau ^2 u_h(t_n)$ is given by

$$\begin{aligned} w_h(t_n) = \frac{1}{2\pi i} \int _{\varGamma _{\theta ,\delta }}\delta _\tau (e^{-z\tau })^2e^{zt_n}K(z)v_h \mathrm{d}z,\quad \text{ with } K(z) = z^{\alpha -1}(z^\alpha +A_h)^{-1}. \end{aligned}$$

We split the contour $\varGamma _{\theta ,\delta }$ into $\varGamma _{\theta ,\delta }^\tau$ and $\varGamma _{\theta ,\delta } {\setminus }\varGamma _{\theta ,\delta }^\tau$, and denote the corresponding integral by $\mathrm{I}$ and $\mathrm{II}$, respectively. We discuss the cases $v\in L^2(\varOmega )$ and $v\in D(A)$, separately.

Case (i): $v\in L^2(\varOmega )$. By (B.1) and (B.3), $\Vert K(z)\Vert \le c$ for $z\in \varGamma _{\theta ,\delta }^\tau$. Then choosing $\delta =c/t_n$ in $\varGamma _{\theta ,\delta }^\tau$ gives

$$\begin{aligned} \Vert \mathrm{I}\Vert _{L^2(\varOmega )} \le c \Vert v_h \Vert _{L^2(\varOmega )} \Big (\int _{\frac{c}{t_n}}^{\frac{\pi \sin \theta }{\tau }} \rho e^{t_n\rho \cos \theta } \,\mathrm{d}\rho + \int _{-\theta }^\theta t_n^{-2} \,\mathrm{d}\varphi \Big ) \le c t_n^{-2} \Vert v_h \Vert _{L^2(\varOmega )}. \end{aligned}$$

For any $z=\rho e^{\pm \mathrm {i}\theta }\in \varGamma _{\theta ,\delta }{\setminus }\varGamma _{\theta ,\delta }^\tau$, by the estimates (B.2) and (B.3), $\Vert K(z)\Vert \le ce^{2\rho \tau }.$ By choosing $\theta \in (\pi /2,\pi )$ sufficiently close to $\pi$, we deduce

$$\begin{aligned} \Vert \mathrm{II}\Vert _{L^2(\varOmega )} \le c \Vert v_h\Vert _{L^2(\varOmega )} \int _{\frac{\pi \sin \theta }{\tau }}^\infty e^{\rho (\cos \theta t_n+2\tau )} \rho \,\mathrm{d}\rho \le ct_n^{-2} \Vert v_h\Vert _{L^2(\varOmega )}. \end{aligned}$$

Thus, $\Vert {\bar{\partial }}_\tau ^2u_h(t_n)\Vert \le c t_n^{-2}\Vert v_h\Vert _{L^2(\varOmega )}.$ Next, by the identity $A_h(z^\alpha +A_h)^{-1}=I-z^\alpha (z^\alpha +A_h)$ and (B.3), $\Vert A_hK(z)\Vert \le |z|^{\alpha -1}$ for $z\in \varSigma _\theta$. Then repeating the argument gives

$$\begin{aligned} \tau ^\alpha \Vert A_h {\bar{\partial }}_\tau ^2 u_h(t_n)\Vert \le c \tau ^\alpha t_n^{-2-\alpha } \Vert v_h\Vert _{L^2(\varOmega )} \le ct_n^{-2} \Vert v_h \Vert _{L^2(\varOmega )}. \end{aligned}$$

Then the assertion for the case $v\in L^2(\varOmega )$ follows from the triangle inequality.

Case (ii): $v\in D(A)$. Simple computation gives the identity $K(z)v_h=z^{\alpha -1}(z^\alpha +A_h)^{-1}v_h = z^{-1}v_h - z^{-\alpha }(z^\alpha +A_h)^{-1}A_hv_h$. Thus, we have

$$\begin{aligned} w_h(t_n) = - \frac{1}{2\pi i} \int _{\varGamma _{\theta ,\delta }}e^{zt_n} \delta _\tau (e^{-z\tau })^2 z^{-\alpha }K(z)A_h v_h \mathrm{d}z, \end{aligned}$$

in which we split the contour $\varGamma _{\theta ,\delta }$ into $\varGamma _{\theta ,\delta }^\tau$ and $\varGamma _{\theta ,\delta }{\setminus } \varGamma _{\theta ,\delta }^\tau$, and accordingly the integral. Then the rest of the proof follows from the estimates (B.1), (B.2) and (B.3) as before.$\square$

Last, we prove Lemma 5.

Proof of Lemma 5

By Laplace transform and its discrete analogue, we have

$$\begin{aligned} \partial _t^\alpha y_h(t_n)-{\bar{\partial }}_\tau ^\alpha y_h(t_n)&= \frac{1}{2\pi \mathrm {i}} \int _{\varGamma _{\theta ,\delta }^\tau } e^{zt_n} K(z) A_h v_h\,\mathrm{d}z + \frac{1}{2\pi \mathrm {i}} \int _{\varGamma _{\theta ,\delta }{\setminus }\varGamma _{\theta ,\delta }^\tau } e^{zt_n} K(z) A_h v_h\,\mathrm{d}z\\&:=\mathrm{I}+\mathrm{II}, \end{aligned}$$

with $K(z)=(\delta _\tau (e^{-z\tau })^{\alpha }-z^\alpha )z^{ -1} (z^\alpha +A_h)^{-1}$. Recall the following estimate:

$$\begin{aligned} | \delta _\tau (e^{-z\tau })^{\alpha }-z^\alpha | \le c \tau z^{1+\alpha },\quad \forall \; z\in \varGamma _{\theta ,\delta }^\tau . \end{aligned}$$

(B.4)

Then by choosing $\delta =c/t_n$ in the contour $\varGamma _{\theta , \delta }^\tau$ and the resolvent estimate (B.3), we obtain

$$\begin{aligned} \Vert \mathrm{I}\Vert _{L^2(\varOmega )}&\le c \tau \Vert A_h v_h \Vert _{L^2(\varOmega )} \Big (\int _{\frac{c}{t_n}}^{\frac{\pi \sin \theta }{\tau }} e^{-c\rho t_n} \,\mathrm{d}\rho + \int _{-\theta }^\theta ct_n^{-1} \,\mathrm{d}\varphi \Big ) \le c\tau t_n^{-1} \Vert A v \Vert _{L^2(\varOmega )}. \end{aligned}$$

Further, by (B.2), for any $z=\rho e^{\pm \mathrm {i}\theta }\in \varGamma _{\theta ,\delta }{\setminus }\varGamma _{\theta ,\delta }^\tau$ and choosing $\theta \in (\pi /2,\pi )$ close to $\pi$,

$$\begin{aligned} |e^{zt_n} (\delta _\tau (e^{-z\tau })^{\alpha }-z^\alpha )z^{ -1}|&\le e^{t_n\rho \cos \theta }(c|z|^\alpha e^{\alpha \rho \tau }+|z|^\alpha )|z|^{-1}\le c|z|^{\alpha -1}e^{-c\rho t_n}. \end{aligned}$$

Then we deduce

$$\begin{aligned} \Vert \mathrm{II}\Vert _{L^2(\varOmega )}&\le c \Vert A_hv_h\Vert _{L^2(\varOmega )} \int _{\frac{\pi \sin \theta }{\tau }}^\infty e^{-c\rho t_n} \rho ^{-1}\,\mathrm{d}\rho \le c\tau t_n^{-1} \Vert Av\Vert _{L^2(\varOmega )}. \end{aligned}$$

Thus, we show the assertion for $\beta =0$. For the case $\beta =1$, the identity $A_h(z^\alpha +A_h)^{-1}=I- z^\alpha (z^\alpha +A_h)$, (B.3) and (B.4) give

$$\begin{aligned} \Vert A_h\mathrm{I}\Vert _{L^2(\varOmega )}&\le c \tau \Vert A_h v_h\Vert _{L^2(\varOmega )} \Big (\int _{\frac{c}{t_n}}^{\frac{\pi \sin \theta }{\tau }} e^{-c\rho t_n}\rho ^\alpha \,\mathrm{d}\rho + \int _{-\theta }^\theta ct_n^{-1-\alpha } \,\mathrm{d}\varphi \Big )\\&\le c\tau t_n^{-1-\alpha } \Vert A v\Vert _{L^2(\varOmega )}, \end{aligned}$$

and the bound on $\Vert A_h\mathrm{II}\Vert _{L^2(\varOmega )}$ follows analogously, completing the proof for $\beta =1$. Then the case $\beta \in (0,1)$ follows by interpolation. This shows part (i). The proof of part (ii) is similar and applies the $L^2(\varOmega )$ stability of $P_h$, and hence the detail is omitted. $\square$