Error estimate and superconvergence of a high-accuracy difference scheme for 2D heat equation with nonlocal boundary conditions

1.   Introduction

In recent years, nonclassical boundary and initial-boundary value problems have garnered significant attention across diverse disciplines such as physics, biology, ecology, chemistry, and beyond. Among these, parabolic partial differential equations (PDEs) with nonlocal initial and/or boundary conditions have emerged as powerful tools for modeling a wide array of phenomena. These include, but are not limited to, heat conduction^[1], thermoelasticity^[2], biotechnology^[3], electrochemistry^[4], population dynamics ^[5], and petroleum exploration ^[6]. The incorporation of nonlocal conditions into these PDEs allows for a more nuanced and realistic representation of the complex interactions and dynamics at play within these systems.

Let $Q_T = \Omega\times I$ be the computational domain, where $\Omega = (0, 1)^2$ and $I = (0, T)$ represent the spatial domain and the time domain, respectively, and $T$ is a positive constant. Here, we consider the following 2D parabolic problem to find a high-accuracy numerical scheme and obtain its theoretical error estimates:

which is subject to the initial conditions

the Dirichlet boundary conditions

and the nonlocal boundary conditions

where $u(x, y, t)$ is the unknown function, $g(x, y)$, $\mu_i(y, t)$ $(i = 1, 2)$ and $\mu_j(x, t)$ $(j = 3, 4)$ are known functions, and $a$ is a positive constant.

These two nonlocal boundary conditions (1.5) and (1.6) are often be used to describe the correlation of a physical quantity across two parallel boundaries in a physical system, as well as the situation where the normal derivative at the boundaries is controlled by external factors, which is commonly used to simulate the interactions between boundaries and boundary effects in processes such as heat conduction and fluid flow.

If the exact solution $u$ of problems (1.1)–(1.6) satisfies certain smootheness conditions, then the compatible condition is deduced as follows: $\forall (x, y)\in \Omega$, the following relations hold:

The analytical frameworks and numerical techniques employed in tackling parabolic problems with nonlocal conditions have aroused the concern of many scholars. Pertaining to the crucial aspects of convergence and stability for such problems, we acknowledge the foundational work presented in ^[7,8,9], as well as the extensive references cited therein. Among the prevalent numerical methodologies, finite difference methods (FDM) stand out prominently, with notable contributions from studies such as ^{[7,10,11,12,13]}. Additionally, finite element methods (FEM) have garnered substantial attention, exemplified by works cited in ^[14,15]. Furthermore, the realm of numerical solutions encompasses innovative approaches like Adomian expansions ^[16], the local coordinates method ^[17], and the utilization of reproducing kernel spaces ^[18], each offering unique insights and advancements in this field.

It is widely acknowledged that two-dimensional parabolic partial differential equations (PDEs), characterized by their two spatial variables, pose significant challenges for theoretical analysis, particularly in the realms of convergence analysis and error estimation. The dimensionality of these variables often complicates the mathematical treatment, necessitating innovative strategies. One promising approach to mitigate these difficulties is the utilization of the discrete fourier transform (DFT) method, which offers advantages in reducing the complexity of self-variables during convergence analysis. In this study, we build upon our previous work ^[19,20] by extending the numerical schemes and integrating the DFT method on the spatial variable $x$ for error estimation within the context of a two-dimensional parabolic PDE subject to a nonlocal boundary condition.

However, a major obstacle arises from the complex boundary condition imposed on the spatial variable $y$. This condition presents a challenge to traditional DFT methods, which are inherently designed to preserve some boundary conditions. To overcome this limitation, we propose a novel transformation tailored specifically to handle this periodic boundary scenario. Furthermore, we contribute by deriving formulas for the solution derivatives and rigorously proving that these formulas enable us to achieve optimal asymptotic error estimates in the maximum norm. This achievement underscores the effectiveness and applicability of our proposed methodology in accurately approximating and analyzing solutions to two-dimensional parabolic PDEs with intricate boundary conditions.

This paper is organized as follows. In Section 2, the backward Euler difference scheme for the solution of problems (1.1)–(1.6) is presented. Then, in Section 3, we utilize the DFT and develop a new transformation to analyze the error estimate for the corresponding difference equation. The superconvergence for the derivative and its theoretical results are also considered. Finally, some numerical experiments are presented in Section 5.

2.   Finite difference discretization

Now, we use the FDM to discretize problems (1.1)–(1.6). The domain $\overline{Q}_T$ is discretized by the uniformly distributed grid points $(x_i, y_j, t_n)$, where

where $\tau$ is time stepsize, and $h$ is space stepsize along both $x$ and $y$ directions.

Define a function space by

and its norm by

where $m$ and $s_i$ $(i = 1, 2, 3)$ are given nonnegative integers.

The key to seeking a numerical solution for problems (1.1)–(1.6) lies in how to discretize the nonlocal boundary conditions (1.6). Suppose $u\in C^{4}(\overline{Q}_T)$, using the Taylor formula, we have

Using (1.1), we have

Moreover, we obtain

Therefore, with (1.6), we obtain

Substituting (1.6), (2.2), and (2.3) into (2.1), we have

i.e.

where

From the derivation process described above, the discretization of (1.6) is converted to discretizing (2.4).

Let $u^n_{i, j}$ and $U^n_{i, j}$ be the exact value and the approximation of $u(x, y, t)$ at grid point $(x_i, y_j, t_n)$, respectively. Let $f^n_{i, j} = f(x_i, y_j, t_n)$, $g_{i, j} = g(x_i, y_j)$, $(\mu_m)^n_j = \mu_m(y_j, t_n)$ ($m = 1, 2$), $(\mu_3)^n_i = \mu_3(x_i, t_n)$ and $(\widetilde{\mu}_4)^n_i = \widetilde{\mu}_4(x_i, t_n)$.

Then, (2.4) is approximated by the following difference equations:

Also, we obtain the difference equations of (1.1)

Let $\alpha^n_{i, 0}$ be the local truncature error of (2.6). When $u\in C^4(\overline{Q}_T)$, using the Taylor formula, we can easily deduce that

Similarly, when $u\in C^4(\overline{Q}_T)$, it holds that

where $\alpha^n_{i, j}$​ is the local truncation error of (2.7).

Moreover, we obtain

From the above, we obtain the backward Euler difference scheme of problems (1.1)–(1.6).

3.   Error estimate

Let $e^n_{i, j} = u^n_{i, j}-U^n_{i, j}$ be the error of the approximation solution $U$ at the grid point $(x_i, y_j, t_n)$, and $\mu = \frac{\tau}{h^2}$ be the grid ratio. Then, the error equations of (2.9a)–(2.9f) are

Given the complexity of the above error equations, the key to obtaining an error estimate lies in finding transformations that separate the index variables $i$, $j$, and $n$.

Since the error sequence $\{{e}^n_{i, j}\}$ satisfies (3.1c) and (3.1d), applying the DFT to $\{e^n_{i, j}\}$ with respect to $i$, we obtain

Similarly, applying the DFT to $\{{\alpha}^n_{i, j}\}$ with respect to $i$, we obtain

It follows from (2.8) and (3.3) that

Substituting (3.2) and (3.3) into (3.1a), we obtain

Utilizing the properties of the DFT and (3.5), we obtain

Similarly, substituting (3.2) and (3.3) into (3.1f), we have

Substituting (3.2) into (3.1b) and (3.1e), we deduce that

and

Given that the sequence $\{\widehat{e}^n_{k, j}\}$ adheres to the condition specified in (3.9), the conventional DFT is found to be inadequate for our analytical needs. In pursuit of a suitable tool for analysis, we aspire for a novel transformation that not only fulfills the criteria outlined in (3.8) but also possesses the property of invertibility. Drawing inspiration from the formulation of the DFT, we introduce a fresh transformation tailored specifically for the sequence $\{\widehat{e}^n_{k, j}\}$ with respect to $j$ in the following way, aiming to address the aforementioned limitations and meet our analytical needs.

where

It is straightforward to verify Lemma 3.1.

Lemma 3.1. The sequence $\{T_l(y_j)\}$ has the following properties.

(1) $T_l(y_0) = \begin{cases} 1, & l = 0, 1, \cdots, N, \\ 0, & l = N+1, N+2, \cdots, 2N-1. \end{cases}$

(2) $T_l(y_1)-T_l(y_0) = \begin{cases} (\cos{(2l\pi h)}-1)T_l(y_0), & l = 0, 1, \cdots, N, \\ h\sin{(2l\pi h)}, & l = N+1, N+2, \cdots, 2N-1. \end{cases}$

(3) For any $0\leq l, j\leq 2N-1$,

(4) $T_l(y_{2N-j}) = \begin{cases} \cos{(2l\pi y_j)}, & l = 0, 1, \cdots, N, \\ (y_j-1)\sin{(2l\pi y_j)}, & l = N+1, N+2, \cdots, 2N-1. \end{cases}$

For the sake of simplicity in the subsequent analysis, we introduce

and consider the orthogonality relation of the polynomials $P_l(y_i)$ and $P_l(y_j)$.

Lemma 3.2. Given that $i, j = 0, \cdots, N$, we have the following identity:

where

Proof. Using (3.14) and (3.12), and noticing $2Nh = 1$ and $y_i = ih$, we have

For $0\leq m\leq 2N$, we have

i.e.,

If $\cos{(2l\pi y_m)}\neq 1$, then

Now we focus on the case $i\neq j$. Since $0\leq i, j\leq N$, it follows that $0 < y_{i+j} < 1$ and $0 < \left\vert{y_{i-j}}\right\vert < 1$. Furthermore, when $l$ ranges from $1$ to $N−1$, we observe that $\cos{(2l\pi y_{i-j})} \neq 1$ and $\cos{(2l\pi y_{i+j})}\neq 1$.

Thus, with (3.16), and observing the same parity of $i+j$ and $i-j$, we obtain

Substituting the above equality into (3.15), we deduce that

In the next, we consider the case $i = j$. Noting that when $1\leq l\leq N-1$, $\cos{(2l\pi y_{2i})}$ is equal to 1 only in $i = 0$ and $i = N$. Therefore, by utilizing (3.16), we obtain

Substituting (3.18) into (3.15), we deduce that

Furthermore, with (3.14), (3.17), and (3.19), we arrive at the conclusion stated in (3.13). □

Similar to Lemma 3.2, we can derive the subsequent lemma as well.

Lemma 3.3. Given that $i, j = N+1, \cdots, 2N-1$, we have the following identity

Therefore, we can conclude the following lemma.

Lemma 3.4. Suppose

Then

Proof. Using (3.21), (3.12), and Lemma 3.3, we obtain

The proof is finished. □

Similar to Lemma 3.4, we obtain

Lemma 3.5. Suppose

Then

Based on Lemmas 3.4 and 3.5, we obtain the invertible transformation of (3.25).

Lemma 3.6. Suppose

Then,

where

Proof. Using (3.25), (3.12), and Lemma 3.1, we have

and

From the two equalities above, it follows that

and

Moreover, using Lemmas 3.5 and 3.4, we arrive at the conclusion stated in (3.26). □

Similar to (3.10), we use the same transformation to $\{\widehat{\alpha}^n_{k, j}\}$ with respect to $j$ in the following way:

Using (3.4), (3.27), and Lemma 3.6, and noting that $0\leq y_j\leq 1$ $(j = 0, 1, \cdots, 2N)$, we can deduce

Substituting (3.10) and (3.28) into (3.6), we obtain

Using Lemma 1, (3.30) can be rewritten as

Let $l: = 2N-l$ in $\sum\limits^{2N-1}_{l = N+1} \widetilde{e}^{n}_{k, l}T_{2N-l}(y_j)\sin{(2l\pi h)}$, the above equalities have the following form:

Substitute (3.10) into (3.7), then

Moreover, using Lemma 3.1, we obtain

Through comparing with the above equality and (3.31), we find that (3.31) also holds for $j = 0$. Therefore,

Using Lemma 3.6 to perform an invertible transformation on (3.32), we obtain

Let

Obviously,

Using (3.34), (3.33) can be rewritten as

Substituting (3.10) into (3.8), and using Lemma 3.6, we can easily deduce that

Using (3.36) and (3.37), we obtain the following recursive formula for $\left\{\widetilde{e}^n_{k, l}\right\}$:

In order to estimate $\left\{\widetilde{e}^n_{k, l}\right\}$, we first prove the following estimation

In fact, from (3.34) and (3.35), we can derive that

For $0\leq l\leq N$, we have $l\pi h\in [0, \frac{\pi}{2}]$. Observe that $\frac{k\pi h}{2}\in (0, \frac{\pi}{2}) \; (1\leq k\leq 2N-1)$. Therefore, (3.40) can be rewritten as

For $N\leq l\leq 2N-1$, from $2Nh = 1$ and $0 < (2N-l)\pi h\leq \frac{\pi}{2}$, we have

Moreover, (3.40) can be written as

Therefore, combining (3.41) with (3.42), (3.39) holds.

Now, we give the estimation of $\left\{\widetilde{e}^n_{k, l}\right\}$ in three cases.

Case 1. $N\leq l\leq 2N-1$

With (3.38), (3.35), (3.29), and (3.39), we have

Case 2. $1\leq l\leq N-1$

Using (3.43), we obtain

From the above inequality, and using (3.38), (3.35), (3.39), (3.29), and $\mu = \frac{\tau}{h^2}$, we obtain

Case 3. $l = 0$

Observing that $\omega_{k, 0} = \frac{1}{1+4a^2\mu\sin^2{\frac{k\pi h}{2}}}$, similar to deduce (3.43), we obtain

From (3.10), (3.43), (3.45), and (3.46), and noticing that $T_l(y_j)$ is bounded for any $0\leq l, j\leq 2N-1$, we have

Furthermore, from (3.2), we obtain

Since $\frac{1}{x^2+y^2}$ is increasing monotonically with respect to variables $x$ and $y$ for $x, y > 0$, respectively, it follows that

where $\Omega_h = [h, 1]\times [h, 1]$.

Using (3.49) and (3.48), and noticing (3.1c)–(3.1e), we can obtain the following error estimation theorem.

Theorem 3.1. Suppose $u\in C^{4}(\overline{Q}_T)$. For any postive integer $1\leq n \leq M$, the following estimates for (2.9a)–(2.9f)

hold.

4.   Superconvergence anaylysis

In this section, we present the approximation formulas for the partial derivatives of $u$ with respect to two spatial variables, which exhibit superconvergence under certain smooth conditions.

Let $U_x$ and $U_y$ be the approximation functions for the partial derivatives $u_x$ and $u_y$, respectively. For any $t_n$ ($1\leq n\leq M$), we introduce the following approximation formulas for $u_x$ and $u_y$ at the grid point $(x_i, y_j, t_n)$, respectively:

and

Before exploring the superconvergence of (4.1) and (4.2), we first present the following lemma.

Lemma 4.1. Suppose that the function $p(x)\in C^1[0, 1]$ satisfies

where $M$ is a positive constant. If

then

Proof. Let $\theta_k = \frac{k\pi h}{2}$. Obviously, $\theta_k\in (0, \frac{\pi}{2})$. We can easily verify the following equality:

Noting that

the following equality is also verified:

Subtracting (4.7) from (4.8), we have

Thus, using (4.3), we obtain

Using (4.6), (4.3), and (4.9), and noting $\theta_k\in (0, \frac{\pi}{2})$, we obtain

Therefore, the lemma is proved with (4.4). □

Next, we study the superconvergence of (4.1).

Theorem 4.1. Suppose $u\in C^{5}(\overline{Q}_T)$. Then, for any integer $1\leq n\leq M$,

Proof. From (4.1) and $u\in C^{5}(\overline{Q}_T)$, we obtain

Thus, in order to prove this theorem, it suffices to prove the following inequality, i.e., for any given integer $1\leq n\leq M$,

From (3.2) and (3.10), we have

Since $u\in C^5(\overline{Q}_T)$, using Lemma 4.1, (2.8), and (3.3), we obtain

Correspondingly, (3.29) is written as

For this, by modifying (3.43), (3.45), and (3.46), we obtain

Using (4.12) and (4.15), and given $\left\vert{T_l(y)}\right\vert\leq 1$ for $y\in [0, 1]$, we obtain

From this, using (3.49), we prove (4.11). Therefore, (4.10) holds. □

In the following, we discuss the superconvergence properties of (4.2).

Theorem 4.2. Suppose $u\in C^{5}(\overline{Q}_T)$. Then, for any integer $1\leq n\leq M$,

hold.

Proof. From (4.2) and $u\in C^{5}(\overline{Q}_T)$, we obtain

From (4.18), in order to prove this theorem, we only need to prove that for any integer $1\leq n\leq M$,

Using (3.2) and (3.10), and noting that $T_0(y_{j+1})-T_0(y_{j-1}) = 0$, we find that

Using (3.11), when $N\leq l\leq 2N-1$, we have

Furthermore, we also deduce that

Using (4.20), (4.15), (4.21), and (4.22), we obtain

Upon observing

by substituting this result into (4.23), we obtain

This result confirms (4.19), thereby establishing the validity of (4.17). □

5.   Numerical experiments

In this section, we present two numerical examples to validate the theoretical results and investigate the efficiency and the superconvergence properties of the numerical schemes. Our aim is to demonstrate the practical implications of the theoretical findings and assess the performance of the proposed methods. Let

Example 5.1. In (1.1)–(1.6), take

The exact solution is $u = e^{x+y+2t}$ which can be easily verified.

The results are reported in Tables 1–3. From Table 1, we can observe that in the cases of $\tau = h^2$ and $\tau = h$, the error $\|U-u\|_{\infty}$​ is approximately of the order $O(h^2)$ and $O(h)$, respectively. This observation verifies the correctness of Theorem 3.1.

Furthermore, from Tables 2 and 3, it is evident that when $\tau = h^2$, both $\|U_x-u_x\|_{\infty}$ and $\|U_y-u_y\|_{\infty}$ are close to the order $O(h^2)$. On the other hand, when $\tau = h$, $\|U_x-u_x\|_{\infty}$ and $\|U_y-u_y\|_{\infty}$ approach the order $O(h)$. These findings support the theoretical expectations regarding the convergence rates of the spatial derivatives. Therefore, the correctness of Theorems 4.1 and 4.2 is verified.

Example 5.2. In problems (1.1)–(1.6), take

It is easily verified that its exact solution is $u = (1+y)e^{x+t}$.

Numerical results for Example 5.2 are reported in Tables 4–6. These results verify the correctness of Theorems 3.1, 4.1 and 4.2 again.

6.   Conclusions

This work focuses on a heat conduction problem with nonlocal boundary conditions. We develop an implicit Euler scheme and demonstrate that it achieves asymptotic optimal order with the DFT. Furthermore, we introduce two approximation formulas that exhibit superapproximation for first-order partial derivatives along the $x$ and $y$ directions of the exact solution, respectively. In the future, we plan to extend this work to other difference schemes for parabolic problems with nonlocal boundary conditions, such as the explicit Euler scheme, the Crank-Nicolson scheme, and other schemes. Additionally, we aim to consider heat conduction problems with different nonlocal boundary conditions.

Author contributions

Liping Zhou: Conceptualization, methodology, formal analysis, writing-original draft, validation; Yumei Yan: Editing, software; Ying Liu: Writing-review and editing. All authors contributed equally to the manuscript. All authors have read and approved the final version of the manuscript for publication.

Acknowledgments

This work is partially supported by the National Natural Science Foundation of China (No.12101224), the Hunan Provincial Natural Science Foundation of China (No. 2022JJ30271, 2024JJ7203) and the Key Project of Hunan Provincial Education Department of China (No. 23A0577).

Conflict of interest

The authors declare no conflicts of interest.