The alternating direction implicit difference scheme and extrapolation method for a class of three dimensional hyperbolic equations with constant coefficients

1.   Introduction

The class of three-dimensional hyperbolic partial differential equations (PDEs) with constant coefficients studied in this paper is as follows:

where $a$, $b$, and $c$ are constants, $\Omega = (0, L_{1})\times(0, L_{2})\times(0, L_{3})$, $\Gamma$ is the boundary of $\Omega$, and when $(x, y, z)\in\Gamma$, $\alpha(x, y, z, 0) = \varphi(x, y, z)$, $\partial\alpha(x, y, z, 0)/\partial t = \psi(x, y, z)$.

The hyperbolic partial differential equations (1.1)–(1.3) hold substantial theoretical and practical significance, serving as a cornerstone for both the evaluation of numerical methods and the simulation of complex natural and physical phenomena. These models provide a rigorous framework for assessing the feasibility and performance of numerical approaches ^[1,2,3], as well as for modeling intricate dynamic processes observed in nature ^[4,5].

The applications of hyperbolic models span a wide range of disciplines, extending far beyond mathematics. In biology, for instance, the hyperbolic equations are extensively employed to describe population dynamics and the propagation of neural signals ^[6,7,8,9]. The models play a pivotal role in the analysis of vibrational systems and the exploration of fundamental physical principles. For example, Fedotov et al. ^[10] utilized hyperbolic and pseudo-hyperbolic equations to investigate vibration theory, offering insights into the characterization and analysis of vibrational phenomena. Paul ^[11] examined the application of Huygens' principle to hyperbolic equations, while Alinhac ^[12] delved into the blowup behavior of nonlinear hyperbolic equations. The hyperbolic models find significant applications in engineering ^[13,14], where they are used to address a variety of practical problems.

This problem (1.1)–(1.3) is formulated for a three-dimensional wave equation with boundary conditions of the first kind, which has been a subject of considerable research interest in the mixed problem for a hyperbolic equation ^[15,16]. Gordienko ^[17] studied the mixed problem for real wave equations satisfying uniform Lopatinskii conditions, systematically establishing all viable reduction methods to transform it into a mixed problem for symmetric hyperbolic systems with dissipative boundary conditions. Malyshev ^[18] researched initial-boundary value problems for second-order hyperbolic PDEs with complex first-order boundary conditions. These investigations collectively highlight the significance of this research topic in contemporary mathematical physics.

The ADI method has long been recognized as a powerful and efficient numerical technique in scientific computing and engineering, particularly for solving PDEs ^{[19,20,21,22]}. Its unique ability to decompose high-dimensional problems into a series of one-dimensional problems has made it a preferred choice for tackling complex systems, especially those involving high dimension problems in space. The model in Eqs (1.1)–(1.3) has been widely investigated using ADI methods ^{[23,24,25,26,27]}, demonstrating both theoretical significance and ADI unique advantages for high-dimensional simulations of the model. Nie and Cao ^[28] developed an improved ADI scheme based on a staggered grid to address transient heat conduction problems in mixed representations. Zhao ^[29] proposed a novel Douglas ADI method for solving two-dimensional heat equations with interfaces, further demonstrating the versatility of ADI techniques in handling complex boundary conditions.

Over the years, the difference and ADI difference method has garnered significant attention from researchers due to its computational efficiency and stability, which are critical for solving large scale problems in various fields ^{[30,31,32,33,34]}. Zhou et al. ^[35] employed an efficient ADI approach to solve three-dimensional heat equations with irregular boundaries and interfaces. Shen et al. ^[36] provided the high-order ADI difference method and extrapolation method for 2D nonlinear parabolic evolution equations. Chen et al. ^[37] investigated numerical methods for two-dimensional integro-differential equations with two fractional Riemann-Liouville (R-L) integral kernels. Liu et al.^[38] utilized the ADI method to transform three-dimensional problems into one-dimensional problems, enabling the efficient solution of three-dimensional integro-differential equations with weakly singular kernels. Moreover, the ADI method has found extensive use in solving nonlocal evolution equations ^{[39,40,41,42,43]}, it demonstrates the wide adaptability of ADI across different models.

The Richardson extrapolation scheme plays a significant role in refining numerical solutions and achieving higher-order accuracy, a widely applied technique in the numerical solution of linear differential equations ^{[44,45,46,47]} and nonlinear differential equations ^{[48,49,50,51]}. In the present study, we explore the ADI method and its Richardson extrapolation scheme for solving a class of three-dimensional constant-coefficient hyperbolic partial differential equations. By leveraging the ADI method within the framework of the original hyperbolic equation, we decompose the three-dimensional problem into a sequence of one-dimensional subproblems, solving them iteratively.

The main contributions of this article are as follows.

● We prove strictly the solvability, stability and convergence of the proposed ADI difference schemes for the three-dimensional problem, which is the first time to provide a strict theoretical analysis of the three-dimensional problem with constant coefficients. Although there were some papers on the ADI method, but they did not provide strict theoretical analysis for three-dimensional problems.

● We construct the Richardson extrapolation for the original model, present the theoretical analysis, and compare the extrapolated results with the initial data, which fill the gap in the ADI extrapolation method for three-dimensional problems.

The paper is organized as follows. Section 2 introduces notations and lemmas. Section 3 presents the discretization and finite difference scheme. Section 4 proves the solvability, stability and convergence. Section 5 discusses the Richardson extrapolation scheme. Section 6 provides numerical experiments to validate the method. Finally, Section 7 concludes with a summary.

2.   Notations and lemmas

Take positive integers $m_{1}$, $m_{2}$, $m_{3}$, and $n$. Denote $h_{1} = L_{1}/m_{1}$, $h_{2} = L_{2}/m_{2}$, $h_{3} = L_{3}/m_{3}$, and $\tau = T/n$. Define $x_{i} = ih_{1}(0\leq i\leq m_{1})$, $y_{j} = jh_{2}(0\leq j\leq m_{2})$, $z_{k} = kh_{3}(0\leq k\leq m_{3})$, and $t_{l} = l\tau(0\leq l\leq n)$. Define

Let $u = \{u_{ijk}^{l}\mid0\leq i\leq m_{1}, 0\leq j\leq m_{2}, 0\leq k\leq m_{3}, 0\leq l\leq n\}\}$ be a grid function defined on $\Omega_{h}\times\Omega_{\tau}$, and we hereby introduce the following notations:

Let $v^{l} = \{v_{ijk}^{l}\mid(i, j, k)\in\omega\}$ be a grid function on $\Omega_{h}$. Assume the spatial grid $V_{h} = \{v\mid v = \{v_{ijk}^{l}\mid(i, j, k)\in\bar{v}\}\in\Omega_{h}$, $\mathring{V}_{h} = \{v\mid v = \{v_{ijk}\mid(i, j, k)\in\bar{\omega}\}$ and $v_{ijk} = 0$ if $(i, j, k)\in\gamma\}$.

Let $\iota, \mu\in \mathring{V}_{h}$ and introduce the following inner product and norm:

Lemma 2.1. ^[52,53,54] For arbitrary $\iota, \mu\in \mathring V_{h}$, it always follows that

Lemma 2.2. ^[53] Suppose that $g\in C^{(4)}[t_{l}, t_{l}+\tau]$. Then, we have

Lemma 2.3. ^[53] If the function $g\in C^{4}[t_{l}-\tau, t_{l}+\tau]$, then

Lemma 2.4. ^[53] Assume $v\in\mathring{V}_{h}$. It follows that

wherein, $C$ denotes a positive constant that may differ across various scenarios. It invariably hinges upon the solution and the prescribed data, while remaining independent of the time step $\tau$ and the grid spacing $h$.

Lemma 2.5. ^[53] (Gronwall's inequality-E) Let $\{F^{k}\}_{k = 0}^{\infty}$ be a non-negative sequence, $\{g^{k}\}_{k = 0}^{\infty}$ be non-negative and monotonically increasing (not necessarily strictly monotonically increasing), satisfying the condition that

Then, we have

Lemma 2.6. Assume that $v\in \mathring{V_{h}}$. Then, we conclude that

Proof. For $1\leq i\leq m_{1}-1, 1\leq j\leq m_{2}-1, and 1\leq k\leq m_{3}-1$, it holds that

Square the above two equations respectively, and then apply the Cauchy-Schwarz inequality. We can obtain that

Multiply $(2.1)$ by $L_{1}-x_{i}$, multiply $(2.2)$ by $x_{i}$, and then add the obtained results together to get

By applying

it can be derived that

and, following the same argumentation, one arrives at the conclusion that

From $|v|_{1} = \sqrt{\|\delta_{x}v\|^{2}+\|\delta_{y}v\|^{2}+\|\delta_{z}v\|^{2}}$ and the above three equations, we can obtain that

That is,

Moreover, because

then

□

3.   Establishment of the ADI difference scheme

Define

Consider Eqs (1.1)–(1.3) on $(x_{i}, y_{j}, z_{k}, t_{1})$,

From the combination of (1.1) and (1.2), we get

Based on Lemma 2.2, it follows that

where $\psi_{ijk} = \psi(x_{i}, y_{j}, z_{k})$, $\rho_{ijk} = \rho(x_{i}, y_{j}, z_{k})$, $\zeta_{ijk}\in(0, 1)$.

Using Lemma 2.3, we have

Next, substitute the above four formulas into $(3.1)$ and add $M_{1}U_{ijk}^{1}$ on both sides at the same time, so there is

where

Then, there exists a constant $c_{1}$ such that

Consider Eqs (1.1)–(1.3) on $(x_{i}, y_{j}, z_{k}, t_{l})(1\leq{l}\leq{n-1})$,

From Lemma 2.3, it can be concluded that

Then, substitute the above four formulas into (3.4), and add $M_{1}U_{ijk}^{\bar{l}}$

where

A constant $c_{2}$ exists such that

and

Notice the initial boundary value condition (1.1)–(1.3):

Ignoring the minimum terms $(R_{1})_{ijk}^{0}$ and $(R_{1})_{ijk}^{l}$, the following ADI difference scheme is established for the solution of problem (1.1)–(1.3):

Equation (3.8) can be transformed into

Letting

then

For any fixed $j(1\leq j\leq m_{2}-1)$ and $k(1\leq k\leq m_{3}-1)$, take the boundary conditions

Using (3.15) and (3.16), solve the equation

from which we obtain $\{\bar{u}_{ijk}\mid1\leq i\leq m_{1}-1\}$.

For arbitrary fixed $i(1\leq i\leq m_{1}-1)$ and $k(1\leq k\leq m_{3}-1)$, take the boundary conditions

Using (3.17) and (3.18), we solve the equation

from which we obtain $\{u_{ijk}^{*}\mid1\leq j\leq m_{2}\}$.

For fixed $i(1\leq i\leq m_{1}-1)$ and $j(1\leq j\leq m_{2}-1)$, with the boundary conditions

we solve the equation

from which we obtain $\{u_{ijk}^{1}\mid1\leq k\leq m_{3}-1\}$.

Also, Eq (3.9) can be rewritten as

Letting

then

When the values of the $l$-th layer and the $l-1$-th layer are known, for any fixed $j(1\leq j\leq m_{2}-1)$ and $k(1\leq k\leq m_{3}-1)$, take the boundary conditions

Using (3.24) and (3.25), we solve the equation

from which we obtain $\{\hat{u}_{ijk}\mid1\leq i\leq m_{1}-1\}$.

With $i(1\leq i\leq m_{1}-1)$ and $k(1\leq k\leq m_{3}-1)$ being any fixed values, take the boundary conditions

Using (3.26) and (3.27), we solve the equation

from which we obtain $\{\hat{u}_{ijk}^{*}\mid1\leq j\leq m_{2}-1\}$.

For fixed $i(1\leq i\leq m_{1}-1)$ and $j(1\leq j\leq m_{2}-1)$, with the boundary conditions

we solve the equation

from which we obtain $\{u_{ijk}^{\bar{l}}\mid1\leq k\leq m_{3}-1\}$. Then, we utilize

from which we obtain $\{u_{ijk}^{l+1}\mid(i, j, k)\in\omega\}$.

4.   Solvability, stability and convergence

Theorem 4.1. (Solvability) The solution of the ADI difference scheme (3.8)–(3.11) exists and is unique.

Proof. Note $u^{l} = \{u_{ijk}^{l}\mid(i, j, k)\in\omega\}$. It is known from (3.10) and (3.11) that $u^{0}$ has been given. The difference scheme regarding $u^{1}$ is

Consider its homogeneous system of equations

Using $u^{1}$ and (4.1) to calculate the inner product, we obtain

so

and it is easy to get

Now that $u^{l-1}, u^{l}(1\leq{l}\leq{n-1})$ have been determined, the difference scheme for $u^{l+1}$ is

Consider its homogeneous system of equations

Using the equation with $u_{ijk}^{l+1}$ and above for the inner product, we can obtain

Thus,

It is easy to conclude

Based on the induction principle, the system of difference equations (3.8)–(3.11) is uniquely solvable. □

The following is the analysis and proof process regarding stability.

Theorem 4.2. (Stability) Let $\{u_{ijk}^{l}|(i, j, k)\in\bar{\omega}\}$ be the solution of the following system of difference equations:

Then, there is

where

Proof. Taking the inner product of both sides of Eq $(4.2)$ and $\delta_{t}u^{\frac{1}{2}}$, we obtain

From the above equation, we obtain

When $u\in{V_{h}}$, we have

because

Then,

thus $|u|_{A1}^{2}$ and $|u|_{1}^{2}$ are equivalent. If $C_{d}^{2} = \frac{C^{2}}{C_{3}^{2}}$ is set, then it can be inferred from Lemma 2.4 that

After substituting $|u|_{A1}^{2}$ into (4.3) and simplifying both sides of the equation, we obtain

Then, let

so

Taking the inner product of both sides of Eq (4.2) with respect to $\Delta_{t}u^{l}$, we have

through calculation and simplification, we can obtain

thus, it follows that

and when $l = 1$,

so

When $l\geq2$, replace $m$ with $l$ in ($4.5$), and sum $m$ from 1 to $l$. Then, we obtain

observing that

we can obtain

Furthermore, we can obtain

According to Lemma 2.5, we can obtain

moreover, since

then

By noting the definition of $E^{1}$ and ($4.4$), we can see that the above formula also holds for $l = 0$. The proof of stability is completed. □

The rigorous proof of convergence is derived as follows.

Theorem 4.3. (Convergence) Let $\{u(x, y, z, t)\mid(x, y, z)\in\bar{\Omega}, 0\leq t\leq T\}$ be the solution of the problem (1.1)–(1.3), and $\{u_{ijk}^{l}\mid(i, j, k)\in\bar{\omega}, 0\leq l\leq n\}$ be the solution of the difference scheme (3.8)–(3.11). Denote

Then, we have

Proof. Subtracting (3.2), (3.5), and (3.7) from $(3.8)$–$(3.9)$, we can obtain the following error equation system:

From the conclusions in the proof of stability, along with (4.8) and (4.9), we can obtain

Thus

Then, according to Lemma 2.6, we can obtain

The proof of convergence is completed. □

5.   Richardson extrapolation scheme

Theorem 5.1. Assume that the problems

and

have smooth solutions, where

Then we have

Proof. By subtracting $M_{1}e_{ijk}^{1}$ and $M_{1}e_{ijk}^{l}$ from both sides of the error equation system (4.6)–(4.9), we can obtain

where

Denote

Multiply $(5.1)$ by $\tau^{2}$, $(5.2)$ by $h_{1}^{2}$, $(5.3)$ by $h_{2}^{2}$, and $(5.4)$ by $h_{3}^{2}$, respectively. Then, add the resultant expressions to $(5.5)$, yielding

As can be seen from the proof of stability in Subsection 4.2,

that is

By transposing the terms, we can obtain

By the same token, we have

Multiply both sides of (5.6) by $4/3$ and both sides of (5.7) by $1/3$, respectively. Subsequently, subtract the expression obtained from the latter operation from that of the former. As a result, we arrive at

The proof of the theorem is completed. □

6.   Numerical example

In this section, we present two numerical examples to compute the errors and convergence orders using the ADI finite difference scheme. The overall convergence accuracy under extrapolation methods is also calculated, which validates the theoretical analysis.

We take $h_{1} = h_{2} = h_{3} = h$, and denote the maximum absolute error as

The maximum absolute error of the Richardson extrapolation method is

6.1.   Example 1

We consider

The exact solution of this problem is $u(x, y, z, t) = e^{\frac{x+y+z}{2}-t}$.

Tables 1 and 2 show the maximum errors and convergence orders for different step sizes using the finite difference method and for Richardson extrapolation scheme.

As is evident from the aforementioned table, with the progressive refinement of the mesh partitioning, the maximum error exhibits a gradual decline. The precision of the initial finite difference scheme is of the second order. In contrast, following the implementation of the Richardson extrapolation technique, the precision is elevated to the third order.

Figures 1–4 presented herein are the numerical solutions corresponding to diverse step sizes, the exact solutions, and the error graphs derived from the numerical experiments of Example 1 with $1/h = 1/\tau = 16, 32, 64,128$.

6.2.   Example 2

We consider

The exact solution is $u(x, y, z, t) = t^{2}(sinx+siny+sinz)$.

The following two tables shows the maximum error and convergence order of computational Example 2 under different step sizes, as well as the maximum error and convergence order for the extrapolation scheme.

As is discernible from the aforementioned Tables 3 and 4, for Example 2, the errors and convergence orders corresponding to diverse meshing scenarios are elucidated. Evidently, within the context of the original finite difference scheme, the error ratio hovers approximately around order 2. Subsequent to the utilization of the Richardson extrapolation technique, the convergence order approximates order 3. In addition, paralleling the situation in Example 1, as the mesh refinement progresses, the maximum error consistently exhibits a decreasing trend.

In Figures 5–8 the numerical solutions, exact solutions, and error graphs of Example 2 are presented under different step size ratios with $1/h = 1/\tau = 16, 32, 64,128$.

7.   Conclusions

We investigate the ADI difference scheme and the Richardson extrapolation scheme for a class of three-dimensional hyperbolic equations. By employing Lemma 2 and the central difference method, we discretize the equation (1.1)–(1.3) and construct an ADI finite difference scheme. We prove solvability, stability, and convergence, and it is demonstrated that the convergence order is $O(\tau^{2}+h_{1}^{2}+h_{2}^{2}+h_{3}^{2})$. By further studying the Richardson extrapolation scheme for the model, the accuracy is significantly enhanced, and the convergence order is improved.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors are grateful to the editor and the anonymous reviewers for their careful readingand many patient checking of the whole manuscript. The work was supported by National Natural Science Foundation of China Mathematics Tianyuan Foundation (12226337, 12226340), Scientific Research Fund of Hunan Provincial Education Department (24A0422), and Hunan Provincial Natural Science Foundation of China (2024JJ7146).

Conflict of interest

The authors declare there is no conflicts of interest.