A new three-term conjugate gradient algorithm with modified gradient-differences for solving unconstrained optimization problems

1.   Introduction

In this paper, the following one-sided tempered fractional diffusion equations are considered:

and

where $1 < \alpha < 2$, $\lambda\geq0,$ the diffusion coefficient K is positive, and $f(x, t)$ is the source term. $D^{\alpha, \lambda}_{a, x}u(x, t)$ and $D^{\alpha, \lambda}_{x, b}u(x, t)$ represent the left and right Riemann-Liouville tempered fractional derivatives, respectively, and are defined as

The concept of fractional derivatives appeared almost at the same time as integer derivatives, but the lack of application background of fractional derivatives at the beginning of their appearance made fractional models not widely developed until recent decades. Fractional models are widely used in physics ^[1,2,3,4], finance ^[5], biology ^[6], and hydrology ^[7,8,9]. Recently, the study of fractional diffusion equations has attracted a lot of attention, in which the integral derivatives of diffusion equations are replaced by the fractional derivatives to obtain fractional diffusion equations. Usually the time fractional derivatives describe anomalous sub-diffusion, and the space fractional derivatives describe anomalous super-diffusion ^[10].

The analytical solutions of fractional differential equations usually cannot be obtained because of the nonlocality of fractional derivatives. Therefore, a lot of attention has been paid to the development of high-precision numerical methods, and a lot of results have been obtained ^{[11,12,13,14,15]}. For the space fractional advection-dispersion equations, Meerschaert and Tadjeran ^[11] point out that the standard Grünwald difference operator to approximate the Riemann-Liouville fractional derivative is unconditionally unstable regardless of the implicit and explicit Euler methods, so a modified Grünwald difference operator, called the shifted Grünwald difference operator, is proposed to solve this problem. Based on the shift Grünwald difference operator and its idea, more research has been done ^{[16,17,18,19,20]}. While developing numerical methods for fractional diffusion equations, some researchers have found that if the power law (waiting time or jump length) is tempered by an exponential factor, it has practical advantages ^[21,22,23]. Therefore, the time and space tempered fractional derivatives are obtained, and the time tempered fractional derivative yields the time tempered fractional diffusion equation ^[22]. Luo et al. ^[24] proposed an effective Lagrange-quadratic spline optimal collocation method to solve it. The space tempered fractional derivative yields the space tempered fractional diffusion equation ^[23]. It is very important to develop numerical methods for tempered fractional diffusion equations.

A common way to solve the space tempered fractional diffusion equation is to develop a modified Grünwald difference operator based on the idea of the shifted Grünwald difference operator, called the tempered shifted Grünwald difference operator by Baeumer and Meerschaert ^[25], which is used to solve the tempered fractional advection-dispersion equation, but the space direction has only first-order precision. Based on the idea of the tempered shifted Grünwald difference operator, Li and Deng ^[26] obtained a class of second-order approximations of left and right Riemann-Liouville tempered fractional derivatives by using the weighted idea, called the tempered weighted and shifted Grünwald difference operators, which combined with the Crank-Nicolson method is used to solve the two-sided tempered fractional diffusion equations. See more studies on space-time tempered fractional differential equations ^[27,28], time tempered fractional differential equations ^[29,30], and space tempered fractional differential equations ^{[31,32,33,34,35]}. Similar to the quasi-compact scheme of fractional diffusion equations ^[36,37,38], it is natural to think of applying the quasi-compact technique to the numerical solutions of two-sided tempered fractional diffusion equations. However, due to the incompatibility of the left and right Riemann-Liouville tempered fractional quasi-compact operators, Yu et al. ^[39] considered the quasi-compact technique on the one-sided tempered fractional diffusion equations and obtained that the numerical schemes are stable and convergent with order $O(\tau+h^3)$ and the schemes don't depend on time and space steps. The high-order quasi-compact scheme of the two-sided tempered fractional diffusion equation needs to be studied.

The novelty of this paper is that the fourth-order numerical schemes of the one-sided tempered fractional diffusion equations are given, and the convergence accuracy is one order higher than that of the existing literature. In the process of proposing the numerical schemes, the left and right third-order Riemann-Liouville tempered derivatives will be generated, so it is necessary to study their effective second-order numerical approximations, and we obtain them through the definition of Riemann-Liouville tempered integer derivatives. Finally, the effective fourth-order numerical schemes for solving these two classes of problems are obtained. In addition, there is a simpler numerical scheme that will be available to solve the two-sided fractional advection-diffusion equations ^[38] by using the idea of this paper.

The remaining sections of this paper are arranged as follows: Section 2 gives the fourth-order quasi-compact approximations of tempered fractional derivatives. In Section 3, the numerical scheme for solving problems (1.1) and (1.2) are derived. The stability and convergence of the numerical schemes are proved in Section 4. Section 5 shows by examples that the numerical schemes are effective. Finally, a brief summary of this work is given in Section 6.

2.   Fourth-order quasi-compact approximations of tempered fractional derivatives

Define a fractional Sobolev space $S_{\lambda}^{n+\alpha}(\mathbb{R})$,

where ${\hat{\nu}}(w) = \int_{\mathbb{R}}\nu(x)e^{-iwx}dx$ is the Fourier transform of $\nu(x).$

Lemma 2.1. ^[25,26] Let $1 < \alpha < 2$, $\lambda\geq0$. The shift number p is an integer, $h$ is the step size, $\nu(x)$ is defined on the bounded interval $[a, b]$, and belongs to $S_{\lambda}^{n+\alpha}(\mathbb{R})$ after zero extension on the interval $x\in(-\infty, a)\cup(b, +\infty)$. The tempered and shifted Grünwald type difference operators are defined as

then

where $g_{k}^{(\alpha)} = (-1)^{k}(_{k}^{\alpha}) (k\geq0)$ denotes the normalized Grünwald weights, that is,

$c_{k}^{\alpha, p}$ are the power series expansion coefficients of the function $W_{p}({s}) = e^{ps}\frac{(1-e^{-s})^{\alpha}}{s} = \sum\limits_{k = 0}\limits^{+\infty}c_{k}^{\alpha, p}s^k,$ and the first four coefficients are given as

Remark 2.1. Let the equations be

then

and there are the following approximations:

where

in particular, when $\alpha = 1,$ the following approximations are given as

For the first left and right Riemann-Liouville tempered derivatives $D^{1, \lambda}_{a, x}\nu(x)$ and $D^{1, \lambda}_{x, b}\nu(x)$, we know

so the normalized left and right Riemann-Liouville tempered fractional derivatives can be rewritten as

and

Let $\Lambda^{\alpha}_{x}: = (I+ch^2D^{2, \lambda}_{a, x})$ and $\Delta^{\alpha}_{x}: = (I+ch^2D^{2, \lambda}_{x, b})$, where $c = \frac{-\alpha^2+\alpha+4}{24}\in(\frac{1}{12}, \frac{1}{6}]$. Now, we give the following theorem as the contribution of this paper.

Theorem 2.1. Let $\nu(x)\in S_{\lambda}^{4+\alpha}(\mathbb{R}), \; 1 < \alpha < 2$, $\lambda\geq0$, and the continuous operators $\Lambda^{\alpha}_{x}$ and $\Delta^{\alpha}_{x}$ are given to operate on the normalized left and right Riemann-Liouville tempered fractional derivatives, respectively. Then,

and

have a fourth-order approximation, respectively. The details are as follows:

and

where

Proof. From Remark $2.1$, the conclusion of $(i)$ and $(ii)$ in Eqs $(2.12)$ and $(2.13)$ can be obtained.

For the conclusion $(iii)$ in Eqs $(2.12)$ and $(2.13)$, it is easy to obtain the following:

Finally, $D^{2, \lambda}_{a, x}(D^{1, \lambda}_{a, x}\nu(x))$ and $D^{2, \lambda}_{x, b}(D^{1, \lambda}_{x, b}\nu(x))$ can be represented as

For derivatives of order 1 to 3, we adopt the following discretizations, respectively,

Thus, the whole conclusion of the theorem is proved. □

3.   Derivation of the numerical schemes

In this section, we consider the numerical scheme for solving problems $(1.1)$ and $(1.2)$. Here, we always assume that the function $u(x, \cdot)$ in problem $(1.1)$ and $(1.2)$ belongs to $S_{\lambda}^{4+\alpha}(\mathbb{R})$ after zero extension.

We make space grid ${\{x_i = a+ih\}}_{0}^{M}$ and time grid ${\{t_n = n\tau\}}_{0}^{N}$, where $h = \frac{b-a}{M}$ and $\tau = \frac{T}{N}$ represent the space step size and the time step size, respectively. Let $u_{i}^{n} = u(x_i, t_n)$, $f_{i}^{n} = f(x_i, t_n)$, and $U_i^{n}$ represent the numerical solution at the point $(x_i, t_n)$.

For problem $(1.1)$, the operator ${\Lambda^{\alpha}_{x}}$ is applied to both the left and right ends of equation. We obtain

then the backward Euler method is carried out at the point $(x_i, t_{n})$ to discrete the time partial derivative, and by Theorem 2.1, we get

where $R_{i}^{n} = O(\tau+h^4)$ is the local truncation error.

Removing the local truncation error in Eq $(3.2)$ yields the numerical scheme as

the matrix form of the numerical scheme $(3.3)$ can be written as

where $P = B^{(\alpha)}+(\alpha-1)\lambda^{\alpha}C^{(\alpha)}-\alpha\lambda^{\alpha-1}B^{(1)}+(\frac{1}{6}-c)\alpha\lambda^{\alpha+2}h^{2}E$, $Q = (\frac{1}{6}-c)\alpha\lambda^{\alpha-1}h^{2}(3\lambda^{2}D^{(1)}+3\lambda D^{(2)}+D^{(3)})$, $U^n = (U_1^n, U_2^n, ..., U_{M-2}^n, U_{M-1}^n)^T$, $f^{n} = (f_1^n, f_2^n, ..., f_{M-2}^n, f_{M-1}^n)^T$, E is the identity matrix of order $M-1$,

and

The same idea is used to solve problem $(1.2)$. The operator ${\Delta^{\alpha}_{x}}$ is applied to both the left and right ends of the equation, and we obtain

then the backward Euler method is carried out at the point $(x_i, t_{n})$ to discrete the time partial derivative. By Theorem $2.1$, we get

where $\hat{R}_{i}^{n} = O(\tau+h^4)$ is the local truncation error.

Removing the local truncation error in Eq $(3.14)$ yields the numerical scheme as

the matrix form of the numerical scheme $(3.15)$ can be written as

where $\hat{P} = \hat{B}^{(\alpha)}+(\alpha-1)\lambda^{\alpha}\hat{C}^{(\alpha)}-\alpha\lambda^{\alpha-1}\hat{B}^{(1)}+(\frac{1}{6}-c)\alpha\lambda^{\alpha+2}h^{2}E$, $\hat{Q} = (\frac{1}{6}-c)\alpha\lambda^{\alpha-1}h^{2}(-3\lambda^{2}D^{(1)}+3\lambda D^{(2)}-\hat{D}^{(3)})$, $\hat{B}^{(\alpha)} = (B^{(\alpha)})^{T}$, $\hat{B}^{(1)} = (B^{(1)})^{T}$, $\hat{C}^{(\alpha)} = (C^{(\alpha)})^T$, $\hat{P} = P^T$, $\hat{D}^{(3)} = -(D^{(3)})^T$, $\hat{Q} = Q^T$,

4.   Stability and convergence analysis of the numerical schemes

In this section, we use the energy method to prove in detail that numerical schemes $(3.3)$ and $(3.15)$ are valid. Before doing so, we first introduce some lemmas that will be used.

Lemma 4.1. ^[40] A real matrix $A$ of order $M$ is positive definite if and only if $D = \frac{A+A^T}{2}$ is positive definite.

Lemma 4.2. ^[41] Let $T$ be a Toeplitz matrix with the generating function $f(x)$ being a $2\pi$-periodic continuous real-valued function. Denote $\lambda_{min}(T)$ and $\lambda_{max}(T)$ as the smallest and largest eigenvalues of $T$, respectively. Then, we have

where $f_{min}(x)$ and $f_{max}(x)$ denote the minimum and maximum values of $f(x),$ respectively. In particular, if $f(x)$ is a nonpositive function and is not always zero, and $f_{min}(x)\neq f_{max}(x)$,

Lemma 4.3. ^[41] (Weyl's theorem). Let $A, E \in \mathbb{C}^{n\times n}$ be a Hermitian matrix and the eigenvalues $\lambda_{i}(A), \lambda_{i}(E), \lambda_{i}(A + E)$ be arranged in an increasing order. Then, for each k = 1, 2, ..., n, we have

Lemma 4.4. For $1 < \alpha < 2, \; 0 < \lambda h\leq1,$ then the matrices P, Q, and $C^{(\alpha)}$ in Eq $(3.4)$ have the following properties:

and the matrix P is negative definite.

$(ii).$ The matrix Q is negative definite.

$(iii).$ For all given nonzero real column vectors $\epsilon$, $C^{(\alpha)}$ satisfies that

Proof. $(i).$ By simple calculation, it is easy to obtain the properties of the elements of the matrix P; From the Gerschgorin disk theorem ^[42], we get that the matrix P is negative definite.

$(ii).$ Note that $Q = (\frac{1}{6}-c)\alpha\lambda^{\alpha-1}h^{2}(3\lambda^{2}D^{(1)}+3\lambda D^{(2)}+D^{(3)}),$ so we just consider $3\lambda^{2}D^{(1)}+3\lambda D^{(2)}+D^{(3)}.$ First, it is easy to show that the matrix $3\lambda^{2}D^{(1)}+3\lambda D^{(2)}$ is negative definite; and then for the matrix $D^{(3)}$, let $D = \frac{D^{(3)}+(D^{(3)})^{T}}{2},$ which is a Toeplitz matrix, and the generating function ^[41] is

From Lemma $4.2$, we can see that the matrix $D$ is negative definite. Further from Lemma $4.1$, the matrix Q is negative definite.

$(iii).$ Let $H = \frac{C^{(\alpha)}+(C^{(\alpha)})^{T}}{2},$ and the generating function of the matrix $H$ is

From Lemma $4.2$, we obtain $1-c(e^{\lambda h}+e^{-\lambda h}+2)\leq\lambda(H)\leq1+c(e^{\lambda h}+e^{-\lambda h}-2),$ and it is easy to check $\frac{1}{12} < \lambda(H) < \frac{4}{3},$ which means $\frac{1}{12}\epsilon^T \epsilon < \epsilon^TH\epsilon < \frac{4}{3}\epsilon^T \epsilon$, that is, $\frac{1}{12}\epsilon^T \epsilon < \epsilon^T{C^{(\alpha)}} \epsilon < \frac{4}{3}\epsilon^T \epsilon$.

Thus, all the conclusions of the lemma are proved. □

Lemma 4.5. ^[43] Assume that $\{k_n\}$ and $\{p_n\}$ are nonnegative sequences, and the sequence $\{\phi_n\}$ satisfies

where $g_{0}\geq0.$ Then, the sequence $\{\phi_n\}$ satisfies

To prove the stability and convergence of numerical schemes by the energy method, we define $\textbf{U}_h = \{\textbf{u}|\textbf{u} = \{\textbf{u}_i\}$ as a grid function defined on $\{x_i = a+ih\}_{i = 1}^{M-1} \; \},$ for $\textbf{u}\in \textbf{U}_h,$ and the corresponding discrete $L_2$-norm is defined as $\|\textbf{u}\|_{L_2} = (h\sum\limits_{i = 1}^{M-1}\textbf{u}_i^2)^{1/2}$. Next, we present the theoretical analysis.

Theorem 4.1. For $\alpha\in(1, 2)$, let $0 < \lambda h\leq1$, then the numerical scheme $(3.3)$ is stable for solving problem $(1.1)$.

Proof. Let $U_{i}^{n}$ and $V_{i}^{n}$ represent the solution obtained by solving problem $(1.1)$ using scheme $(3.3)$ from different initial values. Denoting $\varepsilon_i^n = U_i^n-V_i^n, \; \varepsilon^n = (\varepsilon_1^n, \varepsilon_2^n, ..., \varepsilon_{M-1}^n)^T,$ it can be seen from Eq $(3.3)$ that

Multiplying the left side of both sides of Eq $(4.1)$ by $h(\varepsilon^n)^T$, we obtain

From Lemma $4.3$, we know that matrix $P+Q$ is negative definite, so

which implies that

Finally, it follows from Lemma $4.4$ that

that is,

At this point, we have proved that the numerical scheme $(3.3)$ is stable. □

Theorem 4.2. For $\alpha\in(1, 2)$, let $0 < \lambda h\leq1$, then the numerical scheme $(3.3)$ is convergent for solving problem $(1.1)$, that is, the following relationship is satisfied:

where $e^n = (e_1^n, e_2^n, ..., e_{M-1}^n)^T$, $e_{i}^n = u_{i}^n-U_{i}^n$, and $C_1$ is a constant independent of n, $\tau$, and h.

Proof. Subtracting $(3.3)$ from $(3.2)$ and from scheme $(3.4)$, we obtain

where $R^n = (R_1^n, R_2^n, ..., R_{M-1}^n)^{T}$.

Let's multiply both sides of $(4.6)$ by $h(e^{n})^T,$ which can be written as

Denoting $E^{n} = h(e^{n})^TC^{(\alpha)}e^{n}$, similar to the proof of Theorem $4.1$,

Noticing that $\|e^0\|_{L_2}^2 = 0,$ the Eq $(4.8)$ can be rearranged as

By the Lemma $4.5$ (discrete $Gronwall$ inequalities), we obtain

that is,

Thus, the theorem is proved. □

By the same idea, we can obtain the following theorem, whose proof we skip here.

Theorem 4.3. For $\alpha\in(1, 2)$, let $0 < \lambda h\leq1$, then the numerical scheme $(3.15)$ is stable for solving problem $(1.2)$.

Theorem 4.4. For $\alpha\in(1, 2)$, let $0 < \lambda h\leq1$, then the numerical scheme $(3.15)$ is convergent for solving problem $(1.2)$, that is, the following relationship is satisfied:

where $e^n = (e_1^n, e_2^n, ..., e_{M-1}^n)^T$, $e_{i}^n = u_{i}^n-U_{i}^n$, and $C_2$ is a constant independent of n, $\tau$, and h.

5.   Numerical experiments

In this section, we present some numerical experiments to verify the convergence accuracy of the numerical schemes. Let

be the order of observation.

Example 5.1. Consider the following left Riemann-Liouville tempered fractional diffusion equation with initial-boundary value problem

where $1 < \alpha < 2$, and the linear source term is

and the exact solution is $u(x, t) = e^{t-\lambda x}x^{6}.$

Set the relationship between time step size and space step size as $\tau = h^4$ for all numerical experiments. Choose different $\alpha$, $\lambda$, and space step size $h$, and use numerical scheme $(3.3)$ to solve Example 5.1. The error and observation order results obtained are shown in Table 1. From Table 1, it can be seen that the space convergence order reaches the fourth order, and the numerical results are in perfect agreement with the theoretical analysis.

Example 5.2. Consider the following right Riemann-Liouville tempered fractional diffusion equation with initial-boundary value problem

where $1 < \alpha < 2$, and the linear source term is

and the exact solution is $u(x, t) = e^{t+\lambda x}(1-x)^{8}.$

The numerical scheme $(3.15)$ is used to solve Example 5.2, and the obtained numerical results are presented in Table 2 for different $\alpha, \; \lambda$, and space step size $h$. From Table 2, we find that the space convergence order is fourth order, and the numerical results are consistent with the theoretical results.

6.   Conclusions

In this paper, we focus on fourth-order numerical algorithms for one-sided tempered fractional diffusion equations. By using the quasi-compact technique, the fourth-order quasi-compact approximations of the tempered fractional derivatives and the effective second-order approximations of the third-order tempered derivatives are given. It is proved that the numerical schemes are stable and convergent. Finally, some experiments are given to demonstrate the effectiveness of the numerical schemes.

Acknowledgments

The author would like to express thanks to the editors and anonymous reviewers for their valuable comments and suggestions. This research is supported by the Xinjiang University of Political Science and Law.

Conflict of interest

The author has no conflict of interest.