A fast linearized finite difference method for the nonlinear multi-term time-fractional wave equation

TL;DR

This paper introduces a fast, linearized finite difference method for solving nonlinear multi-term time-fractional wave equations, achieving second-order convergence and computational efficiency.

Contribution

It develops a novel discretization for multi-term Caputo derivatives and constructs a fully linearized scheme that simplifies solving nonlinear problems.

Findings

Second-order convergence in discrete H1-norm.

Efficient solution of nonlinear multi-term fractional wave equations.

Numerical results confirm the method's accuracy and efficiency.

Abstract

In this paper, we study a fast and linearized finite difference method to solve the nonlinear time-fractional wave equation with multi fractional orders. We first propose a discretization to the multi-term Caputo derivative based on the recently established fast L2-1{\sigma} formula and a weighted approach. Then we apply the discretization to construct a fully fast linearized discrete scheme for the nonlinear problem under consideration. The nonlinear term, which just fulfills the Lipschitz condition, will be evaluated on the previous time level. Therefore only linear systems are needed to be solved for obtaining numerical solutions. The proposed scheme is shown to have second-order unconditional convergence with respect to the discrete H1-norm. Numerical examples are provided to justify the efficiency.

Equations

$$\sum _{r=0}^m{\lambda }_r{{}_0^{C}D}_t^{{\alpha }_r}u=\frac{{\partial }^{2}u}{\partial x^2}+f(u)+p(x,t),x\in \Omega ,t\in (0,T],$$

$$u(x,t)=0,x\in \partial \Omega ,t\in (0,T],$$

$$u(x,0)=\phi (x),u_t(x,0)=\psi (x),x\in \bar{\Omega },$$

$${{}_0^{C}D}_t^{\delta }u(t)=\frac{1}{\Gamma (n-\delta )}{\int }_0^t\frac{u^{(n)}(s)}{(t-s{)}^{\delta +1-n}}ds,n=\lceil \delta \rceil .$$

$$|f({\varphi }_1)-f({\varphi }_2)|\le L|{\varphi }_1-{\varphi }_2|,\forall {\varphi }_1,{\varphi }_2\in {\Omega }_f,$$

$${\delta }_t{u}^{n+\frac{1}{2}}=\frac{u^{n+1}-u^n}{\tau },{\delta }_{\hat{t}}{u}^n=\frac{u^{n+1}-u^{n-1}}{2\tau }.$$

$${}_{\lambda }{D}_t^{\alpha }u(t)=\sum _{r=0}^m{\lambda }_r{{}_0^{C}D}_t^{{\alpha }_r}u(t).$$

$${{}_{\lambda }D}_t^{\alpha }u(t)={{}_{\lambda }D}_t^{\beta }v(t).$$

$$a_0^{(\beta )}={\sigma }^{1-\beta },a_l^{(\beta )}=(l+\sigma {)}^{1-\beta }-(l-1+\sigma {)}^{1-\beta },{\mu }^{(\beta )}=\frac{{\tau }^{-\beta }}{\Gamma (2-\beta )},$$

$$b_l^{(\beta )}=\frac{1}{2-\beta }[(l+\sigma {)}^{2-\beta }-(l-1+\sigma {)}^{2-\beta }]-\frac{1}{2}[(l+\sigma {)}^{1-\beta }+(l-1+\sigma {)}^{1-\beta }],l\ge 1,$$

$$g_{k}^{(n+1,\beta)}=\mu^{(\beta)}\cdot\left\{\begin{array}[]{ll}a_{n}^{(\beta)}-b_{n}^{(\beta)},&k=0,\\ a_{n-k}^{(\beta)}+b_{n-k+1}^{(\beta)}-b_{n-k}^{(\beta)},&1\leq k\leq n-1,\\ a_{0}^{(\beta)}+b_{1}^{(\beta)},&k=n,\end{array}\right.\quad n\geq 1.$$

$${{}^{\mathcal{H}}\hat{\mathcal{D}}}_t^{\beta }v(t_{n+\sigma })\triangleq$$

$$\displaystyle+\int_{t_{n}}^{t_{n+\sigma}}(\Pi_{1,n}v(s))^{\prime}(t_{n-1+\sigma}-s)^{-\beta_{r}}ds\bigg{]}$$

$$=$$

$${\hat{g}}_k^{(n+1)}=\sum _{r=0}^m{\lambda }_r{g}_k^{(n+1,{\beta }_r)}.$$

$$\left|{{}_{\lambda }D}_t^{\beta }v(t_{n+\sigma })-{{}^{\mathcal{H}}\hat{\mathcal{D}}}_t^{\beta }v(t_{n+\sigma })\right|\le C\sum _{r=0}^m{\lambda }_r{\tau }^{3-{\beta }_r},$$

$$F(\sigma)=\displaystyle\sum_{r=0}^{m}\frac{\lambda_{r}\sigma^{1-\beta_{r}}}{\Gamma(3-\beta_{r})}\Big{[}\sigma-(s-\displaystyle\frac{\beta_{r}}{2})\Big{]}\tau^{2-\beta_{r}},$$

$$\left|t^{-\beta }-\sum _{j=1}^{N^{(\beta )}}{\omega }_j^{(\beta )}{e}^{-s_j^{(\beta )}t}\right|\le ϵ,\forall t\in [\hat{\tau },T],$$

$$N^{(\beta )}=\mathcal{O}\left(\log \frac{1}{ϵ}\left(\log \log \frac{1}{ϵ}+\log \frac{T}{\hat{\tau }}\right)+\log \frac{1}{\hat{\tau }}\left(\log \log \frac{1}{ϵ}+\log \frac{1}{\hat{\tau }}\right)\right).$$

$${}_0^C{D}_t^{\beta }v(t_{n+\sigma })$$

$$\approx \sum _{j=1}^{N^{(\beta )}}{\hat{\omega }}_j{V}_j^{(n,\beta )}+{\mu }^{(\beta )}{a}_0^{(\beta )}(v^{n+1}-v^n)\triangleq {}^{\mathcal{FH}}{\mathcal{D}}_t^{\beta }{v}^{n+\sigma },$$

$${\hat{\omega }}_j^{(\beta )}=\frac{{\omega }_j^{(\beta )}}{\Gamma (1-\beta )},V_j^{(n,\beta )}=e^{-s_j^{(\beta )}\tau }{V}_j^{(n-1,\beta )}+A_j^{(\beta )}(v^n-v^{n-1})+B_j^{(\beta )}(v^{n+1}-v^n),$$

$${}^{\mathcal{FH}}{\mathcal{D}}_t^{\beta }{v}^{n+\sigma }=\sum _{k=0}^n{}^{\mathcal{F}}{g}_k^{(n+1,\beta )}(v^{k+1}-v^k),$$

$$^{\cal F}g_{k}^{(n+1,\beta)}=\left\{\begin{array}[]{ll}\displaystyle\sum_{j=1}^{N^{(\beta)}}{\widehat{\omega}_{j}^{(\beta)}}e^{-(n-1)s_{j}^{(\beta)}\tau}A_{j}^{(\beta)},&k=0;\\ \displaystyle\sum_{j=1}^{N^{(\beta)}}{\widehat{\omega}_{j}^{(\beta)}}\Big{[}e^{-(n-k-1)s_{j}^{(\beta)}\tau}A_{j}^{(\beta)}+e^{-(n-k)s_{j}^{(\beta)}\tau}B_{j}^{(\beta)}\Big{]},{}{}&1\leq k\leq n-1;\\ \displaystyle\sum_{j=1}^{N^{(\beta)}}{\widehat{\omega}_{j}^{(\beta)}}B_{j}^{(\beta)}+\mu^{(\beta)}a_{0}^{(\beta)},&k=n;\end{array}\right.\quad n\geq 1.$$

$${}^{\mathcal{F}}{\hat{g}}_k^{(n+1)}=\sum _{r=0}^m{\lambda }_r{}^{\mathcal{F}}{g}_k^{(n+1,{\beta }_r)},$$

$${}^{\mathcal{FH}}{\hat{\mathcal{D}}}_t^{\beta }{v}^{n+\sigma }=\sum _{k=0}^n{}^{\mathcal{F}}{\hat{g}}_k^{n+1}(v^{k+1}-v^k),0\le n\le N-1.$$

$${{}_{\lambda }D}_t^{\beta }v(t_{n+\sigma })={{}^{\mathcal{FH}}\hat{\mathcal{D}}}_t^{\beta }{v}^{n+\sigma }+\mathcal{O}({\tau }^{3-{\beta }_0}+ϵ),1\le n\le N,$$

$${{}_{\lambda }D}_t^{\beta }v(t_{\sigma })={{}^{\mathcal{FH}}\hat{\mathcal{D}}}_t^{\beta }{v}^{\sigma }+\mathcal{O}({\tau }^{3-{\beta }_0}).$$

$${{}_{\lambda }D}_t^{\beta }v(t_{n+\sigma })=$$

$$=$$

$$\displaystyle-\sum\limits_{j=1}^{N^{(\beta_{r})}}\omega_{j}^{(\beta_{r})}e^{-s_{j}^{(\beta_{r})}(t_{n+\sigma}-s)}\bigg{]}ds\Bigg{\}}+{\cal O}(\tau^{3-\beta_{0}})$$

$$=$$

$$f(t_n)=(1-\sigma )f(t_{n+\sigma })+\sigma f(t_{n-1+\sigma })+\mathcal{O}({\tau }^2).$$

$${}^{\mathcal{F}}{\hat{g}}_k^{(n+1)}=$$

$${}^{\mathcal{F}}{\hat{g}}_1^{(n+1)}=$$

Full text

Key words: nonlinear fractional equation; multi-term derivative; fast algorithm; linearized method; second-order scheme

1 Introduction

Fractional derivatives were found to be more accurate tools to describe diverse materials and processes which present memory and hereditary properties, thus it gives rise to great interest in the study of fractional differential equations (FDEs). However, the closed-form solution of most of FDEs are hardly to be obtained. Investigating the efficient numerical methods for FDEs becomes a popular and also urgent topic.

Recently, the multi-term time-fractional diffusion and wave equations, which is the single fractional order in the FDEs generalized to multi-orders [5, 23, 25], were successfully applied to model various types of visco-elastic damping [27] and to describe the phenomenon of subdiffusion of oxygen in both transverse and longitudinal directions [28]. Jin et al. studied the multi-term fractional diffusion equation by using the Galerkin finite element method where its extraordinary capability of modeling anomalous diffusion phenomena in highly heterogeneous aquifers and complex viscoelastic materials were mentioned [7]. The theoretical and numerical methods to the solution of linear multi-term time-fractional equations were studied by many researchers [2, 3, 14, 15]. But the studies of numerical solutions to the nonlinear multi-term time-fractional wave equations are still limited.

In this paper, we consider an efficient numerical method to the nonlinear multi-term time-fractional wave equation:

[TABLE]

where $\Omega=(x_{L},x_{R})$, $1<\alpha_{m}<\cdots<\alpha_{0}<2$, and the positive weights $\lambda_{r}$ fulfills $\displaystyle\sum_{r=0}^{m}\lambda_{r}\leq C$. Here and hereafter, $C$ always denotes a generic positive constant. Moreover, ${{}_{0}^{C}D}_{t}^{\delta}$ denotes the Caputo derivative of order $\delta$:

[TABLE]

The nonlinear function $f$ in (1.1) satisfies the Lipschitz condition:

[TABLE]

where $\Omega_{f}$ is a suitable domain, and $L$ is a positive constant only depends on the domain $\Omega_{f}.$

The nonlocal dependence of fractional derivative is inherited by its discretizations. Thus it should be necessary to develop high accurate and/or fast algorithms to FDEs in order to save the computational costs and the memory storage. Moreover, the problem under consideration is a nonlinear time-fractional equation, and classical nonlocal numerical schemes blending with iterative methods would take more expensive computation and have more complexity in the analysis. So it is really competitive and also necessary to construct linearized numerical methods for nonlinear time-fractional problems [12, 16, 17, 18, 33]. It is worth to mention that the solution of many time-fractional differential equations typically displays a weak singularity near the initial time [8, 9, 10, 19, 20, 21, 26, 29, 30], which leads to the loss of time accuracy for many related high-order numerical methods. The nonuniform grids technique [10, 29] is a very popular method to recover the full accuracy very recently. Since our proposed linearized method for the considered nonlinear time-fractional problem is based on a weighted approach, the analysis will be very difficult if the grids are nonuniform. In view of the facts, based on the fast ${\cal L}2$-$1_{\sigma}$ formula (named ${\cal FL}2$-$1_{\sigma}$) investigated in [34], we construct a fast linearized finite difference scheme on uniform grids to solve the nonlinear multi-term time fractional wave equation. This may offset the possible accuracy loss in the computational sense.

Traditional direct numerical methods for time-fractional PDEs require ${\cal O}(MN)$ memory and ${\cal O}(MN^{2})$ work, where $M$ and $N$ denote the total number of space steps and time steps, respectively. Based on an efficient sum-of-exponentials (SOE) approximation for the kernel function $t^{-\beta-1}$ on the interval $[\tau,T]$, where $\beta\in(0,1)$ and $\tau$ is the time step size, and combined with the classical ${\cal L}1$ discretization [13, 22, 32], Jiang et al. [6] proposed an efficient fast evaluation to the Caputo derivative, and then applied it to solve time-fractional diffusion equations. The corresponding fast numerical algorithm requires only ${\cal O}(MN_{\exp})$ memory and ${\cal O}(MNN_{\exp})$ work, where $N_{\exp}$ is the number of exponentials and of order ${\cal O}(\log N)$. Based on the ${\cal L}2$-$1_{\sigma}$ discretization [1], Yan et al. [34] then constructed a fast and second-order ${\cal FL}2$-$1_{\sigma}$ numerical method for time-fractional diffusion equations. The computational cost just require ${\cal O}(MN\log^{2}N)$ and the overall storage is ${\cal O}(M\log^{2}N)$.

In this paper, based on the ${\cal FL}2$-$1_{\sigma}$ formula in [34] and the weighted idea in [18], we first approximate the multi-term derivative in (1.1) at time point $t_{n}$ ($n\geq 1$) by a weighted discretization which will solve the function to time level $t_{n+1}$, and give a fitted time approximation on the diffusion term. Then we apply the proposed approximation to construct a linearized and fast finite difference scheme to solve the nonlinear equation (1.1)–(1.3), which will evaluate the nonlinear term at previous time level. Thus only linear systems are needed to be solved for obtaining approximated solutions. With some important rigorously verified properties of the discrete coefficients, we show that our proposed fast linearized numerical scheme is unconditionally convergent with the order ${\cal O}(h^{2}+\tau^{2}+\epsilon)$, where $h$ is the spatial step size and $\epsilon$ is the tolerance error in the approximation of SOE to the kernel.

The structure of the paper is as follows. In section 2, based on the ${\cal FL}2$-$1_{\sigma}$ formula and a weighted approach, we propose a discretization to approximate the multi-term Caputo derivative at time grid point $t_{n}$ ($n\geq 1$). With some necessary properties of the corresponding coefficients being verified, we give the truncation error analysis for the proposed discretization. Then applying the discretization and using a special approximation for first time level solution, we construct a fully fast-linearized finite difference scheme to solve the considered problem (1.1)–(1.3). In section 3, we first estimate some important properties of the discrete coefficients, and then we show that our proposed fast-linearized scheme is unconditionally convergent with second-oder accuracy with respect to discrete $H^{1}$-norm. In section 4, numerical examples are carried out to confirm the efficiency of the numerical scheme. A brief conclusion is followed in section 5.

2 The fast linearized numerical method

Notations for clarifying some coefficients and parameters:

• •

$g_{k}^{(n+1,\beta)}$ —– the ${\cal L}2$-$1_{\sigma}$ type coefficients;

• •

${\widehat{g}}_{k}^{(n+1)}$ —– the ${\cal L}2$-$1_{\sigma}$ type coefficients of multi-term Caputo derivative;

• •

$N^{(\beta)},s_{j}^{(\beta)},\omega_{j}^{(\beta)},{\widehat{\omega}}_{j}^{(\beta)}$ —– the corresponding parameters of SOE approximation;

• •

${}^{\cal F}g_{k}^{(n+1,\beta)}$ —– the ${\cal FL}2$-$1_{\sigma}$ type coefficients;

• •

${}^{\cal F}{\widehat{g}}_{k}^{(n+1)}$ —– the ${\cal FL}2$-$1_{\sigma}$ type coefficients of multi-term Caputo derivative;

• •

${\bf g}_{k}^{(n+1)}$ —– the refined ${\cal FL}2$-$1_{\sigma}$ type coefficients of multi-term Caputo derivative.

2.1 Preliminary

We first introduce some temporal notations. For a given positive integer $N$, denote the uniform time step size by $\tau=\frac{T}{N}$, and take $t_{n}=n\tau~{}(0\leq n\leq N),~{}t_{n+\sigma}~{}(\sigma\in[0,1])$. Let {\cal W}_{\tau}=\big{\{}w^{n}|0\leq n\leq N\big{\}}. For any $u^{n}\in{\cal W_{\tau}},$ denote

[TABLE]

For simplicity, we denote the multi-term derivative by

[TABLE]

In the rest of the paper, we take $v=u_{t}$, $\alpha_{r}=\beta_{r}+1$ or $\alpha=\beta+1$ when it is required. Then

[TABLE]

Our approximation to fractional derivative is based on the ${\cal L}2$-$1_{\sigma}$ discretization in [1, Lemma 2], we first introduce the corresponding coefficients. For $\beta\in(0,1)$, denote

[TABLE]

where the particular $\sigma$ (here and hereafter) is generated in Lemma 2.1, which may be different from the $\sigma$ in [1]. Further denote $g_{0}^{(1,\beta)}=\mu^{(\beta)}a_{0}^{(\beta)}$ and

[TABLE]

The high-order approximation proposed in [4] to multi-term Caputo derivative is:

[TABLE]

where $\Pi_{2,k}$ and $\Pi_{1,n}$ are the quadratic and linear polynomials [4], respectively; and

[TABLE]

The next lemma is given with respect to the truncation error between ${{}_{\lambda}D}_{t}^{\beta}v(t_{n+\sigma})$ and ${\cal{}^{H}{\widehat{D}}}_{t}^{\beta}v(t_{n+\sigma})$, where the basic way to approximate the multi-term derivative is shown.

Lemma 2.1**.**

([4])* Suppose $v(t)\in{\cal C}^{3}[0,t_{n+1}]$, it holds that*

[TABLE]

where \sigma\in\big{[}1-\frac{\beta_{0}}{2},1-\frac{\beta_{m}}{2}\big{]} is the root of equation

[TABLE]

generated by the method of Newton iteration.

To construct fast numerical scheme, we introduce a ${\cal FL}2$-$1_{\sigma}$ formula, proposed in [34], which is based on the ${\cal L}2$-$1_{\sigma}$ formula and a SOE approximation to the kernel function $t^{-\beta}$ in the Caputo derivative. The SOE approximation reads as:

Lemma 2.2**.**

([34])* For any $\beta\in(0,1)$, tolerance error $\epsilon$, cut-off time step size ${\widehat{\tau}}$ and final time $T$, there are some positive integer $N^{(\beta)}$, positive points $s_{j}^{(\beta)}$ and corresponding positive weights $\omega_{j}^{(\beta)},~{}(j=1,2,\ldots,N^{(\beta)})$ satisfying*

[TABLE]

and the number of exponentials needed is of the order

[TABLE]

Taking ${\widehat{\tau}}\leq 10^{-1}$, then we have N^{(\beta)}={\cal O}\Big{(}\log\displaystyle\frac{1}{\widehat{\tau}}\log\displaystyle\frac{T}{\widehat{\tau}}+\log\displaystyle\frac{1}{\widehat{\tau}}\log\displaystyle\frac{1}{\widehat{\tau}}\Big{)}. Then let ${\widehat{\tau}}=\sigma\tau$ and $\tau<\min\{\sigma,\max\{\sigma^{2},T\sigma^{2}\}\}$, we have N^{(\beta)}={\cal O}\Big{(}\log\displaystyle\frac{1}{\tau}\log\displaystyle\frac{T}{\tau}+\log\displaystyle\frac{1}{\tau}\log\displaystyle\frac{1}{\tau}\Big{)}.

For a single-term Caputo derivative of order $\beta\in(0,1)$ at time $t_{n+\sigma}$: ${{}_{0}^{C}}D_{t}^{\beta}v(t_{n+\sigma})$, the ${\cal FL}2$-$1_{\sigma}$ type formula here is estimating ${(t_{n+\sigma}-s)^{-\beta}}$ $(s\in(0,t_{n}))$ by Lemma 2.2 and $v^{\prime}(s)$ $(s\in(t_{n},t_{n+\sigma}))$ by linear polynomial, and then estimating the history part of the integral by using a recursive relation and quadratic interpolation, that is

[TABLE]

where

[TABLE]

with $A_{j}^{(\beta)}=\int_{0}^{1}(\frac{3}{2}-s)e^{-s_{j}^{(\beta)}\tau(\sigma+1-s)}ds$, $B_{j}^{(\beta)}=\int_{0}^{1}(s-\frac{1}{2})e^{-s_{j}^{(\beta)}\tau(\sigma+1-s)}ds$, and $V_{j}^{0}=0$.

In fact, the ${\cal FL}2$-$1_{\sigma}$ type discretization can be rewritten as (see [34])

[TABLE]

where ${}^{\cal F}g_{0}^{(1,\beta)}=\mu^{(\beta)}a_{0}^{(\beta)}$ and

[TABLE]

Denote ${{}^{\cal FH}}{\cal\widehat{D}}_{t}^{\beta}v^{n+\sigma}=\displaystyle\sum_{r=0}^{m}\lambda_{r}{{}^{\cal FH}}{\cal D}_{t}^{\beta_{r}}v^{n+\sigma}$ and

[TABLE]

then the multi-term ${\cal FL}2$-$1_{\sigma}$ type discritization can be rewritten as

[TABLE]

We have the following lemma which presents the error of the above fast approximation to multi-term Caputo derivative.

Lemma 2.3**.**

For $v(t)\in{\cal C}^{3}[0,t_{n}]$, it holds that

[TABLE]

Proof.

For $1\leq n\leq N$, it follows from Lemmas 2.1 and 2.2 that

[TABLE]

where $\epsilon_{r}$ is the tolerance related to the SOE approximation of the kernel $(t_{n-1+\sigma}-s)^{-\beta_{r}}$, and $\epsilon=\max\limits_{0\leq r\leq m}\{\epsilon_{r}\}$. As ${{}^{\cal FH}{\cal{\widehat{D}}}}_{t}^{\beta}v^{\sigma}$ $=$ ${{}^{\cal H}{\cal{\widehat{D}}}}_{t}^{\beta}v^{\sigma}$, the second conclusion can be directly obtained from Lemma 2.1. ∎

2.2 Fast approximation to multi-term derivative

2.2.1 A weighted method to approximate ${{}_{\lambda}D}_{t}^{\alpha}u(t_{n})$

Now we go to study a second-order approximation to the multi-term derivative ${{}_{\lambda}D}_{t}^{\alpha}u(t_{n})$ based on the method of order reduction and a weighted approach. The following lemma provides an important weighted approach for discretizing ${{}_{\lambda}D}_{t}^{\alpha}u(t_{n})~{}(1\leq n\leq N-1)$ that leads to the linearized approximation to our considered nonlinear problem.

Lemma 2.4**.**

([18])* For any $f(t)\in{\cal C}^{2}[t_{n-1+\sigma},t_{n+\sigma}]$, it holds that*

[TABLE]

Before applying Lemma 2.4, we first observe that (2.3) and (2.4) can yield ($n\geq 1$)

[TABLE]

where $b_{n}=\displaystyle\sum_{r=0}^{m}\lambda_{r}\displaystyle\sum_{j=1}^{N^{(\beta_{r})}}{\widehat{\omega}_{j}^{(\beta_{r})}}e^{-(n-1)s_{j}^{(\beta_{r})}\tau}B_{j}^{(\beta_{r})}$. Then, for $n\geq 1$,

[TABLE]

For simplicity of representation, we take

[TABLE]

Invoking the weights investigated in Lemma 2.4, we have

[TABLE]

However, we find that the last coefficient ${{}^{\cal F}{\widehat{g}}_{0}^{(n+1)}}-\displaystyle\frac{\sigma}{1-\sigma}b_{n}$ is not always positive, which may cause difficulty in our analysis. In view of this, we regroup the last term

[TABLE]

where

[TABLE]

Thus

[TABLE]

Note that the function $v$ is still not discretized by $u$ in the above derivations, to obtain a fully discretization for ${{}_{\lambda}D}_{t}^{\alpha}u(t_{n})~{}(1\leq n\leq N-1)$, we next first verify some necessary properties of the coefficients.

Lemma 2.5**.**

The coefficients in (2.9) and $b_{n}$ satisfy $(n\geq 1)$

[TABLE]

where $\alpha_{i_{1}},\alpha_{i_{2}}\in(\alpha_{m},\alpha_{0})$.

Proof.

$\bullet$ Estimation of $(a)$: Since

[TABLE]

then

[TABLE]

Hence, there exist $\alpha_{i_{0}},\alpha_{i_{1}}\in(\alpha_{m},\alpha_{0})$ such that

[TABLE]

$\bullet$ Estimation of $(b)$: Note that

[TABLE]

that is $B_{j}^{(\beta_{r})}<2\left(A_{j}^{(\beta_{r})}-\frac{1}{2}B_{j}^{(\beta_{r})}\right)$, then it follows from (2.11) that

[TABLE]

$\bullet$ Estimation of $(c)$: We first notice that

[TABLE]

Then, for $1\leq k\leq n-1$ we have

[TABLE]

and for $k=n$

[TABLE]

Therefore

[TABLE]

Thus $(c)$ can be verified by the fact

[TABLE]

$\bullet$ Estimation of $(d)$: As $A_{j}^{(\beta_{r})}>0$, for $n\geq 2$ we have

[TABLE]

for $n=1$, noticing (2.12), we have

[TABLE]

Thus by $(b)$ and (2.14)–(2.15), we get the verification of $(d)$:

[TABLE]

∎

Taking

[TABLE]

From the subsection 2.2 in [18], we know that

[TABLE]

where the truncation errors satisfy

[TABLE]

provided $u(t)\in{\cal C}^{3}[0,T]$. Moreover,

[TABLE]

provided $u(t)\in{\cal C}^{4}[0,T]$.

We are now ready to show a fully discretization to ${{}_{\lambda}D}_{t}^{\alpha}u(t_{n})~{}(1\leq n\leq N-1)$ and its truncation error.

Lemma 2.6**.**

Suppose $u(t)\in{\cal C}^{4}[0,T]$. Denote

[TABLE]

Then

[TABLE]

Proof.

By (2.10) and (2.19)–(2.20), we can derive

[TABLE]

where

[TABLE]

With the help of Lemma 2.5 (b), (c) and (d), we have

[TABLE]

in which ${\bf g}_{0}^{(n+1)}{\tilde{b}^{n}}=\frac{3\sigma-1}{2(1-\sigma)}b_{n}\leq C{\bf g}_{0}^{(n+1)}$ has been used.

Applying the relation (2.1), Lemmas 2.3 and 2.4, we have

[TABLE]

Therefore,

[TABLE]

∎

2.2.2 The first time level approximation

Note that the discretization (2.20) will be used to solve the grid function $u^{n+1}~{}(n\geq 1)$. For the first level grid function, denote ${\bf g}_{0}^{(1)}=\frac{{{}^{\cal F}{{\widehat{g}}_{0}^{(1)}}}}{1-\sigma}$ and

[TABLE]

Then it can be observed that

[TABLE]

By the relation (2.1) and Lemma 2.3 , we further get

[TABLE]

where

[TABLE]

2.3 The fast linearized scheme

Before proposing our numerical scheme, we need an important lemma which provides a fitted approximation (as it plays a significant role in our later convergence analysis) to the time discretization on diffusion term. The lemma reads as:

Lemma 2.7**.**

([18])* Suppose $u(t)\in{\cal C}^{2}[0,T]$. For $n\geq 1$, it holds that*

[TABLE]

where

[TABLE]

For a given positive integer $M$, denote $h=\frac{b-a}{M}$ be the spatial step size, and ${\bar{I}}=\{i|~{}i=0,1,\ldots,M\}$ and $I=\{i|~{}i=1,2,\ldots,M-1\}$ be the index spaces, and take $x_{i}=ih~{}(i\in{\bar{I}})$ be the spatial uniform partition. Denote $\varphi_{i}=\varphi(x_{i})$ and $\psi_{i}=\psi(x_{i})$. Suppose {\cal V}_{h}=\big{\{}u|u=\{u_{i}|0\leq i\leq M\},~{}u_{0}=u_{M}=0\big{\}}, then for any $u,~{}v\in{\cal V}_{h}$, the discrete operators are needed:

[TABLE]

Denote the numerical solution at the point $(x_{i},t_{n})$ by $u_{i}^{n}$, and $p(x_{i},t_{n})$ by $p_{i}^{n}$.

Considering the equation (1.1)–(1.3) on the grid point $(x_{i},t_{n})$:

[TABLE]

Applying Lemmas 2.6 and 2.7, and standard approximation on space derivative, the above equation yields

[TABLE]

where

[TABLE]

The equation (2.23) solves the solution $u_{i}^{n+1}$, thus it is a linearized equation since the nonlinear term is stayed on the previous time level $t_{n}$.

For the first level solution, considering the equation (1.1)–(1.3) on the grid point $(x_{i},t_{\sigma})$:

[TABLE]

With the help of Taylor expansion, we use the initial data to approximate $u(x_{i},t_{\sigma})$:

[TABLE]

Then for $\Omega_{f}\supset[\inf\{\varphi+\sigma\tau\psi\},\sup\{\varphi+\sigma\tau\psi\}]$, it follows from (2.24), (2.21) and (2.25) that

[TABLE]

where ${\tilde{R}}_{i}^{\sigma}={\cal O}(\tau^{3-\beta_{0}}+\tau^{2})={\cal O}(\tau^{2})$. We can see that equation (2.26) is also a linearized approximation.

Omitting $R_{i}^{n+1}$ in (2.23), and ${\tilde{R}}_{i}^{\sigma},R_{i}^{1}$ in (2.26), we obtain the following fully fast linearized scheme for solving the nonlinear problem (1.1)–(1.3):

[TABLE]

Remark 2.8**.**

In order to interpret the fast algorithm of the proposed scheme (2.27)–(2.30), we display the efficient computation for the fractional discretization ${{}^{\cal FH}}{\cal{\widehat{D}}}_{t}^{\alpha}u_{i}^{n+1}$ based on the recursive relation (2.2) and the refined grid functions (2.16)–(2.18):

For $n=0,$ we can obtain $u_{i}^{1}$, $v_{i}^{1}$ (combining (2.5)) by solving (2.28).

For $n=1,$ we can solve the following system to obtain $u_{i}^{2},~{}v_{i}^{2}$:

[TABLE]

where ${\bf a_{0}}=\displaystyle\sum_{r=0}^{m}\lambda_{r}\mu^{(\beta_{r})}a_{0}^{(\beta_{r})},~{}V_{j}^{1}=A_{j}^{(\beta_{r})}(v_{i}^{1}-v_{i}^{0})+B_{j}^{(\beta_{r})}(v_{i}^{2}-v_{i}^{1}),$ and $v_{i}^{2-\sigma},~{}v_{i}^{1-\sigma}$ are defined by (2.8).

For $2\leq n\leq N-1$,

[TABLE]

where a recursive relation

[TABLE]

with V_{j}^{1-\sigma}=(1-\sigma)\Big{[}A_{j}^{(\beta_{r})}(v_{i}^{1}-v_{i}^{0})+B_{j}^{(\beta_{r})}(v_{i}^{2}-v_{i}^{1})\Big{]}.

3 Analysis of the proposed scheme

3.1 Some necessary lemmas

To show the convergence of proposed scheme, we need some important lemmas.

Lemma 3.1**.**

([34])* For the sequence $\{{{}^{\cal F}{g_{k}^{(n+1,\beta)}}}\}~{}(k=0,\ldots,n)$, it holds that*

[TABLE]

and

[TABLE]

where $C^{\cal F}=\min\{\mu^{(\beta)}a_{0}^{(\beta)},{{}^{\cal F}}g_{0}^{(n+1,\beta)}\}.$

Lemma 3.2**.**

([4])* There exists a number $\tau_{0}>0$, when $\tau\leq\tau_{0}$, it holds that*

[TABLE]

where

[TABLE]

Lemma 3.3**.**

There exists a number $\tau_{0}>0$, when $\tau\leq\tau_{0}$, it holds that ($n\geq 1$)

[TABLE]

for a sufficiently small $\epsilon$.

Proof.

From the relation (2.13) and Lemma 3.1, we easily get

[TABLE]

Notice that

[TABLE]

with property $(a)$ in Lemma 2.5, (3.1) is obtained.

For $n\geq 2$, then by Lemma 3.1, we have

[TABLE]

By Lemma 3.2, we reach

[TABLE]

Combining the above inequality with (3.1), we can conclude that

[TABLE]

holds as long as

[TABLE]

For $n=1$, using the similar way as $n\geq 2$, we have

[TABLE]

Next we show $b_{1}>0$. Note that

[TABLE]

For the integral term

[TABLE]

Recursively, the above integral has the lower bound

[TABLE]

It is easy to check that

[TABLE]

hence $b_{1}>0$.

Therefore, from (3.4), it holds that $(2\sigma-1){\bf g}_{1}^{(2)}-\sigma{\bf g}_{0}^{(2)}>0$ when

[TABLE]

Thus in general, (3.2) holds if

[TABLE]

∎

In the following parts, we always suppose conditions in Lemma 3.3 are fulfilled when required.

Lemma 3.4**.**

([1])* If the positive sequence $\{d_{k}^{(n+1)}|~{}0\leq k\leq n,n\geq 1\}$ is strictly decrease for $k$, and satisfies $(2\sigma-1)d_{n}^{(n+1)}-\sigma d_{n-1}^{(n+1)}>0$ for a constant $\sigma\in(0,1)$ . Then*

[TABLE]

We now present a particular form of Lemma 3.4, which will be used in our analysis.

Lemma 3.5**.**

For any real sequence $F^{n}$, the following inequality holds:

[TABLE]

Proof.

By Lemma 3.3, incorporating $d_{k}^{(n+1)}={\bf g}_{k}^{(n+1)}$, $y^{n}={\widehat{v}}^{n-\sigma}$ $(n\geq 1)$ and $y^{0}={\widehat{v}}^{-\sigma}+\frac{1}{{\bf g}_{0}^{(n+1)}}F^{n}$ with Lemma 3.4, and noticing the properties (2.6)–(2.7), we can get

[TABLE]

∎

Lemma 3.6**.**

([24])(Gronwall’s inequality)* Let $\{G_{{j}}\}$ and $\{k_{{j}}\}$ be nonnegative sequences satisfying*

[TABLE]

where $K\geq 0$. Then

[TABLE]

Another important lemma concerns some properties about the coefficients for analysis:

Lemma 3.7**.**

The coefficients ${\bf g}_{0}^{(k+1)}~{}(0\leq k\leq n)$ and ${\bf g}_{k}^{(n+1)}$ satisfy

[TABLE]

Proof.

Since ${\bf g}_{0}^{(1)}=\frac{{}^{\cal F}{\widehat{g}}_{0}^{(1)}}{1-\sigma}$ and

[TABLE]

then $(a)$ can be easily obtained.

By Lemma 2.5 $(a)$, we have

[TABLE]

uniting $\tau{\bf g}_{0}^{(1)}\leq C\tau^{2-\alpha_{0}}$, $(b)$ is verified .

Applying Lemma 2.5 $(a)$,

[TABLE]

so $(c)$ can be obtained by combining $(a)$.

From Lemma 3.3, we know that ${\bf g}_{k}^{(n+1)}>{\bf g}_{0}^{(n+1)}$ for $1\leq k\leq n$. Then applying Lemma 2.5 $(a)$ again, we get

[TABLE]

so (d) is proved. ∎

3.2 The unconditional convergence

For two mesh functions $v_{h},w_{h}\in\mathcal{V}_{h}$, we define the inner product and norms (the discrete $L_{2}$-norm and a semi-norm)

[TABLE]

Furthermore, we introduce the discrete $H^{1}$-norm $\|v\|_{H^{1}}=\sqrt{\|v\|^{2}+|v|_{1}^{2}}.$

Referring to [11, 31], we know that $\|v\|\leq(x_{R}-x_{L})/{\sqrt{6}}|v|_{1}$, then

[TABLE]

Now we denote the errors $e_{i}^{n}=u(x_{i},t_{n})-u_{i}^{n}$, $i\in{\bar{I}}$ and $0\leq n\leq N$, and denote

[TABLE]

in which

[TABLE]

Take

[TABLE]

From subsection 2.3, we easily obtain the following error system:

[TABLE]

Next theorem will show that the error $e_{i}^{n}$ is bounded in the $H^{1}$-norm unconditionally.

Theorem 3.8**.**

Let $u(x,t)$ be the solution of the problem (1.1)–(1.3). Assume $u(x,t)\in{\cal C}^{4}(\Omega)\cap{\cal C}^{4}[0,T]$. Let $\{u_{i}^{n},i\in{\bar{I}},0\leq n\leq N\}$ be the solutions of the scheme (2.27)–(2.30). Then the errors $e_{i}^{n}$ satisfy

[TABLE]

Proof.

We finish the proof by using some techniques in Theorem 3.5 in [33].

From (3.11), one has $\|e^{0}\|_{\infty}=0$ . We first utilize mathematical induction to show

[TABLE]

where

[TABLE]

with $C_{\alpha}=Ct_{n+\sigma}^{2-\alpha_{i_{1}}}$.

It follows from (3.5)–(3.6) that

[TABLE]

Then by (2.22) and Lemma 3.7 $(a)$, we have

[TABLE]

Hence (3.12) holds for $n=0$ and $n=1$ ($\tau\leq 1$).

Suppose (3.12) is valid for $1\leq n\leq q$ ($1\leq q\leq N-1$), that is

[TABLE]

Before proving that (3.12) is valid for $n=q+1$, we show that the numerical solutions $u^{n}~{}(1\leq n\leq q)$ are uniformly bounded based on the inductive assumption. By using Lemma 3.6 (Gronwall’s inequality) and Lemma 3.7 $(c)$ on (3.15), we can obtain

[TABLE]

With the smooth assumption on the exact solution, which yields $\|U^{n}\|_{\infty}\leq C_{u}$ for a positive constant $C_{u}$, it follows

[TABLE]

Hence we can take $\Omega_{f}=[-C_{u}-1,C_{u}+1]\cup[\inf\{\varphi+\sigma\tau\psi\},\sup\{\varphi+\sigma\tau\psi\}]$ in the rest proof.

We now verify that (3.12) is valid for $n=q+1$.

Taking the inner product of (3.8) with

[TABLE]

we have

[TABLE]

where {\tilde{F}}_{i}^{n}=\big{[}f(u(x_{i},t_{n}))-f(u_{i}^{n})\big{]}+R_{i}^{n+1}.

With the boundary values being zero, it is easy to verify that

[TABLE]

Utilizing Lemma 3.5, we get

[TABLE]

Substituting (3.17) and (3.18) into (3.16), we obtain

[TABLE]

where

[TABLE]

Summing up (3.19) for $n$ from $1$ to $p$ yield

[TABLE]

that is

[TABLE]

It can be verified by using Cauchy-Schwarz inequality and Lemma 3.7 $(d)$ that

[TABLE]

Furthermore, the inequality $2(y^{2}+z^{2})\geq(y+z)^{2}$ gives

[TABLE]

Consequently, it follows from (3.20)–(3.22) that

[TABLE]

where

[TABLE]

We further note that the term on the left hand side of (3.23) satisfies

[TABLE]

Combining (3.23) and (3.25), we get

[TABLE]

Referring to the proof of Theorem 3.5 in [33], we discuss the following two cases:

Case (I) \|e^{q+1}\|+\sqrt{2C_{\alpha}}\big{|}\sigma e^{q+1}+(1-\sigma)e^{q}\big{|}_{1}\leq\|e^{q}\|+\sqrt{2C_{\alpha}}\big{|}\theta e^{q}+(1-\sigma)e^{q-1}\big{|}_{1}.

Similar to that in [33], we can obtain $S^{q+1}\leq S^{q}$ in this case, so (3.12) follows directly.

Case (II) \|e^{q+1}\|+\sqrt{2C_{\alpha}}\big{|}\sigma e^{q+1}+(1-\sigma)e^{q}\big{|}_{1}\geq\|e^{q}\|+\sqrt{2C_{\alpha}}\big{|}\sigma e^{q}+(1-\sigma)e^{q-1}\big{|}_{1}.

In this situation, we have

[TABLE]

With the above inequality and (3.26), it follows

[TABLE]

We next estimate $B^{q}$ term by term. Firstly, with (3.14), we have

[TABLE]

Combining (3.7) and (3.14), we get

[TABLE]

From (3.5), (3.9) and (2.22), we have

[TABLE]

Then with Lemma 3.7 $(a)$,

[TABLE]

Observing (2.3) and (2.9), we know that ${\bf g}_{0}^{(n)}={{}^{\cal F}{{\widehat{g}}_{0}^{n}}}-\frac{1}{2}b_{n-1}$ for $n\geq 2$, ${\bf g}_{0}^{(1)}={{}^{\cal F}{{\widehat{g}}_{0}^{1}}}$, and ${\bf g}_{1}^{(n+1)}={{}^{\cal F}{\widehat{g}}_{1}^{n+1}}={{}^{\cal F}{\widehat{g}}}_{0}^{n}+b_{n}$ for $n\geq 1$. Then we get

[TABLE]

where Lemma 2.5 (b) and Lemma 3.7 (b) have been used in the last two inequalities. Therefore,

[TABLE]

From (3.6) and Lemma 2.5 (b),

[TABLE]

thus it follows

[TABLE]

Note that

[TABLE]

where Lemma 3.7 $(b)$ and $(c)$ have been used.

Thus, (3.24) and (3.27)–(3.33) yield

[TABLE]

Similar discussion with Case (II) in Theorem 3.5 of [33], we can get

[TABLE]

Hence, the inequality (3.12) is clarified according to (3.34)–(3.35).

Consequently, we can apply Lemma 3.6 (Gronwall’s inequality) and Lemma 3.7 $(c)$ on (3.12) to conclude

[TABLE]

With (3.13), we then obtain the desired result.∎

Remark 3.9**.**

In a way similar to the proof of convergence, we can show that the numerical scheme (2.27)–(2.30) is unconditionally stable with respect to discrete $H^{1}$-norm. Readers can refer to [33] for more details.

4 Numerical experiments

In this section, we carry out numerical experiments for the proposed finite difference schemes (2.27)-(2.30) to illustrate our theoretical statements. The $H^{1}$-norm errors between the exact and the numerical solutions

[TABLE]

are shown in the following tables and the convergence rates defined by

[TABLE]

are also recorded.

In the following tests, scheme1 stands for the scheme (2.27)-(2.30), and scheme2 represents the direct scheme which is similar to the scheme1 but without using the SOE.

We consider the problem for $x\in[0,1],~{}T=1$ and the forcing term

[TABLE]

is chosen to such the exact solution $u(x,t)=\sin(\pi x)t^{4}$, where

[TABLE]

In Table 1, taking serval sets of $\lambda$ and $\alpha$, with three different modeling cases, the $E_{1}(h,\tau)$ and the temporal convergence of scheme (2.27)-(2.30) are presented, which confirms the second-order convergence of the difference scheme with respect to temporal direction. In Table 2, the numerical results are shown in temporal direction where the numerical results of three cases are also demonstrated and one can realize that scheme2 is of second-order convergence. For the spatial direction, in Table 3, the second-order accuracy of scheme (2.27)-(2.30) for three cases are also shown with fixed $\tau=\frac{1}{1000}.$ Above all, the chosen $\epsilon_{r}$ are $\tau^{4-\alpha_{r}}\times 10^{-3}$, $r=0,1,2.$ Table 4 demonstrates CPU time in seconds of scheme1 and scheme2. Figure 1 presents the comparison between two schemes about the memory in bytes with different temporal step $\tau=\frac{1}{10000}$ to $\frac{1}{40000}$, and actually the results of three cases are the same. One can check that the scheme1 do require less time and memory and if $\tau$ is small enough, the comparison can be more evident.

5 Conclusion

We considered a fast and linearized finite difference method for solving the nonlinear multi-term time-fractional wave equation. The proposed scheme based on the fast ${\cal L}2$-$1_{\sigma}$ discretization, the multi-term ${\cal L}2$-$1_{\sigma}$ type discretization and a weighted approach. By showing some important properties of the refined coefficients of fully discretization, we obtained the truncation error of our proposed weighted discretization to the multi-term Caputo derivative, and we displayed the unconditional convergence rigorously. The accuracy and efficiency of proposed method are well demonstrated by several numerical tests.