On the impact of quantum biology and relativistic time dilation in autism

1.   Introduction

Let $ \Omega = $ $ [b_{1}, b_{2}] $ $ \times [d_{1}, d_{2}] $, where $ b_{1} $, $ b_{2} $, $ d_{1} $ and $ d_{2} $ are constant numbers and belong to $ \mathbf{R} $. Then we denote $ \Omega = (b_{1}, b_{2})\times(d_{1}, d_{2}) $ $ \subset \mathbf{R}^{2} $ and $ {\mathbf x} = (x, y) $. In this study, we are concerned with the numerical solutions of the following nonlinear wave equations with time delay ^{[1,2,3,4,5,6,7,8]}

by devising efficient numerical methods. Here, $ s > 0 $ is a constant delay, $ x $ and $ y $ represent spatial variables in $ x $- and $ y $-directions, respectively, and $ t $ denotes temporal variable. Denote $ (x, y) $ by $ {\mathbf x} $ for brevity. The exact solution to the problem (1.1a)–(1.1c) is denoted by $ u({\mathbf x}, t) $. Furthermore, we assume that $ f(u({\mathbf x}, t), u({\mathbf x}, t-s), {\mathbf x}, t) $, $ \phi({\mathbf x}, t) $ and $ \varphi({\mathbf x}, t) $ are sufficiently smooth to make the regularity of the exact solutions satisfied, and our algorithms achieve claimed accuracies.

Delayed partial differential equations (DPDEs) are extensively utilized to simulate many natural, social and engineering phenomena such as population ecology, cell biology and control theory ^[9,10]. For example, the following nonlinear delayed reaction-diffusion equations

is used to describe many practical problems. Also, in Eqs (1.2a)–(1.2c), $ s > 0 $ is a time delay, $ u(x, t) $ denotes the exact solution to the IBVPs (1.2a)–(1.2c), and $ x $ and $ t $ represent spatial and temporal variables, respectively. As

Eq (1.2a) is referred to as Hematopoiesis model, and applied to describe the dynamics of blood cell production (cf. ^[11,12,13]). Denote the density of mature stem cells by $ u(x, t) $ in blood circulation. The time delay between the production of immature stem cells in the bone marrow and their maturation for release in the circulating blood stream is represented by $ s \ge 0 $. Parameters $ \alpha $, $ \beta $, $ \theta $, $ \gamma $ and $ d\in (1, \infty) $ are positives constants. Besides, as

Eq (1.2a) is known as diffusive Nicholson's blowflies equation (cf. ^[14,15,16]). Here, $ u(x, t) $, $ \varrho $, $ 1/\xi $, $ \nu $ and $ s $ denote the size of the population at time $ t $, represents maximum per capita daily egg production rate, the size at which the population reproduces at its maximum rate, the per capita daily adult death rate and the generation time, respectively. Finally, as

Eq (1.2a) is well-known Hutchinson equation. It is proposed by American biologist Hutchinson in 1948 (cf. ^[17,18,19]). Here, $ u(x, t) $ denotes the current population density of mammals. For much more detailed mathematical models with delays, readers are referred to ^[9,10]. Hence, it is important and necessary to comprehensively research this kinds of equations from the theoretical point of view and from numerical analysis because they possess strong applied background. In these three examples, $ u(x, t) $ should be more than zero according to its meaning.

Theoretical studies on the exact solutions of the DPDEs, such as the unique existence, asymptotic behavior and blow up, have been attracting a lot of attention. Considerably theoretical findings can be found in ^[9,10] and the related references. However, for DPDEs with arbitrary initial-boundary conditions, it is difficult to obtain the exact solutions of them. Thus, it is of great significance to explore high-performance numerical algorithms used to solve them.

In the past two decades, there have been great advance in numerical researches for delayed parabolic equations (DPEs). For example, a FDM combined with domain decomposion technique was derived for solving one-dimensional (1D) nonlinear delayed DPEs in ^[20]. It is second-order accurate in both time and space. To improve spatial accuracy, a compact FDM, which is second-order and fourth-order accurate in time and space, respectively, was developed for 1D and two-dimensional (2D) DPEs in ^[21]. In ^[22], it is shown that the use of fourth-order compact FDM to the spatial discretization leads to a reduction of the asymptotic stability region of the original DPEs. Recently, a class of $ \theta $-schemes combined with Lagrange interpolation has been suggested for nonlinear parabolic equations with variable delays in ^[23]. A Crank-Nicolson alternating directional implicit (ADI) FDM ^[24], which is second-order accurate in both time and spaces, was first developed for 2D DPEs in ^[24]. A Crank-Nicolson ADI compact FDM, which is second-order and fourth-order accurate in time and spaces, respectively, was developed in ^[25]. By using second-order backward differentiation formula (BDF2) to discrete temporal derivatives, a multistep ADI compact FDM, which is second-order and fourth-order accurate in time and spaces, respectively, was developed for 2D nonlinear delayed parabolic equations with constant coefficients in ^[26]. Similarly, two multistep ADI compact FDMs, which are second-order and fourth-order accurate in time and spaces, respectively, were developed for 2D nonlinear delayed parabolic equations with different variable coefficients in ^[27] and ^[28], respectively. Other efficient numerical methods including locally discontinuous Galerkin method ^[29], waveform relaxation methods combined with spectral collocation methods ^[30] and finite volume element method ^[31] were developed for solving 1D parabolic equations with constant delays. Furthermore, a box scheme with second-order temporal and spatial accuracies was derived for delayed convection-diffusion equations with Riemannian boundary conditions in ^[32]. A compact FDM, which is second-order and fourth-order accurate in time and space, respectively, was developed for delayed convection-diffusion equations in ^[33]. Moreover, Legendre spectral Galerkin method was also established for solving 2D convection-diffusion equations with delays in ^[34]. Similarly, FDMs ^[35,36], multi-domain Legendre spectral collocation method ^[37] and Galerkin methods ^[38] were constructed for the neutral DPEs. More recently, in ^[39], a compact difference method with second-order temporal accuracy and fourth-order spatial accuracy, which was used along with a Richardson extrapolation method (REM) to obtain fourth-order temporally and spatially accurate solutions, was devised for Sobolev equations with constant delay. A adaptive FDM based on Bakhvalov mesh has been constructed for a singularly perturbed Sobolev equations with constant delay in ^[40]. An exponentially-fitted difference scheme on a uniform mesh was designed for a singularly perturbed Sobolev equations with constant delay in ^[41].

Also, much attention has been paid on the theoretical studies of the delayed wave equations. A number of theoretical results, such as the well posedness, stability, long time behavior and existence of the exact solutions, can be found in ^[1,2,3,4,5]. Moreover, separation of variables and Lie group method were applied to find the exact solutions of them in ^[6,7] and ^[8], respectively. In contrast, there have been only a few numerical studies on nonlinear wave equations with delays. A three-level compact FDM, which is second-order and fourth-order accurate in time and space, respectively, has been derived for 1D delayed wave equation in ^[42]. Subsequently, a three-level ADI compact FDM, which is also second-order and fourth-order accurate in time and spaces, respectively, has been derived for 2D delayed wave equation in ^[43]. To promote the numerical studies for delay wave equations, two explicit finite difference methods (EFDMs) and corresponding REMs are designed to solve the delayed wave equations.

As we know, the stable requirement of the standard EFDM for linear wave equation is accepted. Optimally convergent solutions can be obtained as long as spatial and temporal meshsizes are reasonably chosen in accordance with the stable requirement. Besides, Du Fort-Frankel scheme proposed by Du Fort and Frankel for 1D linear parabolic equations with periodic boundary conditions is unconditionally Riemann stable, and fast convergent as consistent condition is satisfied (cf. ^[44]). However, there has been no numerical studies for nonlinear wave equation with delays by Du Fort-Frankel-type schemes. In this paper, inspiring by these two types of EFDMs, two EFDMs, which can be used along with REMs to get high-order accurate numerical solutions of nonlinear wave equation with delay, have been designed. By using the discrete energy methods, error estimations are deduced, rigorously. There are some advantages of our algorithms listed as follows. (1) The algorithms are both very easy to be implemented because they are explicit schemes. (2) The explicit methods combined with REMs possess low computational burden because of explicitness and high-order accuracy. (3) They are easily generalized to multi-dimensional problems with various kinds of delays, such as multiple delays, variable delays and distributed delays. What's more, a main contribution of this study is that Du Fort-Frankel scheme and a new REM are constructed for nonlinear wave equation with delay. They possess good stability, high-order accuracy and low computational burden.

The basic organization of the article is as follows. Partition and notations are given in Section 2. The construction and theoretical analyses of the the first EFDM and corresponding REM are given in Section 3. Section 4 is devoted to the establishment and theoretical analyses of the the second EFDM and corresponding REM. The extensions of the proposed methods to wave equation with other delays are discussed in Section 5. Section 6 concentrates on numerical experiments. A conclusion is summarized in Section 7.

2.   Partition and Notations

To begin with, we divide the spatial domain $ \Omega = $ $ [b_{1}, b_{2}] $ $ \times [d_{1}, d_{2}] $ $ \subset $ $ \mathbf{R}^{2} $, where $ b_{1} $, $ b_{2} $, $ d_{1} $ and $ d_{2} $ are constant numbers. Let $ h_x = (b_2-b_1)/m_1 $, $ h_y = (d_2-d_1)/m_2 $, $ (m_1, m_2\in Z^+) $ be spatial meshsizes in $ x $- and $ y $-directions, respectively.

Let $ T $ and $ s $ be positive constants. Then the temporal domain and delay are denoted by $ [0, T] $ and $ s $, respectively. As ^{[20,21,22,23]}, the restricted grid with meshsize $ \tau $ $ (\tau = s/n_1 = T/n) $ $ (n_1, n\in Z^+) $ is used. If the restricted grid is not used, Lagrange interpolation formula will be applied to approximate the delay term $ u({\mathbf x}, t-s) $. The use of Lagrange interpolation formula brings about computational complexity and the failure of improving temporal accuracy by REMs. Thus, for constant time delays, we choose constrained timestep.

By denoting $ t_k = k\tau $, $ h = \max{\{h_{x}, h_{y}\}} $, $ x_{i} = b_1+ih_{x} $, $ y_{j} = d_{1}+jh_{y} $ and $ (x_{i}, y_{j}) $ $ = {\mathbf x}_{i, j} $, the domain $ \Omega\times[0, T] $ is covered by $ \Omega_{h\tau} $ $ = \bar{\Omega}_h\times\Omega_\tau $, where

Besides, we introduce

On $ \Omega_{h\tau} $, further define

Then, on $ \Omega_{h}\times\Omega_{\tau} $, we further introduce the following difference notations

Finally, $ \forall $ $ U $, $ V $ $ \in S^{0}_{h} $, we define the norms as follows

3.   The construction and analysis of the first EFDM

3.1.   The construction of the first EFDM

Define grid functions $ u_{i, j}^k = u({\mathbf x}_{i, j}, t_k) $, $ {\mathbf x}_{i, j}\in\Omega_h $, $ -n_1\leq k\leq n $. The approximation of $ u^{k}_{i, j} $ is denoted by $ U^{k}_{i, j} $.

Applying second-order centered difference formulas to approximate temporal and spatial derivatives at point $ ({\mathbf x}_{i, j}, t_k) $ yields that

where

Omitting the small term $ (R_1)_{i, j}^k $, and then replacing $ u_{i, j}^k $ with $ U_{i, j}^k $ in Eq (3.1), we get the first EFDM as follows

Obviously, the EFDM (3.3a)–(3.3c) is uniquely solvable because it is an explicit scheme.

3.2.   The convergence analysis of the first EFDM

Besides, we give an assumption for $ f(\mu $, $ \upsilon $, $ {\mathbf x}, t) $ as follow.

Assumption A. Let $ u({\mathbf x}, t) $ be the exact solution to the problem (1.2a)–(1.2c). Then, there are positive constants $ c_1 $ and $ \varepsilon_0 $, and $ \forall $ $ \varepsilon_{l} $ $ (l = 1, \; 2) $ satisfying $ |\varepsilon_l| < \varepsilon_0 $, $ (l = 1, 2) $, it holds that

Moreover, denote $ e_{i, j}^k = u_{i, j}^k-U_{i, j}^k $, $ 0\leq i\leq m_1 $, $ 0\leq j\leq m_2 $, $ -n_1\leq k\leq n $ and $ h = \max{\{h_x, h_y\}} $, and give several lemmas used later.

Lemma $ 3.1 $ ^[45] Let $ V\in S_h^{0} $, we have

Lemma $ 3.2 $ ^[26] Let both A and B be positive constants, $ \{F^k|k\geq0\} $ be a non-negative sequence and satisfy the following inequality $ F^{k+1}\leq A+B\tau\sum\limits_{k = 0}^KF^k $. Then we can obtain $ \max\limits_{0\leq k\leq K+1}F^k\leq A\exp(B(K+1)\tau) $. Moreover, if $ F^{k+1}\leq A+B\tau\sum\limits_{k = 0}^{K+1}F^k $, we have $ \max\limits_{0\leq k\leq K+1}F^k\leq A\exp(2B(K+1)\tau) $ as long as $ \tau $ is sufficiently small such that $ B\tau\leq\frac{1}{2} $.

Lemma $ 3.3 $ ^[46] Let $ r_\chi = (|a|\tau)/h_{\chi} $ ($ \chi = x $ or $ y $), $ E^k = ||\delta_te^{k+\frac{1}{2}}||^2+a^2\langle e^k, e^{k+1}\rangle_{1, x}+a^2\langle e^k, e^{k+1}\rangle_{1, y} $, $ c_2 = 1-r_x^2-r_y^2 $. Then when $ r_x^2+r_y^2 < 1 $, the following inequalities

hold.

Proof. By using the equality $ \alpha\beta = [(\alpha+\beta)^{2}-(\alpha-\beta)^{2}]/4 $ and simple computations, we have

which is used along with $ r_x^2+r_y^2 < 1 $ to yield Eq (3.4a).

According to $ \delta_xe_{i+\frac{1}{2}, j}^{k+1} = \frac{\tau}{2} \delta_x\delta_te_{i+\frac{1}{2}, j}^{k+\frac{1}{2}}+ \delta_xe_{i+\frac{1}{2}, j}^{k+\frac{1}{2}} $, by simple computation, it is easy to find that

Multiplying $ a^2h_xh_y $ to both sides of Eq (3.5) and summing $ i $ from $ 1 $ to $ m_{1}-1 $ and $ j $ from $ 1 $ to $ m_{2}-1 $, we obtain

which shows that

Also, noting $ \delta_ye_{i, j+\frac{1}{2}}^{k+1} = \frac{\tau}{2} \delta_y\delta_te_{i, j+\frac{1}{2}}^{k+\frac{1}{2}}+ \delta_ye_{i, j+\frac{1}{2}}^{k+\frac{1}{2}} $ and utilizing the analytical methods similar to Eq (3.6) get that

Adding Eq (3.6) to Eq (3.7) and using Eq (3.4a) lead to

which is used along with Lemma 3.1 to find that

The proof is finished.

Let exact solution $ u({\mathbf x}, t)\in C^{4, 4}(\Omega\times[-s, T]) $. Then there exists positive constant $ c_3 $ such that

Theorem 3.1 Let $ u({\mathbf x}, t) $ $ \in \mathcal {C}^{4, 4}(\Omega\times [-s, T]) $ be the exact solution of (1.2a)–(1.2c), and $ U_{i, j}^k $ be the solution of the scheme (3.3a)–(3.3c) at the time level $ k $. Assuming that there exist positive constants $ \mu $, $ \sigma $ and $ \epsilon $, such that $ \mu h\leq h_x $ $ \le h $, $ \mu h \le $ $ h_y \leq h $ and $ \tau\leq \sigma h^{\frac{1}{2}+\epsilon} $, we can derive that

under the conditions that

and Assumption A are valid. Here,

Proof. For convenience, denote $ f(u_{i, j}^k) = f(u_{i, j}^k, u_{i, j}^{k-n_1}, x_{i}, y_{j}, t_{k}), f(U_{i, j}^k) = f(U_{i, j}^k, U_{i, j}^{k-n_1}, x_{i}, y_{j}, t_{k}) $ without confusion. Subtracting Eq (3.3a) from Eq (3.1), the error equations are derived as follows

As $ e_{i, j}^k = 0 $ ($ -n_1\leq k\leq0 $), it is obvious that Eq (3.10) holds for $ -n_1\leq k\leq 0 $. Assuming that Eq (3.10) is valid for $ -n_{1}\leq k\leq l $, we will prove that it is also true for $ k = l+1 $.

As $ \mu h \le h_{x} \le h $ and $ \mu h \le h_{y} \le h $, it is easy to find that

Thus, using $ \tau\le \sigma h^{\frac{1}{2}+\epsilon} $ and Eq (3.12), we obtain

Besides, applying $ h\leq\min\{((\frac{\varepsilon_0\mu}{2\sigma^2c_4}) $ $ ^{\frac{1}{2\epsilon}} $, $ \frac{\varepsilon_0\sqrt{\mu}}{4c_4})\} $, it is easy to find that

As a result, using the induction hypothesis and Eqs (3.12)–(3.14) yields that

which follows from Assumption A that

Taking inner product with Eq (3.11a) by $ 2\delta_{\hat{t}}e^k $ gives that

In what follows, we estimate each terms of Eq (3.16) step by step. Using the equality $ \alpha^{2}-\beta^2 $ $ = (\alpha-\beta)(\alpha+\beta) $ gives that

Applying the discrete Green formula yields that

Similar to Eq (3.18), we also arrive at

Next, we estimate each terms at the right of Eq (3.16). Applying Eq (3.15), the inequalities $ 2\alpha\beta\leq\frac{1}{\varepsilon}\alpha^2+\varepsilon \beta^2 $ and $ (\alpha+\beta)^2\leq 2(\alpha^2+\beta^2) $, and Lemma 3.3, we have

Using the analytical methods similar to Eq (3.20) yields that

Substituting Eqs (3.17)–(3.21) into Eq (3.16) infers that

Using Eq (3.4c) in Eq (3.22), replacing $ k $ with $ \gamma $ and summing $ \gamma $ from $ 0 $ to $ k $, we have

Hence, applying Lemma 3.2 to Eq (3.23), it follows from $ \tau\{4+\frac{8c_3^2(b_2-b_1)^2(d_2-d_1)^2} {3c_2^2a^2[(b_2-b_1)^2+(d_2-d_1)^2]}\}\leq 1 $ that

In Eq (3.24), setting $ k = l $, applying Lemma 3.1, Eqs (3.4b)–(3.4c) in Lemma 3.3, we can conclude that Eq (3.10) holds for $ k = l+1 $.

Thus, by mathematical induction, this theorem is valid.

3.3.   The improvement of the computational efficiency

This section is suggested for the improvement of the computational efficiency of the EFDM in Eqs (3.3a)–(3.3c) by using REMs. As this EFDM (3.3a)–(3.3c) is conditionally stable, temporal meshsize is often taken smaller than spatial meshsizes. Thus, a REM in spaces is devised to improve the computational efficiency. In this section, $ C $ is a positive constant and independent of grid parameters $ \tau $, $ h_{x} $ and $ h_{y} $.

Theorem 3.2 Suppose that $ u({\mathbf x}, t) $ $ \in \mathcal {C}^{6, 4}(\Omega\times [-s, T]) $ is the exact solution of Eqs (1.2a)–(1.2c), and $ U^{k}_{i, j}(\tau, h_{x}, h_{y}) $ is the numerical solution of the 1st EFDM (3.3a)–(3.3c) at $ ({\mathbf x}_{i, j}, t_{k}) $ using meshsizes $ h_{x} $, $ h_{y} $ and $ \tau $. Define a REM as follows

Then, as $ 4a^2[(\frac{\tau}{h_{x}})^2+(\frac{\tau}{h_{y}})^2] < 1 $, and grid parameters $ \tau $, $ h_{x} $ and $ h_{y} $ are small enough, we obtain

Proof. Let $ r_1({\mathbf x}, t) = \frac{a^2}{12}\frac{\partial^4u}{\partial x^4}({\mathbf x}, t) $ and $ r_2({\mathbf x}, t) = \frac{a^2}{12}\frac{\partial^4u}{\partial y^4}({\mathbf x}, t) $. Then a combination of Taylor expansion formula with (3.2) gives that

By Eq (3.27), we can rewrite the error equations Eqs (3.11a)–(3.11c) as

Furthermore, two auxiliary problems are introduced to derive the asymptotic expansion of the numerical solutions. Namely, let $ v_{1}({\mathbf x}, t) $ and $ v_{2}({\mathbf x}, t) $ be the exact solutions of the following initial-boundary-value problems (IBVPs)

respectively.

Denote $ (v_{1})_{i, j}^{k} $ $ = v_{1}({\mathbf x}_{i, j}, t_{k}) $ and $ (v_{2})_{i, j}^{k} $ $ = v_{2}({\mathbf x}_{i, j}, t_{k}) $. Applying the difference methods similar to those utilized in Eqs (3.3a)–(3.3c) to discrete Eqs (3.29b)–(3.29d) and Eqs (3.30b)–(3.30d) yields that

respectively.

Denote $ p_{i, j}^{k} = e_{i, j}^k+h_{x}^{2}(v_{1})_{i, j}^{k}+h_{y}^{2}(v_{2})_{i, j}^{k} $, $ 0\leq i\leq m_1 $, $ 0\leq j\leq m_{2} $, $ -n_{1}\leq k\leq n $. Then, multiplying Eqs (3.31b)–(3.31d) by $ h_x^2 $ and Eqs (3.32b)–(3.32d) by $ h_{y}^{2} $, respectively, and then adding the obtained results to Eqs (3.28b)–(3.28d), we have

where

and

In Eq (3.33a), the following Taylor expansion formula

is used.

Then, as $ a^2[(\frac{\tau}{h_{x}})^2+(\frac{\tau}{h_{y}})^2] < 1 $, using the same techniques as those proposed in the proof of Theorem $ 3.1 $, we find that

where, $ C $ is a positive constant number and independent of grid parameters. Namely,

On the other hand, from Richardson extrapolation solutions

it follows that

Using the triangle inequality to Eq (3.38) and noting Eq (3.36) yields that

holds as long as $ 4a^2[(\frac{\tau}{h_{x}})^2+(\frac{\tau}{h_{y}})^2] < 1 $. The proof is completed.

Remark I Theorem 3.1 shows that the first EFDM (3.3a)–(3.3c) is conditionally convergent with an order of $ O(\tau^{2}+h^{2}_{x}+h^{2}_{y}) $ as $ r^{2}_{x}+r^{2}_{y} < 1 $. Besides, Theorem 4.1 shows that REM (3.25) is also conditionally convergent with an order of $ O(\tau^{2}+h^{2}_{4}+h^{4}_{y}) $ as $ 4(r^{2}_{x}+r^{2}_{y}) < 1 $ (see Eq (3.39)). Besides, as $ h_{x} = h_{y} = h $ and $ \tau = h^{2} $, Richardson extrapolation solutions have a convergent rate of $ O(h^{4}) $. In Section 6, to obtain numerical solutions with optimal convergent rates, we take temporal and spatial grids according to these conditions.

4.   The development and analysis of the second EFDM

In this section, the famous Du Fort-Frankel scheme is generalized to the numerical solutions of delay nonlinear wave Eqs (1.2a)–(1.2c).

4.1.   The establishment of the second EFDM

Using the following difference formulas

to approximate the spatial derivatives of Eqs (1.2a)–(1.2c) yields that

where

Omitting the small term $ (R_2)_{i, j}^{k} $ and replacing $ u^{k}_{i, j} $ by its approximation $ U^{k}_{i, j} $ in Eq (4.2), we get the second EFDM as follows

4.2.   The convergence analysis of the second EFDM

This subsection is suggested for the convergence of the second EFDM (4.4a)–(4.4c).

Lemma $ 4.1 $ Let $ G^k = (1+r_x^2+r_y^2)\|\delta_te^{k+\frac{1}{2}}\|^2+a^{2}\langle e^k, e^{k+1}\rangle_{1, x}+a^{2}\langle e^k, e^{k+1}\rangle_{1, y} $. Then the following inequalities

are valid.

Proof. Using Lemma 3.1, we obtain

Applying Eqs (4.6a), (4.6b) and the equality $ \alpha\beta = $ $ (\alpha+\beta)^{2}/4 $ $ -(\alpha-\beta)^{2}/4 $, we have that

which shows that Eq (4.5a) is valid.

Using the skills similar to those developed in Eqs (3.6) and (3.7), Lemma $ 3.1 $ and Eq (4.5a) in Lemma $ 4.1 $ gives

Making use of Eqs (4.7a) and (4.7c) directly infers that

The proof is completed.

Assuming $ u({\mathbf x}, t) $ $ \in C^{4, 4}(\Omega\times[-s, T]) $, there exists positive constant $ c_6 $ such that

Theorem 4.1 Let $ u({\mathbf x}, t) $ $ \in \mathcal {C}^{4, 4}(\Omega\times [-s, T]) $ be exact solution of Eqs (1.2a)–(1.2c), $ U_{i, j}^{k} $ be the solution of the scheme Eqs (4.4a)–(4.4c) at node $ ({\mathbf x}_{i, j}, t_{k}) $. Assume that there are constants $ \mu $, $ \sigma $ and $ \epsilon $ such that $ \mu h\leq h_{x}\le h $, $ \mu h\leq h_{y}\le h $, $ \tau\leq\sigma h^{\frac{3}{2}+\epsilon} $. It holds that

under the conditions that

and Assumption A hold. Here,

Proof. Also, let $ f(u_{i, j}^{k}) = f(u_{i, j}^{k}, u_{i, j}^{k-n_{1}}, {\mathbf x}_{i, j}, t_{k}) $, $ f(U_{i, j}^{k}) = f(U_{i, j}^{k}, U_{i, j}^{k-n_{1}}, {\mathbf x}_{i, j}, t_{k}) $ for convenience. Subtracting Eq (4.4a) from Eq (4.2), the error equations are derived as follows

As $ e_{i, j}^k = 0 $ ($ -n_1\leq k\leq0 $), it is obvious that Eq (4.9) holds for $ -n_1\leq k\leq0 $. Assuming that Eq (4.9) is valid for $ -n_1\leq k\leq l $, we will prove that it is also true for $ k = l+1 $.

As $ \mu h \le h_{x} \le h $ and $ \mu h \le h_{y} \le h $, it is easy to find that

Thus, using $ \tau\le \sigma h^{\frac{3}{2}+\epsilon} $ and Eq (4.11), we obtain

Besides, applying $ h\leq\min{((\frac{\mu\varepsilon_0} {3c_7\sigma^2})^{\frac{1}{2(1+\epsilon)}}, \frac{\sqrt{\mu}\varepsilon_0} {6c_7}, (\frac{\mu^3\varepsilon_0}{6c_7\sigma^2})^{\frac{1}{2\epsilon}})} $, it is easy to find that

Thus, by applying the induction hypothesis, Eq (4.12) and Eq (4.13), we can deduce that

Hence, applying Assumption A directly infers that

Acting the inner product with Eq (4.10a) by $ 2\delta_{\hat{t}}e^k $ yields that

Similar to the proof of the Theorem $ 3.1 $, and using Lemma $ 4.1 $, the formula at the left of Eq (4.16) is estimated as follows

In what follows, we estimate each terms at the right of Eq (4.16). Applying the inequalities $ 2\alpha\beta\leq\frac{1}{\varepsilon}\alpha^2+\varepsilon \beta^2 $ and $ (\alpha+\beta)^2\leq2(\alpha^2+\beta^2) $, Eq (4.15) and Lemma $ 4.1 $, we can get that

Using Cauchy-Schwarz inequality and Lemma 4.1, it is easy to find that

Inserting Eq (4.17), Eq (4.18) and Eq (4.19) into Eq (4.16), adopting Lemma 3.1 and Eq (4.5b) in Lemma 4.1, and replacing $ k $ with $ \gamma $, summing $ \gamma $ from $ 0 $ to $ k $, we can arrive at

As $ \tau\{4+\frac{4c_1^2(b_2-b_1)^2(d_2-d_1)^2(1+r_x^2+r_y^2)} {3[(b_2-b_1)^2+(d_2-d_1)^2]a^2}\}\leq 1 $, applying Lemma 3.2 to Eq (4.20) derives that

Taking $ k = l $ in Eq (4.21) and using Eq (4.5c) derives that Eq (4.9) holds for $ k = l+1 $. From induction hypothesis, it follows that the claimed finding holds. The proof is completed.

Remark II From Theorem 4.1, we can find that numerical solutions obtained by the second EFDM (4.4a)–(4.4c) converge to exact solution with an order of $ O(\tau^{2}+h_{x}^{2}+h_{y}^2+\frac{\tau^2}{h_{x}^{2}}+\frac{\tau^2}{h_{y}^{2}}) $ in $ H^{1} $-norm without any restriction no $ r_{x} $ and $ r_{y} $. From truncation error Eq (4.3), it follows that the second EFDM (4.4a)–(4.4c) is conditionally consistent, and $ (R_{2})^{k}_{i, j} $ tends to zero as $ \tau = o(h_{x}) $, $ \tau = o(h_{y}) $, and $ h_{x} $ and $ h_{y} $ tend to zero. What's more, although the condition $ \tau\le\sigma h^{\frac{3}{2}+\epsilon} $ ($ \sigma $ and $ \epsilon $ are both positive constant.) is necessary in Theorem 4.1, we would like to adopt the grid $ \tau = h^{2} $ and $ h_{x} = h_{y} = h $ to get numerical solutions with an optimal convergent order of $ O(h^{2}) $ in practical computations.

4.3.   Richardson extrapolation method

In what follows, a REM in both space and time is developed for the 2nd EFDM (4.4a)–(4.4c).

Theorem 4.2 Assume that $ u({\mathbf x}, t) $ $ \in \mathcal {C}^{4, 4}(\Omega\times [-s, T]) $ is the exact solution of Eqs (1.2a)–(1.2c), and $ U^{k}_{i, j}(\tau, h_{x}, h_{y}) $ is the numerical solution of the second EFDM (4.4a)–(4.4c) at $ ({\mathbf x}_{i, j}, t_{k}) $ using meshsizes $ h_{x} $, $ h_{y} $ and $ \tau $. Define the following REM

Then under the conditions of Theorem 4.1, $ h_{x} = O(h_{y}) $, $ \frac{\tau}{h_{x}} $ and $ \frac{\tau}{h_{y}} $ $ \longrightarrow 0 $ as grid parameters $ \tau $, $ h_{x} $, $ h_{y} $ $ \longrightarrow 0 $, we arrive at

as long as parameter grids $ \tau $ $ h_{x} $ and $ h_{y} $ are small enough.

Proof. Let $ w_1({\mathbf x}, t) = \frac{1}{12}\frac{\partial^4u}{\partial t^4}({\mathbf x}, t) $, $ w_2({\mathbf x}, t) = \frac{a^2}{12}\frac{\partial^4u}{\partial x^4}({\mathbf x}, t) $, $ w_{3}({\mathbf x}, t) = \frac{a^2}{12}\frac{\partial^{4}u}{\partial y^{4}}({\mathbf x}, t) $ and $ w_4({\mathbf x}, t) = a^2\frac{\partial^2u}{\partial t^2}({\mathbf x}, t) $. Then, we define the grid functions $ (w_l)_{i, j}^k $ $ = w_l({\mathbf x}_{i, j}, t_k) $, $ (l = 1, 2, 3, 4) $, $ {\mathbf x}_{i, j}\in \bar{\Omega}_h $, $ -n_1\leq k\leq n $. Thus, it follows from Eq (4.3) that

Furthermore, the error Eqs (4.10a)–(4.10c) can be rewritten as

Suppose that functions $ v_{1}({\mathbf x}, t) $, $ v_{2}({\mathbf x}, t) $, $ v_{3}({\mathbf x}, t) $ and $ v_{4}({\mathbf x}, t) $ are exact solutions of the following IBVPs

respectively, where

Denote $ (v_{l})_{i, j}^k = v_{l}({\mathbf x}_{i, j}, t_k) $, $ l = 1, 2, 3, 4 $. Approximating IBVPs (4.26a)–(4.26c), IBVPs (4.27a)–(4.27c), IBVPs (4.28a)–(4.28c) and IBVPs (4.29a)–(4.29c) by using the numerical methods similar to Eqs (4.4a)–(4.4c) yields that

Denote $ q_{i, j}^k = e_{i, j}^{k}+h_{x}^{2}(v_{1})_{i, j}^{k}+h_{y}^{2}(v_{2})_{i, j}^{k}+(\frac{\tau^2}{h_x^2}+\frac{\tau^2}{h_y^2})(v_3)_{i, j}^k+\tau^2(v_4)_{i, j}^k $, $ 0\leq i\leq m_1 $, $ 0\leq j\leq m_2 $, $ -n_1\leq k\leq n $. Then, multiplying $ h^{2}_{x} $, $ h^{2}_{y} $, $ (\frac{\tau^2}{h_x^2}+\frac{\tau^2}{h_y^2}) $ and $ \tau^2 $ to both sides of Eqs (4.31b)–(4.31d), Eqs (4.32b)–(4.32d), Eqs (4.33b)–(4.33d) and Eqs (4.34b)–(4.34d), respectively, and then adding up the obtained results, we arrive at

where

It is worth mentioning that the following Taylor expansion formula

is applied in Eqs (4.35a)–(4.35c).

Then, using the same techniques as those proposed in the proof of Theorem $ 4.1 $ and using the inequalities $ \alpha\beta\le (\alpha^2+\beta^2)/2 $ and $ \max((\frac{h_{y}}{h_{x}})^2, (\frac{h_{x}}{h_{y}})^2)\le \mu^{-2} $, we have

where $ C $ is a positive constant and independent of grid parameters $ \tau $, $ h_{x} $ and $ h_{y} $. Finally, similar to Eq (3.39), using the triangle inequality to the following equality

and applying the estimation Eq (4.38) yield the claimed results. The proof is completed.

Remark III Like Remark II, it follows from Theorem 4.2 that numerical solutions obtained by the second EFDM (4.4a), (4.4b) combined with REM in Eq (4.22) converge to the exact solution with an order of $ O(\tau^{2}+\tau^{4}+h_{x}^{4}+h_{y}^{4} +\frac{\tau^{4}}{h_{x}^{2}}+\frac{\tau^{4}}{h_{y}^{2}} +(\frac{\tau}{h_{x}})^{4}+(\frac{\tau}{h_{y}})^{4}) $ in $ H^{1} $-norm. According to the convergence rate $ O(\tau^{2}+\tau^{4}+h_{x}^{4}+h_{y}^{4} +\frac{\tau^{4}}{h_{x}^{2}}+\frac{\tau^{4}}{h_{y}^{2}} +(\frac{\tau}{h_{x}})^{4}+(\frac{\tau}{h_{y}})^{4}) $, the optimal estimation $ \|u^k-(U_E)^k\|_1 = O(h^{4}) $, $ 0\leq k \leq n $ can be obtained as an optimally temporal and spatial relationship $ h_{x} = h_{y} = h $ and $ \tau = h^{2} $ is used. Thus this relationship $ h_{x} = h_{y} = h $ and $ \tau = h^{2} $ is applied in practical computations.

5.   Extension to nonlinear wave equation with several delays

Partial differential equations with several delays are extensively applied in scientific and engineering fields. For example, a delayed convective-diffusion equation (see ^[10])

which arises from furnace control, models a furnace used to process metal sheets. Here, $ u $ is the temperature distribution in a metal sheet, $ v $ represents velocity, $ f $ denotes a distributed temperature source function, both $ v $ and $ f $ are dynamically adapted by a controlling device monitoring the current temperature distribution. However, the finite speed of the controller leads to two fixed delays of length $ \tau_{1} $ and $ \tau_{2} $. Thus it is necessary and important to study numerical solutions of partial differential equations with several delays.

This section focuses on the generalization of the proposed numerical algorithms and corresponding analytical results to 2D nonlinear wave equation with several delays as follows

where $ s_{l}(t) > 0 $ $ (l = 1, 2, \ldots, L) $ and $ s = \max_{1\le l \le L}{s_{l}(t)} $.

As $ s_{\kappa}(t) $ $ = s_{\kappa} $ ($ \kappa = 1, 2, \cdots, L $) are all constant delays, i.e. independent of temporal variable $ t $, constrained temporal grid $ \tau = s_{1}/n_{1} = s_{2}/n_{2} = \cdots = s_{L}/n_{L}, (n_{l} \in \mathbb{Z}^{+}) $ is adopted to make $ t = K s_{l} $ $ (K \in \mathbb{Z}^{+}) $ locate the temporal grid nodes, thus avoiding the use of Lagrangian interpolation and reducing computational complexity. In this case, the first and second EFDMs, and their corresponding REMs for Eqs (1.2a)–(1.2c) can be easily adapted to the numerical solution of Eq (5.1) by replacing $ f(U^{k}_{i, j} $, $ U^{k-n_{1}}_{i, j} $, $ {\mathbf x}_{i, j}, t_{k}) $ with $ f(U^{k}_{i, j} $, $ U^{k-n_{1}}_{i, j} $, $ U^{k-n_{2}}_{i, j} $, $ \ldots $, $ U^{k-n_{L}}_{i, j} $, $ {\mathbf x}_{i, j} $, $ t_{k}) $. Also, corresponding theoretical results are derived by using similar analytical methods as long as $ f(u({\mathbf x}, t) $, $ u({\mathbf x}, t-s_{1}) $, $ u({\mathbf x}, t-s_{2}) $, $ \ldots $, $ u({\mathbf x}, t-s_{L}) $, $ {\mathbf x} $, $ t) $ satisfies $ \varepsilon_{0} $-continuous with respect to its first, second, $ \ldots, L+1 $-th variables.

As $ s_{l}(t) $ ($ l = 1, 2, \cdots, L $) are variable delays, we should use non-constrained temporal grid, i.e. $ \tau = T/n $. The linear Lagrange interpolation is used to approximate the nonlinear delayed terms. Let $ s_{l}(t_{k}) $ $ = P_{l}\tau+\theta_{l}\tau $, $ (P_{l}\in \mathbb{Z}^{+} $ $ 0\le\theta_{l}\le 1) $. Then delay term $ u({\mathbf x}, t-s_{l}(t_{k})) $ can be approximated by applying the linear Lagrange interpolation as follows $ u({\mathbf x}_{i, j}, t-s_{l}(t_{k})) $ $ = (1-\theta_{l})u^{k-P_{l}}_{i, j} $ $ +\theta_{l}u^{k-P_{l}-1}_{i, j} $ $ +O(\tau^2) $. Write $ \widehat{U}^{k, l}_{i, j} = $ $ = (1-\theta_{l})U^{k-P_{l}}_{i, j} $ $ +\theta_{l}U^{k-P_{l}-1}_{i, j} $, which is viewed as a second-order approximations to $ u({\mathbf x}_{i, j}, t-s_{l}(t_{k})) $. Consequently, the first and second EFDMs for Eqs (1.2a)–(1.2c) can be easily adapted to the numerical solution of Eq (5.1) by replacing $ f(U^{k}_{i, j}, U^{k-n_{1}}_{i, j}, {\mathbf x}_{i, j}, t_{k}) $ with $ f(U^{k}_{i, j}, \widehat{U}^{k, 1}_{i, j}, \widehat{U}^{k, 2}_{i, j}, \ldots, \widehat{U}^{k, L}_{i, j}, {\mathbf x}_{i, j}, t_{k}) $. Also, the convergence analysis of the modified first and second EFDMs can be constructed using the analytical techniques similar to those developed in Theorem 3.1 and Theorem 4.1.

6.   Numerical experiments

This section is devoted to numerical experiments. Three numerical examples are solved by the proposed methods. For convenience, denote the first EFDM in Eqs (3.3a)–(3.3c) by EFDM-I, and the second EFDM (4.4a)–(4.4c) by EFDM-II, respectively. Meanwhile, REM-I and REM-II are used to represent EFDM-I combined with REM (3.25), and EFDM-II combined with REM (4.22), respectively.

To facilitate numerical calculations, the same meshsizes in the spatial directions are adopted (i.e., $ h_{x} = h_{y} = h $). Besides, as EFDM-I and REM-I are used, we should carefully choose temporal and spatial meshsizes to satisfy $ r_{x}^{2}+r_{y}^{2} < 1 $ for EFDM-I and $ r_{x}^{2}+r_{y}^{2} < \frac{1}{4} $ for REM-I. By Theorem 3.2, Theorem 4.1, Theorem 4.2, and Remarks I–III, $ h_{x} = h_{y} = h $ and $ \tau = h^{2} $ are often applied for REM-I, EFDM-II and REM-II to obtain numerical solutions with optimally convergent rates.

The errors $ ME = \max\limits_{0\le k\le n}\|e^k\|_{\infty} $, $ LE = \max\limits_{0\leq k\leq n}||e^k||^2 $ and $ HE = \max\limits_{0\leq k\leq n}\|e^k\|_1^2 $ and corresponding orders denoted by $ r_{\infty} $, $ r_{L^2} $ and $ r_{H^1} $, respectively, and CPU time in seconds are used to test the accuracy and performance of our algorithms. All computer programmes are carried out by MATLAB R2016a.

Example 6.1. Set $ \Omega = [0, 1]\times[0, 1] $. On $ \Omega\times[-s, T] $, consider the numerical solutions of the following delayed wave equation

with Dirichlet initial-boundary conditions by using our numerical methods. Here, $ f({\mathbf x}, t) = -(\sin{x}+\cos{y})^2(\sin{t})^2(\sin^{t-s})^2 $. Initial-boundary-values conditions (IBVCs) are determined by its exact solution $ u({\mathbf x}, t) = (\sin{x}+\cos{y})\sin{t} $.

In Table 1, $ \tau = \frac{h}{2} $ is adopted to satisfy the condition $ r_{x}^{2}+r_{y}^{2} < 1 $. Table 1 shows that EFDM-I is second-order accurate in time.

In Table 2 and Table 3, timestep $ \tau = 1\times 10^{-4} $ is taken to test the spatial accuracies of EFDM-I and REM-I. Obviously, Table 2 and Table 3 exactly show that EFDM-I and REM-I are second-order and fourth-order accurate in space, respectively.

In Table 4, $ \tau = h^{2} $ is taken. By Theorem 3.1, we have $ \max_{-n_{1}\le k \le n}\{\|u^{k}-(U_{E})^{k}\|_{1}\} $ $ = O(\tau^{2}+h^{4}_{x}+h^{4}_{y}) $ $ = O(h^{4}) $ as $ \tau = h^{2} $. It is clearly observed from Table 4 that numerical solutions obtained by REM-I is convergent with an order of $ O(h^{4}) $ in $ L^{2} $-, $ H^{1} $- and $ L^{\infty} $-norms.

In Table 5 and Table 6, $ \tau = h^{2} $ is used. By Remak II and Theorem 4.1, we can infer $ \max_{-n_{1}\le k \le n}\|e^{k}\|_{1} = O(h^{2}) $ provided by the EFDM-II with $ \tau = h^{2} $. By Remak III and Theorem 4.2, we can derive $ \max_{-n_{1}\le k \le n}\{\|u^{k}-(U_{E})^{k}\|_{1}\} = O(h^{4}) $ yielded by the REM-II with $ \tau = h^{2} $. Table 5 and Table 6 clearly show that numerical solutions obtained by EFDM-II and REM-II possess the convergent orders of $ O(h^{2}) $ and $ O(h^{4}) $ in $ L^{2} $-, $ H^{1} $- and $ L^{\infty} $-norms, respectively.

From Table 7, we can find that to obtain much more accurate solutions, REMs often cost much lower time. For example, to achieve $ ME \approx 1.0\times 10^{-8} $, the CPU times of REM-I and REM-II are both much less than that of EFDM-I. Obviously, REM-I is the highest in terms of computational efficiency, while EFDM-II is the lowest.

Figure 1 displays that numerical solutions provided by EFDM-I and EFDM-II coincide with exact solutions very well. This illustrates our methods are accurate and efficient.

In a word, Tables 1–7 and Figure 1 exactly confirm that theoretical findings agree with numerical findings very well.

Example 6.2. Let $ \Omega = (0, 1)\times(0, 1) $. On $ \Omega\times [-s, T] $, we consider numerical solutions of the following wave equations with three constant delays

where $ s_1 = 0.1875 $, $ s_2 = 0.3125 $, $ s_3 = 0.5625 $, $ s = \max(s_{1}, s_{2}, s_{3}) $,

IBVCs are determined by its exact solution $ u(\mathbf{x}, t) = \exp(x+y)\sin(t) $.

In Table 8, $ \tau = \frac{h}{2} $ is used to satisfy the condition $ r_{x}^{2}+r_{y}^{2} < 1 $. Table 8 shows that the generalized EFDM-I is second-order temporally accurate.

In Table 9 and Table 10, $ \tau = 1.0\times 10^{-4} $ is fixed to test the accuracies of the generalized EFDM-I and generalized REM-I in spaces. Table 8 and Table 9 show that the generalized EFDM-I and generalized REM-I are second-order spatially accurate and fourth-order spatially accurate, respectively.

By Remark I and Theorem 3.2, we can find that the generalized REM-I possess the convergent rate of $ O(h^{4}) $ in $ H^{1} $-norm as $ h_{x} = h_{y} = h $ and $ \tau = h^{2} $. Table 11 exactly confirms that numerical solutions are convergent with an order of $ O(h^{4}) $ in $ L^{2} $-, $ L^{\infty} $- and $ H^{1} $-norms as $ \tau = h^{2} $ which satisfes $ r_{x}^{2}+r_{y}^{2} < 1/4 $.

By Remark II and Theorem 4.1, we can derive that the generalized EFDM-II possesses the convergent rate of $ O(h^{2}) $ in $ H^{1} $-norm as $ h_{x} = h_{y} = h $ and $ \tau = h^{2} $. By Remark III and Theorem 4.2, we can deduce that the generalized REM-II owns the convergent order of $ O(h^{4}) $ in $ H^{1} $-norm as $ h_{x} = h_{y} = h $ and $ \tau = h^{2} $. Table 12 and Table 13 exactly confirm that the generalized EFDM-II and the generalized REM-II have the convergent rates of $ O(h^{2}) $ and $ O(h^{4}) $ with respect to $ L^{2} $-, $ L^{\infty} $- and $ H^{1} $-norms, respectively.

Table 14 confirms that REM-I is much more efficient than EFDM-I, and REM-II is also much more efficient than EFDM-II. This exactly shows REMs can greatly improve the computational efficiency indeed. Also, Table 14 shows that REM-I is the most efficient algorithm, in contrast, EFDM-II is the lowest.

Figure 2 exhibits that both of numerical solutions given by the generalized REM-I and the generalized REM-II agree with exact solutions very well. This confirms that the current REMs are accurate and efficient.

In a word, Tables 8–14 and Figure 2 illustrate that the extensions of EFDM-I, REM-I, EFDM-II and REM-II to solve wave equations with several constant delays are practicable and preserve the accuracies of the original algorithms.

Example 6.3. Let $ \Omega = (0, 1)\times(0, 1) $. To further compare EFDM-I with EFDM-II, they are used to solve wave equations with two variable delays as follows

where

Numerical results are exhibited in Table 15 and Table 16, respectively. The exact solution to this problem is $ u(\mathbf{x}, t) = \exp(x+y+t) $, by which we can get IBVCs. Here, we take the values of the parameters as follows: $ \beta = 0.5 $, $ \theta = 3 $, $ d = 2 $, $ b = 1 $, $ \nu_1(t) = \frac{1}{2}|\sin(t)| $, and $ \nu_2(t) = \frac{1}{4}[1-\cos(t)] $. From Table 15 and Table 16, we can obtain several conclusions as follows. (1) As stable condition $ r_{x}^{2}+r_{y}^{2} < 1 $ is satisfied, numerical solutions obtained by the generalized EFDM-I with linear Lagrangian interpolation are convergent with an order of $ O(h^{2}) $ in $ L^{2} $-, $ H^{1} $- and $ L^{\infty} $-norms in the case of $ a = 1 $ and $ \tau = h^{2} $. Besides, as $ \tau = h^{2} $, numerical solutions obtained by the generalized EFDM-II with linear Lagrangian interpolation converge to exact solution with an order of $ O(h^{2}) $ in the cases of $ a = 1 $, $ a = 10 $, and $ a = 50 $. These results show that by introducing linear Lagrangian interpolation, the extensions of the EFDM-I and EFDM-II to solve wave equations with variable delays are workable and efficient with the same accuracies as the original algorithms. (2) By Theorem 3.1 and Theorem 4.1, we can find that EFDM-II does not have any restriction on $ r_{x} $ and $ r_{y} $, in contrast, $ r_{x}^{2}+r_{y}^{2} < 1 $ is necessary for EFDM-I. This finding is sufficiently confirmed by Table 15 and Table 16. For example, from Table 15 and Table 16, we can find that for larger $ a $, much finer grids are needed to obtain accepted solution for EFDM-I, while, there are no grid requirements for EFDM-II whatever $ a $ is. (3) Using the same meshsizes, which fulfill $ r_{x}^{2}+r_{y}^{2} < 1 $, numerical solutions yielded by EFDM-I more accurate than those provided by EFDM-II. Thus, EFDM-I and EFDM-II own the respective advantages of themselves. We use them according to our demands.

7.   Conclusion

In this paper, two explicit difference schemes have been derived for solving nonlinear wave equation with delays. The first EFDM (3.3a)–(3.3c) has been designed by adapting the standard explicit difference scheme for linear wave equations. It is conditionally convergent with an order of $ O(\tau^{2}+h^{2}_{x}+h^{2}_{y}) $ in $ H^{1} $-norm as $ r^{2}_{x}+r^{2}_{y} < 1 $. Besides, a spatial REM (3.25) has been derived to improve computational efficiency. The REM (3.25) is also conditionally convergent with a convergent rate of $ O(\tau^{2}+h^{4}_{x}+h^{4}_{y}) $ as $ 4(r^{2}_{x}+r^{2}_{y}) < 1 $. Although these two methods are conditionally stable, we can get numerical solutions with optimal convergent rate as long as optimal temporal and spatial grids are taken by stable conditions.

The second EFDM (4.4a)–(4.4c) has been proposed by generalizing the Du Fort-Frankel scheme which initially was proposed for parabolic equation with periodic boundary conditions. It is conditionally consistent. Namely, for its truncation error $ (R_{2})^{k}_{i, j} $ (see Eq (4.3)), as $ \tau = o(h_{x}) $, $ \tau = o(h_{y}) $, and $ h_{x} $ and $ h_{y} $ tend to zero, the truncation error $ (R_{2})^{k}_{i, j} $ tends to zero. Although numerical solutions obtained by the second EFDM (4.4a)–(4.4c) are convergent with an order of $ O(\tau^{2}+h^{2}_{x}+h^{2}_{y} $ $ +(\frac{\tau}{h_{x}})^{2}+(\frac{\tau}{h_{y}})^{2}) $ in $ H^{1} $-norm under the conditions of $ \tau = o(h^{\frac{3}{2}}) $ ($ h = \max(h_{x}, h_{y}) $) and other grids conditions, it is suggested that $ h_{x} = h_{y} = h $ and $ \tau = h^{2} $ are applied to obtain numerical solutions with an optimal convergent rate of $ O(h^{2}) $ $ H^{1} $-norm. Besides, for the second EFDM (4.4a)–(4.4c), a new REM (4.22) has been derived to improve computational accuracy. As far as we know, there has been no numerical studies of nonlinear wave equations with delay by the Du Fort-Frankel scheme combined with REMs. Thus, the derivation and analyses of this new REM is a main contribution of this study. What's more, The second EFDM (4.4a)–(4.4c) and the corresponding REM for Eq (4.22) have no any restrictions on $ r_{x} $ and $ r_{y} $.

Finally, numerical results illustrate the correctness of the theoretical results, the performance of the algorithms and the comparisons among the current algorithms. The current algorithms own the respective advantages of themselves. We use them according to our demands.

Recently, although energy-preserving Du Fort-Frankel schemes have been developed for sine-Gordon equation ^[47], coupled sine-Gordon equations ^[47] and Schrödinger with wave operator ^[48], we noted that there is few structuring-preserving numerical methods for nonlinear wave equations with delay. In future, we will consider the developments and analyses of structure-preserving Du Fort-Frankel schemes for nonlinear wave equations with delays.

Acknowledgment

This work is partially supported by the National Natural Science Foundation of China (Grant No. 11861047), Natural Science Foundation of Jiangxi Province for Distinguished Young Scientists (No. 20212ACB211006), Natural Science Foundation of Jiangxi Province (No. 20202$B$ ABL 201005).

Authors are deeply grateful to Editors, Su Meng Editorial Assistant and the anonymous referees for their insightful comments and helpful suggestions, which have greatly improved this manuscript.