Analysis of a fourth-order compact $\theta$-method for delay parabolic equations

1.   Introduction

Mathematical scientists build partial differential equations to model nature phenomena in science and engineering. To better understand the nature phenomena, not only is a reasonable qualitative analysis about the solutions of the equations necessary (see, e.g., ^[1,2,3,4,5]), but, also, an efficient and accurate numerical algorithm is required (see, e.g., low-dimensional space ^[6,7,8,9,10], high-dimensional space ^{[11,12,13,14]}). Time delay phenomena appear in various fields such as signal transduction, population dynamics, and control systems (see, e.g., ^{[15,16,17,18,19]}). Partial differential equations with delay can describe the dynamical systems more accurately, and thus have been receiving more attention. Generally, obtaining an analytical solution for a delay system is challenging ^[20,21]. Researchers have developed effective methods to solve these equations, including finite difference methods ^{[22,23,24,25,26]}, finite element methods ^{[27,28,29,30,31]} (discontinuous Galerkin method ^{[32,33,34,35]}), and many other methods ^{[36,37,38,39,40]}.

To avoid the propagation of the errors of the numerical solutions and control the numerical process, many studies have been devoted to the stability analysis of numerical methods (see, e.g., ^{[41,42,43,44,45]}). However, the stability analysis of a numerical method for solving partial differential equations with delay has not been well investigated. A numerical method which is unconditionally stable for partial differential equations, is probably no longer unconditionally stable for partial differential equations with delay ^[46,47]. Therefore, the stability analysis of a numerical method for solving partial differential equations with delay is a necessity. Our concern here is to investigate the stability estimates of a fourth-order compact $\theta$-method for the following delay parabolic equations

where $k, r$ are diffusion coefficients, and $\tau$ is delay term. One can refer to Reference ^[48] for the existence, uniqueness and continuation of problem (1.1). The fully discrete fourth-order compact $\theta$-method is constructed by performing the following steps: the spatial discretization is realized by using a compact finite difference method (e.g., ^[49,50]), and the time discretization is realized via a linear $\theta$-method (e.g., ^[51]).

The spectral radius condition is utilized as a tool for the stability of numerical methods. As for problem (1.1), References ^[49,50] investigated the asymptotic stability property by using a spectral radius condition and established the sufficient and necessary condition for the method to be asymptotically stable. For more details, please refer to Theorems 5.6 and 5.7 of Reference ^[49] and Theorems 1 and 2 of Reference ^[50]. The results obtained by using the spectral radius condition are very important in the determination of when numerical methods are asymptotically stable. However, scholars have presented instructive examples to confirm that the spectral radius condition may give an unreliable information about the numerical stability estimates (see, e.g., ^[52,53,54]). Moreover, the examples demonstrate that the norm of the powers of the iteration matrix could be obsessively large in the case that the spectral radius condition is satisfied ^[55]. In view of these facts, it is interesting to consider whether the numerical process will actually behave in a stable way within its stability region ^[56].

The upper bounds for the powers of matrices are intimately connected to the stability analysis of numerical processes for solving initial(-boundary) value problems in ordinary and partial linear differential equations ^[54]. Recently, scholars utilized the Kreiss resolvent condition to establish upper bounds for the growth of errors in the numerical process for solving delay differential equations (see, e.g., ^{[51,53,54,55,56,57,58]}). To the best of our knowledge, stability estimates about upper bounds for the powers of the iteration matrix derived via the compact $\theta$-method for problem (1.1) have not been studied. The main contribution of this paper is that, based on the aforementioned results, we look into the upper bounds for the norm of the nth power of the iteration matrix of the compact $\theta$-method and detect whether propagated errors occurring in the numerical process could be bounded. It is worth noting that our stability results are essentially different from and, in some respects, complementary to those obtained in References ^[49,50].

Reference ^[51] demonstrated that the trivial solution of problem (1.1) is asymptotically stable if $k > r > 0$. Throughout the paper, we investigate the stability of the proposed numerical method under this condition. In addition, we assume that the solution $u(x, t)$ of problem (1.1) satisfies that $u(x, t)\in \mathcal{C}^{(6, 4)} \Big((0, \pi) \times (0, T]\Big)$.

The rest of this paper is organized as follows. In Section 2, we briefly establish the fully discrete compact $\theta$-scheme. In Section 3, we present the solvability, asymptotic stability, and convergence of the compact $\theta$-method. In particular, we derive the upper bounds for the norm of the nth power of the iteration matrix derived via the proposed method, and come to conclusions about the stability properties. In Section 4, we conduct numerical experiments to illustrate the theoretical results. Finally, we summarize the paper in Section 5.

2.   The fully discrete numerical scheme

In this section, we give a brief derivation of the fully discrete compact $\theta$-scheme for solving problem (1.1).

Let $\Delta t = \frac{\tau}{m}$ and $\Delta x = \frac{\pi}{N},$ where $m$ and $N$ are two positive integers. Set $\Omega_{\Delta t} = \{t_n| n \geq\!-m\}$ as a partition on the time interval $[-\tau, \infty),$ where $t_n = n\Delta t \; (n\geq\!-m).$ Set $\Omega_{\Delta x} = \{x_j| 0\leq j\leq N\}$ as a mesh on the space interval $\Omega = [0, \pi],$ where $x_j = j\Delta x \; (0\leq j\leq N)$. Let $u_j^n$ denote the numerical approximation of $u(x_j, t_n)$. We define the following notations:

and introduce the compact operator ^[49,50]

Next, employing the compact operator to approximate the diffusion term, and applying the linear $\theta$-method to the semi-discrete system, the fully discrete compact $\theta$-scheme reads as follows:

Remark 2.1. If the discretization of the spatial diffusion term is achieved by using the standard central difference method, then we get the typical linear $\theta$-method. It is worth noting that the convergence order of the compact $\theta$-method in space is 4, while that for the linear $\theta$-method is 2 ^[49,50]. Moreover, the coefficient matrices for the compact $\theta$-method and linear $\theta$-method are both tridiagonal matrices. Thus, the compact finite difference method shows a better convergence result in space without increasing the computational cost, making it a popular method to approximate the spatial term ^{[59,60,61,62,63]}.

Rewrite the compact $\theta$-method in the following matrix form:

where

and the $(N\!-\!1)$-by-$(N\!-\!1)$ matrix $S$ is given by

with its eigenvalues $\lambda_j = 2 \text{cos}(j\Delta x), \; 1\leq j \leq N\!-\!1.$

3.   Main results

In this section, we will present the solvability, asymptotic stability, and convergence of the compact $\theta$-method. In particular, we will derive the upper bounds for the norm of the nth power of the iteration matrix derived via the compact $\theta$-method. We will prove that, under certain conditions, the numerical process actually behaves in a stable way within its stability region.

3.1.   Consistence and solvability

Theorem 3.1 ^[49,50]. The compact $\theta$-method (2.1) is consistent and has a unique solution.

3.2.   Stability

Definition 3.2 ^[51]. A numerical method applied to problem (1.1) is called asymptotically stable with respect to the trivial solution if its approximate solution $u_j^n$ corresponding to any function $u_0(x, t)$ satisfies

From the recurrence relation in (2.2), we obtain an equivalent vector form

Here

and C is a block matrix of order $(m+1)(N\!-\!1)$, defined as

where the elements of ${\textbf{C}}$ are defined as follows:

and the symbols $I$ and $\mathbb{O}$ stand for the identity and zero blocks of order $N\!-\!1,$ respectively. We shall introduce a vector norm $|\cdot|$ for ${\textbf{U}}_n$ and a matrix norm $\|\cdot\|$ induced by the vector norm $|\cdot|$ for ${\textbf{C}}^n\; (n\geq1)$ (see, e.g., ^[51,55]), defined as

where ${\textbf{U}}_0 = (1, 1, ..., 1)^T\in \mathbb{R}^{(m+1)(N\!-\!1)}.$

In order to investigate the stability property, it is necessary to consider the perturbed problem, i.e., the numerical computation is performed by using a different initial vector $\tilde{{\textbf{U}}}_0.$ Let $\tilde{{\textbf{U}}}_n$ denote the perturbed numerical solutions, and let ${\textbf{W}}_n = \tilde{{\textbf{U}}}_n-{\textbf{U}}_n$ be the propagated errors. Therefore,

Our aim is to seek bounds for the propagated errors, given by $|{\textbf{W}}_{n}|\leq L|{\textbf{W}}_0|\; \; (n\geq1),$ where $|\cdot |$ stands for a vector norm, and $L$ is a constant. Noticing (3.1), we arrive at

If there exist upper bounds on ${ \left\|{\textbf{C}}^n\right\|}$ for every natural number $n$, then the numerical process (2.1) will behave in a stable manner.

To achieve the stability estimate for ${ \left\|{\textbf{C}}^n\right\|}$, we introduce the following:

where the $(m\!+\!1)$-by-$(m\!+\!1)$ matrix $F_j$ is given by

with $\tilde{\lambda}_j = -4 \text{sin}^2(j\frac{\Delta x}{2})\in \sigma(A), \; 1 \leq j\leq N\!-\!1.$

We shall present a useful lemma from Reference ^[55].

Lemma 3.3. For $n, m\geq 1$ and $N\geq 2,$ it holds that

It is observed that $\|F_j^n\|_\infty$ plays a crucial role in the stability estimate. As Zubik-Kowal pointed out in Reference ^[55], the value of $\|F_j^n\|_\infty,$ and thus of $\|{\textbf{C}}^n\|$, might be excessively large, even under the spectral radius condition. We wonder under what conditions the norm of the nth power of the iteration matrix ${\textbf{C}}$ could be bounded. The following analysis will give the answers.

We assume that $(\mu, \nu)$ belongs to the stability region $S_{\theta}$ in the compact $\theta$-method. The stability region $S_{\theta}$ is given by Theorems 1 and 2 in Reference ^[50]:

Theorem 3.4. Assume that $(\mu, \nu)\in S_\theta$ and $\theta = 0, \; 1.$ Then,

Proof. Define

Case 1. For $\theta = 0,$ we have

It follows from

that

Thus,

Noticing that our goal is to show that $|\psi_1^j|+|\psi_m^j|+|\psi_{m+1}^j|\leq1,$ it follows that $\|F_j\|_\infty = 1$ holds. Now, we prove the following two cases.

● If $1+(\frac{1}{12}+\mu)\tilde{\lambda}_j\geq 0,$ then $\psi_1^j\geq0.$ Noting that $\mu > \nu > 0$, we have

● If $1+(\frac{1}{12}+\mu)\tilde{\lambda}_j < 0,$ then $\psi_1^j < 0;$ thus,

For $1\!\leq\! j \!\leq \!N\!-\!1,$ noting that

we find that

For $(\mu, \nu)\in S_\theta$ and $\theta = 0,$ we have

Thus,

Incorporating (3.7) into (3.6), we get

According to the two aforementioned cases, we derive that, for $\theta = 0,$ $\|F_j\|_\infty = 1$ holds for $1\leq j\leq N\!-\!1.$

Case 2. For $\theta = 1,$ we have

Noting that $1+\frac{1}{12}\tilde{\lambda}_j\in (\frac{2}{3}, \; 1),$ $\tilde{\lambda}_j < 0,$ and $\mu > \nu > 0$, we arrive at

Hence, for $\theta = 1$, $\|F_j\|_\infty = 1$ holds for $1\leq j\leq N\!-\!1.$

Now, we derive the conditions such that (3.5) holds for $\theta \in(0, \; 1).$

Theorem 3.5. Suppose that $(\mu, \; \nu)\in S_\theta.$ Then, we can derive the following set of conclusions for $\theta \in(0, \; 1)$.

(ⅰ) For $\theta\in[\frac{\Delta x^2}{12k\Delta t}+\frac{k+r}{2k}, \; 1),$ (3.5) holds for every $N\geq2.$

(ⅱ) For $\theta\in[\frac{1}{2}, \; \frac{\Delta x^2}{12k\Delta t}+\frac{k+r}{2k}),$ and $\Delta t\leq\dfrac{\Delta x^2}{(\frac{\Delta x^2}{6\Delta t}+r)[1+\cos(\Delta x)]},$ (3.5) holds for every $N\geq2.$

(ⅲ) For $\theta\in(0, \; \frac{1}{2}),$ and $\Delta t\leq\dfrac{\Delta x^2}{[\frac{\Delta x^2}{6\Delta t}+(1-2\theta)k+r][1+\cos(\Delta x)]},$ (3.5) holds for every $N\geq2.$

Proof. For fixed $j$ ($1\!\leq\! j\! \leq\! N\!-\!1$), noting the facts that $\tilde{\lambda}_j < 0$ and $1+(\frac{1}{12}-\mu\theta)\tilde{\lambda}_j > 0$, we have

Therefore

Case Ⅰ. If $\psi_1^j\geq0,$ then

Thus, (3.5) holds for every $N\geq 2$.

Case Ⅱ. If $\psi_1^j < 0,$ then

We recall our goal, i.e., $|\psi_1^j|+|\psi_m^j|+|\psi_{m+1}^j|\leq1,$ which yields

i.e.,

● For $\theta\in[\frac{\Delta x^2}{12k\Delta t}+\frac{k+r}{2k}, \; 1),$ we have

Then,

implying that inequality (3.8) holds; hence, Eq (3.5) holds for every $N\geq2.$

● For $\theta\in[\frac{1}{2}, \; \frac{\Delta x^2}{12k\Delta t}+\frac{k+r}{2k}),$ we have

Then,

It follows from

that

Therefore, conclusion (ⅱ) holds.

● For $\theta\in(0, \; \frac{1}{2}),$ we have

Then,

Following from

we have

Therefore, conclusion (ⅲ) holds.

We have discussed the conditions such that the formula $\|F_j\|_\infty = 1$ holds for $1\leq j\leq N\!-\!1$ so far. Noticing Lemma 3.3, we arrive at the following corollary.

Corollary 3.6. Assume that $(\mu, \nu)\in S_\theta$. Then for $n, \; m\geq 1$ and $N\geq 2,$ when

1) $\theta = 0,$ or

2) $\theta\in(0, \; \frac{1}{2}),$ and $\Delta t\leq\dfrac{\Delta x^2}{[\frac{\Delta x^2}{6\Delta t}+(1-2\theta)k+r][1+\cos(\Delta x)]},$ or

3) $\theta\in[\frac{1}{2}, \; \frac{\Delta x^2}{12k\Delta t}+\frac{k+r}{2k}),$ and $\Delta t\leq\dfrac{\Delta x^2}{(\frac{\Delta x^2}{6\Delta t}+r)[1+\cos(\Delta x)]},$ or

4) $\theta\in[\frac{\Delta x^2}{12k\Delta t}+\frac{k+r}{2k}, \; 1],$ it holds that

Remark 3.7. Scholars have illustrated that the spectral radius condition may give unreliable information about the numerical stability estimates (see, e.g., ^[52,53,54]). Reference ^[55] pointed out that the norm of the powers of the iteration matrix could be obsessively large in the case that the spectral radius condition is satisfied. It is interesting to consider whether the numerical process will actually behave in a stable way within its stability region ^[56]. To the best of our knowledge, the upper bounds for the powers of the iteration matrix derived via the compact $\theta$-method have not been studied. Here, we have investigated the conditions for the norm of the nth power of the iteration matrix ${\textbf{C}}$ to be bounded, and such that the numerical process will actually behave in a stable manner.

3.3.   Convergence

Theorem 3.8 (Lax equivalence theorem ^[64]). For the linear finite difference scheme, it is convergent if it is consistent and stable with respect to the initial value.

In light of Subsection 3.1 (consistency) and Subsection 3.2 (stability), and with the help of Lax equivalence Theorem 3.8, we have the following convergent result.

Theorem 3.9 ^[49]. Let $\{u_j^n|, \; j = 1, \; 2, \cdots, N-1, \; n\! = \!-m, \!-m\!+\!1, \cdots\}$ be a numerical solution of the fully discrete scheme (2.1). Denote $U_j^n = u(x_j, t_n)\; (j = 0, 1, \cdots, N, \; n = -m, -m+1, \cdots$). Suppose that the assumptions in Section 1 and Corollary 3.6 hold. Then, for $n = 1, 2, \cdots$, we have

where $e^n = [u_1^n- U_1^n, u_2^n- U_2^n, \cdots, u_{N-1}^n-U_{N-1}^n]^T$ and $\hat{C}$ is a positive constant independent of $\Delta t$ and $\Delta x$.

4.   Numerical simulations

In this section, we present numerical simulations to validate the derived stability results and the convergence results in the maximum norm.

4.1.   Stability test

Let $k = 1, \; r = 0.5, \; \tau = 0.1$ and $\; N = 4.$ To survey whether the norm of the nth power of the iteration matrix ${\textbf{C}}$ could be bounded, we define a function

where $\|{\textbf{C}}^n\|: = \frac{{\left|{\textbf{C}}^n{\textbf{U}}_0\right|}}{{\left|{\textbf{U}}_0\right|}}$ with ${\textbf{U}}_0 = (1, 1, ..., 1)^T\in \mathbb{R}^{(m+1)(N\!-\!1)}.$ If Corollary 3.6 holds, then we have that $\|{\textbf{C}}^n\|\leq \sqrt{\min\{N\!-\!1, m+1\}}$, implying that function $E(m)$ could be bounded. We shall examine the behavior of the function $E(m)$ versus $m$ under different values of $\theta$ in Corollary 3.6. Specifically, we take $\theta_1 = 0,$ $\theta_2 = \frac{1}{4},$ $\theta_3 = \frac{1}{2}$ and $\theta_4 = 1$, corresponding to conclusions (1)–(4) of Corollary 3.6, respectively.

Graphs of the function $E(m)$ versus $m$, for $\theta_1 = 0,$ $\theta_2 = \frac{1}{4},$ $\theta_3 = \frac{1}{2}$ and $\theta_4 = 1$, are shown in Figures 1–4, with $n = 50, \; 75, \; 100,$ respectively. For $\theta_1$, as expected, $E(m)$ is bounded for various values of $n$ in Figure 1. For $\theta_2$, under the condition $\Delta t\leq\frac{\Delta x^2}{[\frac{\Delta x^2}{6\Delta t}+(1-2\theta_2)k+r][1+\cos(\Delta x)]},$ it can be observed in Figure 2 that the function $E(m)$ is bounded for various values of $n$. For $\theta_3,$ under the condition $\Delta t\leq\frac{\Delta x^2}{(\frac{\Delta x^2}{6\Delta t}+r)[1+\cos(\Delta x)]},$ as shown in Figure 3, the function $E(m)$ is also bounded for various values of $n$. Again, $E(m)$ is bounded for various values of $n$ in Figure 4 for $\theta_4.$ Each graph in Figures 1–4 indicates the function $E(m)$ is bounded for various values of $n$, thus the tendency of the function $E(m)$ in each graph looks similar, which is consistent with the results for the $\theta$-method in Reference ^[51]. The graphs demonstrate that, under the conditions of Corollary 3.6, the compact $\theta$-method actually behaves in a stable manner within its stability region.

4.2.   Convergence test

We describe a convergence test for the following equation

where the initial condition and the added term $f(x, t)$ are specified so that the exact solution is $u(x, t) = e^{-t} {\rm{sin}}(x).$

We take the parameters $k = 1, r = 0.5, \tau = 0.5, \theta = 0.5$ and solve the problem on $[0, \pi]\times[0, t]$ with different temporal and spatial step sizes. Let $\Delta t \approx \Delta x^2$ when the compact $\theta$-method is employed to solve the problem. Also, we choose to test the problem by using the linear $\theta$-method. The numerical errors and convergence orders in the maximum norm are listed in Tables 1 and 2, respectively. It is shown that the convergence of the compact $\theta$-method is second-order accurate in time and fourth-order accurate in space, while the linear $\theta$-method has second-order accuracy in time and space. Obviously, the convergence result of the compact $\theta$-method is superior to that of the linear $\theta$-method in the spatial direction.

5.   Conclusions

In this paper, we have discussed the numerical analysis of a fourth-order compact $\theta$-method for solving problem (1.1). We have presented the solvability, asymptotic stability, and convergence of the method. In particular, we have derived the upper bounds for the norm of the nth power of the iteration matrix ${\textbf{C}}$ obtained by applying the method in the corresponding stability region, and we have drawn conclusions about the stability properties, such that the propagated errors occurring in the numerical process could be bounded. Numerical simulations have been conducted to illustrate the theoretical results. There has been much concern regarding the fractional-order model recently since the model could characterize complex nature phenomena more accurately ^{[65,66,67,68,69]}. We hope to apply the method to solve the fractional generalized diffusion equation with delay in the future.

Use of AI tools declaration

The authors declare that they have not used artificial intelligence tools in the creation of this article.

Acknowledgments

The authors are grateful to the editor and the anonymous reviewers for their valuable suggestions and comments, which have significantly improved the quality of the manuscript. This work was supported by the NSFC (Grant No. 12001067) and the PhD Research Startup Fund for Shandong Police College (No. YKYPYZX202307). All authors contributed equally to this paper.

Conflict of interest

The authors declare that there is no conflicts of interest.