The ADI difference and extrapolation scheme for high-dimensional variable coefficient evolution equations

1.   Introduction

Variable-coefficient parabolic partial differential equations (PDEs) are widely used in physics and engineering to model processes such as heat conduction, fluid flow, and diffusion. These equations provide a more accurate representation of real-world phenomena involving spatial and temporal variations in material properties, compared to constant-coefficient equations. With advancements in computational power and numerical methods, the study of these equations has become more efficient ^[1,2].

Parabolic partial differential equations with variable coefficients arise in a wide range of physical and engineering contexts, such as heat conduction in heterogeneous media, diffusion in porous materials, and anisotropic fluid flows. They are fundamental models for describing time-evolving processes involving spatial diffusion and varying material properties.

In this paper, we study the initial boundary value problem for a two-dimensional (2D) variable coefficient parabolic evolution equation

The associated initial and boundary conditions are given by

where $a$ is a positive constant, $\Delta$ is the 2D Laplace operator, $\Omega = \left(0, L_x \right) \times \left(0, L_y \right)$, and $\Gamma$ represents the boundary of $\Omega$. For any point $\left(x, y\right) \in \Gamma$, we have $\psi \left(x, y, 0 \right) = \varphi \left(x, y \right)$. $r\left(x, y, t \right)$, $\varphi \left(x, y \right)$, $\psi \left(x, y, t \right)$ are all smooth functions and $r\left(x, y, t \right) > 0$. Mathematically, the problem (1.1)–(1.2) represents a class of second-order parabolic PDEs with variable coefficients. The variable coefficient $r(x, y, t)$ introduces significant challenges in both theoretical analysis and numerical simulation due to its spatial and temporal dependencies. Solving such problems is essential for accurately modeling real-world phenomena involving nonuniform media and time-dependent material properties.

The variable coefficient term $r(x, y, t)$ introduces spatial and temporal dependencies, which pose significant challenges for both theoretical analysis and numerical simulation. Solving such problems is crucial for accurately modeling physical processes in nonuniform media and time-dependent material properties ^[3,4]. For example, the modeling of fluid dynamics, solitons, and wave propagation ^[5,6], as well as epidemic modeling ^[7,8], rely on such equations.

Variable coefficient parabolic equations play a crucial role in science and engineering, describing physical processes whose behavior changes with both spatial and temporal variables. These equations model complex phenomena such as heat conduction, fluid flow, and diffusion, where coefficients vary over space and time, offering a more accurate representation of real-world non-homogeneous media and transient processes compared to constant-coefficient equations. Recent advancements in computational power and numerical methods have significantly improved our ability to study these equations ^{[9,10,11,12,13]}. Various numerical methods have been developed for solving variable coefficient parabolic equations, including finite element methods ^{[14,15,16,17]}, finite difference methods ^{[18,19,20,21,22]}, finite volume methods ^{[23,24,25,26,27]}, and spectral methods ^{[28,29,30,31]}, each contributing substantially to practical applications.

In recent years, alternating direction implicit (ADI) methods and extrapolation schemes have been applied to solve high-dimensional variable-coefficient evolution equations, significantly improving computational efficiency and accuracy. Moreover, physics-informed neural networks (PINNs) have proven effective for solving these complex problems, especially when working with sparse data ^[32], offering new insights for tackling these challenges ^[33].

The ADI methods are known for their computational efficiency and stability, and are widely applied to high-dimensional evolution equations. For example, Liao et al. ^[34] employed the discrete energy method to demonstrate that the ADI solution is unconditionally convergent with a second-order rate in the maximum norm. Zhou et al. ^[35] proposed efficient ADI schemes for three-dimensional nonlocal evolution equations with weakly singular kernels, and Li et al. ^[36] developed a linearized ADI compact difference method for nonlinear two- and three-dimensional partial integro-differential equations.

Several researchers have extended ADI methods to fractional and time-evolution partial differential equations. Chen et al. ^{[37,38,39,40]} introduced L1-ADI methods for time-fractional diffusion and reaction-subdiffusion equations. Wang et al. ^[41,42] proposed compact ADI schemes for time-fractional integro-differential and diffusion-wave equations, while Huang et al. ^[43,44] focused on Grünwald-Letnikov and L1-ADI schemes for time-fractional reaction-diffusion models.

Qiu et al. ^[45,46] contributed BDF2 and Crank-Nicolson difference schemes for Volterra integrodifferential equations, while Qiao et al. ^[47,48,49] developed fast ADI methods for multidimensional tempered fractional integro-differential equations, including applications involving Brownian motion. Wang et al. ^[50] introduced a high-order compact exponential ADI scheme for solving 2D fractional convection-diffusion equations.

Additionally, Richardson extrapolation techniques ^{[51,52,53,54]} have been combined with ADI methods to further improve accuracy ^[54]. Shen et al. ^[55] combined ADI and Richardson extrapolation for 2D nonlinear parabolic equations. Other related strategies include Shi et al.'s time-space two-grid interpolation method ^[56], Wang et al.'s high-order difference method ^[57,58,59], and the compact predictor–corrector scheme of Jiang et al. ^[60].

However, research on ADI methods for high-dimensional variable-coefficient parabolic equations remains limited. In this work, we perform a theoretical energy analysis to demonstrate the uniqueness and stability of a proposed ADI scheme, and employ Richardson extrapolation to enhance its accuracy.

The main contributions of this work are as follows:

● We rigorously demonstrate the stability of the proposed scheme and establish its a priori bounds.

● We develop and apply an extrapolation method to enhance the accuracy of the numerical solutions. To the best of our knowledge, this is the first instance of applying the ADI extrapolation technique to this equation.

The structure of the paper is organized as follows: The design of the BE-ADI scheme is presented in Section 2. In Section 3, we address the uniqueness, convergence, and stability of the scheme we previously mentioned. Section 4 is dedicated to the formulation of the extrapolation method. Then the numerical results are provided in Section 5. Finally, we draw conclusions in Section 6.

2.   Establishment of the difference scheme

Let $m_1$, $m_2$, and $n$ be positive integers, $h_1 = L_x/m_1$, $x_p = ph_1$, $0 \le p \le m_1$, $h_2 = L_y/m_2$, $y_q = qh_2$, $0 \le q \le m_2$, $\tau = T/n$, $t_w = w\tau$, $0 \le w \le n$. Denote $\Omega_{\tau} = \left\{ t_w \mid 0 \le w \le n \right\}$, $\Omega_h = \left\{ (x_p, y_q) \mid 0 \le p \le m_1, 0 \le q \le m_2 \right\}$, $\omega = \left\{ (p, q) \mid (x_p, y_q) \in \Omega \right\}$, $\gamma = \left\{ (p, q) \mid (x_p, y_q) \in \Gamma \right\}$, $\overline{\omega } = \omega \cup \gamma$. For a grid function $v = \left\{ v^w\mid 0\le w\le n \right\}$ on $\Omega _{\tau}$, we introduce the notations

Let $U_h = \left\{ v\left| v = \left\{ v_{pq}\mid \left(p, q \right) \in \overline{\omega } \right\} \, \, and\, \, v_{pq} = 0\, \, if\, \, \left(p, q \right) \in \gamma \, \, \right. \right\}$. For any grid function $v\in U_{h}$, we denote

Analogous notations may be employed, for instance, $\delta _{y}^{2}v_{pq}$, $\delta _yv_{p, q+\frac{1}{2}}$, $D_yv_{pq}$, $D_{\bar{y}}v_{pq}$, $\Delta _yv_{pq}^{w}$.

For any grid function $v, u\in U_{h}$, the inner product and the norms are defined as follows:

Some grid functions are defined as

Lemma 2.1 (^[11]). For any grid function $u\left(x, y, t \right) \in C^2\left(\left[0, L_x \right] \times \left[0, L_y \right] \times \left[0, T \right] \right)$, we have

Proof. This result can be derived from [11, p. 7]. □

Lemma 2.2. For any grid function $u\left(x, y, t \right) \in C^4\left(\left[0, L_x \right] \times \left[0, L_y \right] \times \left[0, T \right] \right)$, we have

Proof. This result can be derived from [11, p. 7]. □

Lemma 2.3 (^[11], p. 19). For any grid function $v, u\in U_{h}$, we have

Lemma 2.4 (^[11], p. 19). (Gronwall Inequality-E) Suppose $\left\{ F^w \right\} _{w = 0}^{\infty}$ is a non-negative sequence and let $a_1$ and $a_2$ denote two non-negative constants which satisfy

Then, the following inequality holds:

Lemma 2.5 (^[11], p. 72). For any grid function $v \in U_{h}$, ${\left\| v \right\|_\infty } \leqslant \frac{1}{{12}}\sqrt {3\left({\sqrt 2 + 1} \right){L_x}{L_y}} \left\| {{\Delta _h}v} \right\|.$

Next, we construct the BE-ADI difference format of (1.1)–(1.2).

Equation (1.1) can be rewritten in the following form at the grid point $\left(x_p, y_q, t_w \right)$ :

According to Lemmas 2.1 and 2.2, it can be concluded that

in which

It is evident that a constant $c_1$ exists, satisfying

in which

From the initial condition (1.2), there is

By adding a small perturbation term $\frac{a^2\tau ^2}{r_{pq}^{w}}\delta _{x}^{2}\delta _{y}^{2}D_{\overline{t}}U_{pq}^{w}$ on both sides of (2.2), we have

in which $\left(R_2 \right) _{pq}^{w} = \left(R_1 \right) _{pq}^{w}+\left(R_3 \right) _{pq}^{w}$, and

The Taylor formula with an integral remainder gives that

It follows from the (2.3) that there exists a positive constant $c_2$ such that

in which $D_{\bar{t}}\left(R_2 \right) _{pq}^{w} = \frac{1}{\tau}\left[\left(R_2 \right) _{pq}^{w}-\left(R_2 \right) _{pq}^{w-1} \right]$. By neglecting the small term $\left(R_2 \right)_{pq}^{w}$ in (2.5) and substituting the exact solution $U_{pq}^{w}$ with the numerical solution $u_{pq}^{w}$, we obtain the BE-ADI scheme

and the initial and boundary conditions are

Equation (2.7) can be rewritten as

where

By introducing the identity operator $I$, (2.9) can be rewritten as:

By introducing an intermediate variable $u_{pq}^{*} = \left(r_{pq}^{w}L-a\tau \delta _{y}^{2} \right) u_{pq}^{w}, \ \left(p, q \right) \in \omega$, the difference scheme (2.7)–(2.8) can be resolved using two independent sets of one-dimensional equations.

3.   Theoretical analysis of the difference scheme

In this section, we apply the energy analysis method to rigorously establish the uniqueness, convergence, and stability of the solution for the BE-ADI difference scheme.

3.1.   Uniqueness analysis

Lemma 3.1. For any grid function $u\in U_{h}$, define $m = \frac{a\tau}{r_{pq}}$, where $r_{pq} > 0$ is smooth, then the following inequality holds:

with $C = \underset{\left(x, y \right) \in \varOmega, \ t\in \left[0, T \right]}{\max}\left| \frac{\partial m^2\left(x, y, t \right)}{\partial x^2} \right|$.

Proof.

Combine (3.1) and (3.2) to obtain

It can be obtained from $\left| \delta _xm_{p-\frac{1}{2}, q}^{w}-\delta _xm_{p+\frac{1}{2}, q}^{w} \right|\le Ch_1$ that

and similarly,

Thus

□

Theorem 3.1. (Uniqueness) The difference scheme (2.7)–(2.8) has a unique solution when $C < 1$.

Proof. Denote

From (2.8), it is known that $u^0$ is given. Now, assume $u^{w-1}$ has been determined, and then the following difference scheme can be applied for $u^{w}$ :

Consider the homogeneous system of equations

Taking the inner product of $u^w$ on each side of (3.3), we obtain

By Lemmas 2.3 and 3.1, it can be obtained that

If $C < 1$, then ${\left\| {u^w} \right\|^2} = 0$. By applying the principle of induction, we prove that the proposed difference scheme (2.7)–(2.8) is uniquely solvable. □

3.2.   Stability analysis

Next, we prove the $H^2$ stability for the BE-ADI difference scheme.

Theorem 3.2. Let $\left\{ v_{pq}^{w}\left| \left(p, q \right) \in \bar{\omega}, 0\le w\le n \right. \right\}$ be the solution to the difference scheme

and on the boundary $\partial \Omega _h$, let $v_{pq}^{w} = 0$. Then, we have

Proof. Taking the inner product of both sides of (3.5) with $-\Delta _hD_{\overline{t}}v^{w}$, we obtain

Due to

The above three equations are substituted into (3.7), and then we have

Substituting $l$ for $w$ in (3.8) and summing over $l$ from 1 to $w$ derives

Multiplying both sides of the above equation by $\frac{2\tau}{a}$ and rearranging terms, we obtain

Therefore,

It follows from Lemma 2.4 that

The proof has been concluded. □

3.3.   Convergence analysis

In this section, the convergence in the maximum norm of the proposed scheme (2.7)–(2.8) will be considered.

Theorem 3.3. Suppose $u\left(x, y, t \right) \in C^4\left([0, L_x] \times [0, L_y] \times [0, T] \right)$ is the solution of (1.1)–(1.2), and $\left\{ u_{pq}^{w} \right\}$ is the solution of the difference scheme (2.7)–(2.8). Denote

then we have

Proof. Subtracting (2.4), (2.5) from (2.6) and (2.7), we obtain the error equations

From Theorem 3.2 and (2.6), we obtain

Taking the square root of both sides, we obtain the inequality

From Lemma 2.5, we obtain

The proof of the convergence theorem is complete. □

4.   Richardson extrapolation

Theorem 4.1. Assume that the problems

and

admit smooth solutions, where

Then

Proof. Since

thus, the error equations (3.10) can be written as

Denote

Discretizing (4.1), we obtain

Discretizing (4.2) and (4.3), we obtain

Let

Multiply both sides of (4.4) by $\tau$, multiply (4.5) by $h_1$, and multiply (4.6) by $h_2$. Then, add the results to obtain

By Theorem 3.2, we obtain

which implies

Rearranging terms, we get

In fact, by multiplying both sides of (4.4) by $\tau/4$, multiplying (4.5) by $(\frac{h_1}{2})^2$, and multiplying (4.6) by $(\frac{h_2}{2})^2$, we obtain, in a similar way, that

Multiply Eq (4.9) by $\frac{4}{3}$ and Eq (4.8) by $\frac{1}{3}$, and then subtract the resulting equations to derive

□

5.   Numerical experiments

In this section, we will verify the theoretical results of the BE-ADI difference scheme through two numerical examples. The convergence rates are listed to test the accuracy of the difference scheme (2.7)–(2.8). Here, $E_{\infty}\left(h_1, h_2, \tau \right)$ represents the maximum error at the grid nodes for step sizes $h_1$, $h_2$, and $\tau$. For practical implementation, set $h_1 = h_2 : = h$. Before presenting the numerical examples, we denote the maximum errors and global convergence rates as follows:

Example 1: The following problem is considered:

The exact solution is $u(x, y, t) = e^{\frac{1}{2}(x + y) - t}$, and the variable coefficient is $r(x, y, t) = e^{(x + 1)(y + 1)t}$.

Example 2: Consider the following problem:

The exact solution is $u\left(x, \; y, t \right) = \text{e}^{\frac{1}{2}\left(x+y \right) -t}$, and the variable coefficient is $r\left(x, \; y, t \right) = 2.5\sin \left(\pi x \right) \cdot \sin \left(\pi y \right) \sin \left(\pi t \right) +2.5$.

In the two examples presented, the spatial discretization is chosen as $m = 10, 20, 40, 80$, while the temporal discretization varies depending on the difference scheme employed.

Tables 1 and 2 present the spatial and temporal maximum error estimates, their convergence rates, and computational times for the BE-ADI method described by (2.7)–(2.8) in Example 1. Table 3 shows the corresponding errors and convergence rates when the extrapolation method is applied. The error surface of the numerical solutions at $t = 1$ for various step sizes is shown in Figure 1. The results indicate that the proposed scheme maintains second-order accuracy in space and first-order accuracy in time. However, when the extrapolation method is applied, the accuracy improves significantly, enhancing the convergence rates from second to fourth order in space.

Tables 4 and 5 present the spatial and temporal maximum error estimates, along with their convergence rates and computational times, for the BE-ADI method described by (2.7)–(2.8) in Example 2. As observed, the maximum error decreases exponentially as the grid is refined, and the convergence rate indicates that the scheme achieves second-order accuracy. Table 6 presents the maximum errors and convergence rates of the numerical solutions obtained using the extrapolation method, which achieves fourth-order accuracy. Figure 2 illustrates the error surface at $t = 1$ for numerical solutions computed with various step sizes.

Figures 1 and 2 show that the error distribution is uniform, and the peak values clearly decrease as the step size decreases, indicating that this difference scheme achieves higher accuracy on finer grids.

6.   Conclusions

This paper investigates the BE-ADI scheme and its fourth-order high-accuracy Richardson extrapolation scheme for 2D parabolic equations with variable coefficients. First, we presented the BE-ADI scheme with an accuracy of order $O\left(\tau +h_{1}^{2}+h_{2}^{2} \right)$ and proved the uniqueness, convergence, and stability of this scheme using energy analysis methods. To further improve the numerical accuracy, we established a Richardson extrapolation scheme and theoretically demonstrated that it achieves an accuracy of order $O\left(\tau ^2+h_{1}^{4}+h_{2}^{4} \right)$. Two illustrative examples were provided to validate the theoretical findings, supported by error surface plots. The results demonstrate that the proposed difference scheme and Richardson extrapolation method effectively address the challenges of low accuracy and computational complexity in the numerical solution of variable-coefficient parabolic equations.

Moreover, by reducing the three-dimensional equation into three sets of independent one-dimensional equations, the proposed method can be naturally extended to three-dimensional problems, offering an efficient and scalable approach for higher-dimensional simulations. In future work, we aim to explore the application of this approach to more complex nonlinear problems and variable-coefficient equations with non-smooth coefficients. Further improvements on adaptive mesh refinement and parallelization techniques will also be considered to enhance computational efficiency for large-scale 3D problems. By reducing the three-dimensional equation into three sets of independent one-dimensional equations, this method can also be extended to three-dimensional problems.

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 reading and many patient checking of the whole manuscript. The work was supported by the 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 are no conflicts of interest.