Two Time-Stepping Schemes for Sub-Diffusion Equations with Singular Source Terms

Singular source terms in sub-diffusion equations may lead to the unboundedness of solutions, which will bring a severe reduction of convergence order of existing time-stepping schemes. In this work, we propose two efficient time-stepping schemes for solving sub-diffusion equations with a class of source terms mildly singular in time. One discretization is based on the Grünwald-Letnikov and backward Euler methods. First-order error estimate with respect to time is rigorously established for singular source terms and nonsmooth initial data. The other scheme derived from the second-order backward differentiation formula (BDF) is proved to possess second-order accuracy in time. Further, piecewise linear finite element and lumped mass finite element discretizations in space are applied and analyzed rigorously. Numerical investigations confirm our theoretical results.

Enquiries about data availability should be directed to the authors.

Authors and Affiliations

1. Institute of Mathematics, Hebei University of Technology, Tianjin, 300401, China

   Han Zhou

2. Center for Applied Mathematics, Tianjin University, Tianjin, 300072, China

   Wenyi Tian

Corresponding author

Correspondence to Wenyi Tian.

Competing interests

The authors have not disclosed any competing interests.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Han Zhou was partially supported by the National Natural Science Foundation of China (No. 11901151).

Wenyi Tian was partially supported by the National Natural Science Foundation of China (Nos. 11701416 and 12071343).

Appendix

A Proofs of Theorems 3.1 and 3.3

1.1 A.1 Proof of Theorem 3.1

Proof of Theorem 3.1

We obtain from (2.5) and (3.5) that

where

By the similar estimate in the proof of Theorem 2.1 in [23], it holds that \(\Vert {\hat{G}}_h(s)\Vert \le Ch^2.\) Then for \(\varepsilon =t^{-1}\), together with the condition \(\Vert {\hat{f}}(s)\Vert \le c|s|^{-\mu -1}\) in Assumption 1, we have

which completes the proof. \(\square \)

1.2 A.2 Proof of Theorem 3.3

For the convergence analysis of the lumped mass FEM, the quadrature error operator \(Q_h : X_h\rightarrow X_h\) was introduced in [4], which is defined by

It was analyzed in [4] that the quadrature error operator \(Q_h\) due to mass lumping satisfies the following estimates.

Lemma A.1

([4]) Let the operators \({\bar{\Delta }}_h\) and \(Q_h\) defined by (3.10) and (A.4), respectively. Then it holds that

Furthermore, if the mesh is symmetric, then it satisfies

With the quadrature error operator defined by (A.4), it yields from (3.1) and (3.8) that the error \(e_h(t)=u_h(t)-{\bar{u}}_h(t)\) satisfies

Taking the Laplace transform on (A.7) implies that

Since \(\Delta _h\) satisfies the resolvent estimate \(\Vert \big (s-\Delta _h\big )^{-1}\Vert \le M|s|^{-1}\), it is derived from (3.4) that

In addition, the operator \({\bar{\Delta }}_h\) defined by (3.10) also satisfies the resolvent estimate, then it follows from (A.5) that

Therefore, by (A.9), (A.10) and the Cauchy’s theorem, the inverse Laplace transform on (A.8) implies that the error \(e_{h}(t)\) for \(t>0\) can be represented by an integral over \(\Gamma _{\varepsilon }^{\theta }\cup S_{\varepsilon }\) as follows

Now it is ready to establish the error estimate for the lumped mass finite element scheme (3.8). The error is splitted into \(u(t)-{\bar{u}}_h(t)=u(t)-u_h(t)+e_h(t)\) with \(u_h(t)\) being the solution of the standard Galerkin finite element scheme in (3.1). Since the error \(\Vert u(t)-u_h(t)\Vert \) is estimated in Theorems 3.1 and 3.2 , then we next focus on the estimate of \(\Vert e_h(t)\Vert \).

Proof of Theorem 3.3

As \((s^{\alpha }-{\bar{\Delta }}_{h})^{-1}{\bar{\Delta }}_h=s^{\alpha } (s^{\alpha }-{\bar{\Delta }}_{h})^{-1}-I\) with I being the identity operator, it follows from the resolvent estimate of \({\bar{\Delta }}_{h}\) that \(\Vert (s^{\alpha }-{\bar{\Delta }}_{h})^{-1}{\bar{\Delta }}_h\Vert \le M+1\). Then we derive from (A.11) and (A.9) with \(u^0(x)\equiv 0\) that

By the trivial inequality \(\Vert \psi \Vert \le c\Vert \nabla \psi \Vert \) for \(\psi \in X_h\) and the estimate (A.5), it yields

Together with (A.9), we obtain the estimate (3.12) by (3.6) and the following argument

If the quadrature error operator \(Q_h\) satisfies (A.6), then it follows the estimate (3.13) from (3.6) and

This completes the proof. \(\square \)

B Four Lemmas

We provide four preliminary lemmas for the error analysis of the GLBE scheme (4.9) and the FBDF22 scheme (4.33).

Lemma B.1

If \(z\in \Sigma _{\pi /2}\), then \((1-e^{-z})^{\alpha }\in \Sigma _{\alpha \pi /2}\). Otherwise if \(z\in \Sigma _{\theta }\setminus \Sigma _{\pi /2}\) and its imaginary part satisfying \(|\Im (z)|\le \pi \) for \(\theta \in (\pi /2,\pi )\), then \((1-e^{-z})^{\alpha }\in \Sigma _{\alpha \theta }\).

Proof

If \(z=x+iy\in \Sigma _{\pi /2}\), then \(x>0\). This yields that the real part of \((1-e^{-z})\) satisfies

for all \(y\in \mathbb {R}\). Then we obtain \((1-e^{-z})\in \Sigma _{\pi /2}\) and \((1-e^{-z})^{\alpha }\in \Sigma _{\alpha \pi /2}\).

Otherwise if \(z\in \Sigma _{\theta }\setminus \Sigma _{\pi /2}\), then \(x\le 0\) and \(x\tan \theta \le | y|\le \pi \). It suffices to consider the case \(\Re (1-e^{-z})<0\). We define

From

and

it follows that \(f(x,y)\ge f(x,x\tan \theta )\) for \(x\le 0\). Taking the derivative of \({\tilde{f}}(x,{\theta }):=f(x,x\tan \theta )\) with respect to x arrives at

From \(\partial g/\partial x=e^{x}\sin (x\tan \theta )(1+\tan ^{2}\theta )\ge 0\), it follows that \(g(x,\theta )\le g(0,\theta )=0\). This leads to \({\tilde{f}}(x,{\theta })\ge {\tilde{f}}(0,{\theta })=-\tan \theta \) for any \(\theta \in (\pi /2,\pi )\). Therefore, we obtain \((1-e^{-z})\in \Sigma _{\theta }\) and the desired result. \(\square \)

Lemma B.2

If \(z\in {\mathbb {C}}\) and \(|z|\le r\) for finite \(r>0\), then

and

hold for \(0<\beta \le 1\), where C denotes a generic constant dependent on the radius r.

Proof

Using Taylor’s expansion of \(e^{-z}\) at 0, we derive

where the last inequality follows from the relation that \(\frac{e^{|z|}-1}{|z|}=\sum \limits _{n=1}^{+\infty }\frac{|z|^{n-1}}{n!}\le \sum \limits _{n=1}^{+\infty }\frac{r^{n-1}}{n!}=\frac{e^{r}-1}{r}\) for \(|z|\le r\). Similarly, we get

Furthermore, for \(0<\beta <1\), by the result in [16, 23], we have

which completes the proof. \(\square \)

Lemma B.3

Let \(|z|\le r\) for finite r. Then it holds that

for \(0<\beta \le 1\), where \(C=C(r)\).

Proof

Using Taylor’s expansions of \(e^{-z}\) and \(e^{-2z}\) at \(z=0\) arrives at

and

Then the second estimate of (B.3) is derived from the approach proposed in [16, 23]. \(\square \)

Lemma B.4

If \(z\in \Sigma _{\pi /2}\), then \((\frac{3}{2}-2e^{-z}+\frac{1}{2}e^{-2z})^{\alpha }\in \Sigma _{\alpha \pi /2}\). Otherwise if \(z\in \Sigma _{\theta }\setminus \Sigma _{\pi /2}\) for \(\theta \in (\pi /2, {\tilde{\theta }})\) with some \({\tilde{\theta }}\in (\pi /2,\pi )\) and \(|\Im z|\le \pi \), then there corresponds some \(\vartheta \in (\pi /2,\pi )\) such that \((\frac{3}{2}-2e^{-z}+\frac{1}{2}e^{-2z})^{\alpha }\in \Sigma _{\alpha \vartheta }\).

Proof

For \(z=x+iy\in \Sigma _{\pi /2}\), it follows that \(x>0\) and then

This yields \(\left( \frac{3}{2}-2e^{-z}+\frac{1}{2}e^{-2z}\right) ^{\alpha }\in \Sigma _{\alpha \pi /2}\).

If \(z\in \Sigma _{\theta }\setminus \Sigma _{\pi /2}\) and \(|\Im z|\le \pi \), then \(x\tan \theta \le |y|\le \pi \) and \(x\in [\pi /\tan \theta ,0]\). It yields

Set \(g(x,y)=2-e^{-x}\cos y\). We find that \(\partial g/\partial y\ge 0\) for \(y\in [x\tan \theta ,\pi ]\) and \(\partial g/\partial y\le 0\) for \(y\in [-\pi ,-x\tan \theta ]\). This leads to \(g(x,y)\ge g(x,x\tan \theta )\). Then taking the derivative of \(g(x,x\tan \theta )\) with respect to x, we get

for \(\theta \in (\pi /2, {\tilde{\theta }})\), where \({\tilde{\theta }}\) is implicitly determined by \({\tilde{g}}({\tilde{\theta }})=0\) as \(\mathrm {d}{\tilde{g}}/\mathrm {d}\theta <0\) for \(\theta \in (\pi /2,\pi )\).

Next it suffices to consider the case \(\theta \in (\pi /2,{\tilde{\theta }})\) and \(\Re (\frac{3}{2}-2e^{-z}+\frac{1}{2}e^{-2z})<0\). Let

For \(y\in [0, \pi ]\), we obtain

in view of

Together with \(f(x,y)=f(x,-y)\), it holds that

for \(0\le y\le \pi \). Similarly, from the relation \(\frac{\partial {\tilde{f}}(y,\theta )}{\partial \theta } =\frac{\partial f}{\partial x}(y/\tan \theta ,y)\frac{\mathrm {d}(y/\tan \theta )}{\mathrm {d}\theta }<0\) together with \({\tilde{f}}(0,\theta )=-\tan \theta >0\) and \(\Re (3/2-2e^{-z}+1/2e^{-2z})>0\) for \(y=\pi \), it follows that \({\tilde{f}}(y,\theta )>0\) for any \(\theta \in (\pi /2,{\tilde{\theta }})\). Then taking the derivative of \({\tilde{f}}\) with respect to y, we deduce that

where \(y_{\theta }\) satisfies \(\frac{\partial {\tilde{f}}}{\partial y}(y_{\theta },\theta )=0\). Thus there corresponds some \(\vartheta \in (\pi /2,\pi )\), defined by \({\tilde{f}}(y_{\theta },\theta )=-\tan (\vartheta )\) such that \((\frac{3}{2}-2e^{-z}+\frac{1}{2}e^{-2z})\in \Sigma _{\vartheta }\). This completes the proof. \(\square \)

Reprints and permissions