Research article Special Issues

Explicit Richardson extrapolation methods and their analyses for solving two-dimensional nonlinear wave equation with delays

  • In this study, we construct two explicit finite difference methods (EFDMs) for nonlinear wave equation with delay. The first EFDM is developed by modifying the standard second-order EFDM used to solve linear second-order wave equations, of which stable requirement is accepted. The second EFDM is devised for nonlinear wave equation with delay by extending the famous Du Fort-Frankel scheme initially applied to solve linear parabolic equation. The error estimations of these two EFDMs are given by applying the discrete energy methods. Besides, Richardson extrapolation methods (REMs), which are used along with them, are established to improve the convergent rates of the numerical solutions. Finally, numerical results confirm the accuracies of the algorithms and the correctness of theoretical findings. There are few studies on numerical solutions of wave equations with delay by Du Fort-Frankel-type scheme. Therefore, a main contribution of this study is that Du Fort-Frankel scheme and a corresponding new REM are constructed to solve nonlinear wave equation with delay, efficiently.

    Citation: Dingwen Deng, Jingliang Chen. Explicit Richardson extrapolation methods and their analyses for solving two-dimensional nonlinear wave equation with delays[J]. Networks and Heterogeneous Media, 2023, 18(1): 412-443. doi: 10.3934/nhm.2023017

    Related Papers:

    [1] Farman Ali Shah, Kamran, Dania Santina, Nabil Mlaiki, Salma Aljawi . Application of a hybrid pseudospectral method to a new two-dimensional multi-term mixed sub-diffusion and wave-diffusion equation of fractional order. Networks and Heterogeneous Media, 2024, 19(1): 44-85. doi: 10.3934/nhm.2024003
    [2] Rong Huang, Zhifeng Weng . A numerical method based on barycentric interpolation collocation for nonlinear convection-diffusion optimal control problems. Networks and Heterogeneous Media, 2023, 18(2): 562-580. doi: 10.3934/nhm.2023024
    [3] Yaxin Hou, Cao Wen, Yang Liu, Hong Li . A two-grid ADI finite element approximation for a nonlinear distributed-order fractional sub-diffusion equation. Networks and Heterogeneous Media, 2023, 18(2): 855-876. doi: 10.3934/nhm.2023037
    [4] Narcisa Apreutesei, Vitaly Volpert . Reaction-diffusion waves with nonlinear boundary conditions. Networks and Heterogeneous Media, 2013, 8(1): 23-35. doi: 10.3934/nhm.2013.8.23
    [5] Min Li, Ju Ming, Tingting Qin, Boya Zhou . Convergence of an energy-preserving finite difference method for the nonlinear coupled space-fractional Klein-Gordon equations. Networks and Heterogeneous Media, 2023, 18(3): 957-981. doi: 10.3934/nhm.2023042
    [6] Yinlin Ye, Hongtao Fan, Yajing Li, Ao Huang, Weiheng He . An artificial neural network approach for a class of time-fractional diffusion and diffusion-wave equations. Networks and Heterogeneous Media, 2023, 18(3): 1083-1104. doi: 10.3934/nhm.2023047
    [7] Yves Achdou, Victor Perez . Iterative strategies for solving linearized discrete mean field games systems. Networks and Heterogeneous Media, 2012, 7(2): 197-217. doi: 10.3934/nhm.2012.7.197
    [8] Xiaoqian Gong, Alexander Keimer . On the well-posedness of the "Bando-follow the leader" car following model and a time-delayed version. Networks and Heterogeneous Media, 2023, 18(2): 775-798. doi: 10.3934/nhm.2023033
    [9] Kexin Li, Hu Chen, Shusen Xie . Error estimate of L1-ADI scheme for two-dimensional multi-term time fractional diffusion equation. Networks and Heterogeneous Media, 2023, 18(4): 1454-1470. doi: 10.3934/nhm.2023064
    [10] Li-Bin Liu, Limin Ye, Xiaobing Bao, Yong Zhang . A second order numerical method for a Volterra integro-differential equation with a weakly singular kernel. Networks and Heterogeneous Media, 2024, 19(2): 740-752. doi: 10.3934/nhm.2024033
  • In this study, we construct two explicit finite difference methods (EFDMs) for nonlinear wave equation with delay. The first EFDM is developed by modifying the standard second-order EFDM used to solve linear second-order wave equations, of which stable requirement is accepted. The second EFDM is devised for nonlinear wave equation with delay by extending the famous Du Fort-Frankel scheme initially applied to solve linear parabolic equation. The error estimations of these two EFDMs are given by applying the discrete energy methods. Besides, Richardson extrapolation methods (REMs), which are used along with them, are established to improve the convergent rates of the numerical solutions. Finally, numerical results confirm the accuracies of the algorithms and the correctness of theoretical findings. There are few studies on numerical solutions of wave equations with delay by Du Fort-Frankel-type scheme. Therefore, a main contribution of this study is that Du Fort-Frankel scheme and a corresponding new REM are constructed to solve nonlinear wave equation with delay, efficiently.



    Let Ω= [b1,b2] ×[d1,d2], where b1, b2, d1 and d2 are constant numbers and belong to R. Then we denote Ω=(b1,b2)×(d1,d2) R2 and 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]

    2ut2a2(2ux2+2uy2)=f(u(x,t),u(x,ts),x,t),(x,t)Ω×[0,T], (1.1a)
    u(x,t)=ϕ(x,t),(x,t)Ω×[s,0], (1.1b)
    u(x,t)=ψ(x,t),(x,t)Ω×[0,T] (1.1c)

    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 x for brevity. The exact solution to the problem (1.1a)–(1.1c) is denoted by u(x,t). Furthermore, we assume that f(u(x,t),u(x,ts),x,t), ϕ(x,t) and φ(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

    utuxx=f(u(x,t),u(x,ts),x,t),(x,t)Ω×[0,T], (1.2a)
    u(x,t)=ϕ(x,t),(x,t)Ω×[s,0], (1.2b)
    u(x,t)=ψ(x,t),(x,t)Ω×[0,T], (1.2c)

    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

    f(u(x,t),u(x,ts),x,t)=αu(x,t)βθdu(x,t)θd+ud(x,t)+2βθdu(x,ts)exp(γs)θd+ud(x,ts),

    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 s0. Parameters α, β, θ, γ and d(1,) are positives constants. Besides, as

    f(u(x,t),u(x,ts),x,t)=νu(x,t)+ϱu(x,ts)exp(ξu(x,ts)),

    Eq (1.2a) is known as diffusive Nicholson's blowflies equation (cf. [14,15,16]). Here, u(x,t), ϱ, 1/ξ, ν 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

    f(u(x,t),u(x,ts),x,t)=u(x,t)[1u(x,ts)],

    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 θ-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.

    To begin with, we divide the spatial domain Ω= [b1,b2] ×[d1,d2] R2, where b1, b2, d1 and d2 are constant numbers. Let hx=(b2b1)/m1, hy=(d2d1)/m2, (m1,m2Z+) 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 τ (τ=s/n1=T/n) (n1,nZ+) is used. If the restricted grid is not used, Lagrange interpolation formula will be applied to approximate the delay term u(x,ts). 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 tk=kτ, h=max{hx,hy}, xi=b1+ihx, yj=d1+jhy and (xi,yj) =xi,j, the domain Ω×[0,T] is covered by Ωhτ =ˉΩh×Ωτ, where

    ˉΩh={xi,j0im1,0jm2},Ωτ={tkn1kn}.

    Besides, we introduce

    Ωh={xi,j1im11,1jm21},Γh=ˉΩhΩh.

    On Ωhτ, further define

    Sh={U|U={Ui,j0im1,0jm2}},S0h={UU={Ui,j0im1,0jm2}andwhenxi,jΓh,Ui,j=0}.

    Then, on Ωh×Ωτ, we further introduce the following difference notations

    δ2tUki,j=Uk+1i,j2Uki,j+Uk1i,jτ2,δˆtUki,j=Uk+1i,jUk1i,j2τ,δtUk12i,j=Uki,jUk1i,jτ,δ2xUki,j=Uki+1,j2Uki,j+Uki1,jh2x,δ2yUki,j=Uki,j+12Uki,j+Uki,j1h2y,δxUki12,j=Uki,jUki1,jhx,δyUki,j12=Uki,jUki,j1hy,ΔhUi,j=(δ2x+δ2y)Ui,j.

    Finally, U, V S0h, we define the norms as follows

    U,V=hxhym11i=1m21j=1Ui,jVi,j,U,V1,x=hxhym1i=1m21j=1δxUi12,jδxVi12,j,U,V1,y=hxhym11i=1m2j=1δyUi,j12δyVi,j12,U=U,U,U=maxi,jΩh|Ui,j|,δxU=U,U1,x,δyU=U,U1,y,|U|1=δxU2+δyU2,U1=U2+|U|21.

    Define grid functions uki,j=u(xi,j,tk), xi,jΩh, n1kn. The approximation of uki,j is denoted by Uki,j.

    Applying second-order centered difference formulas to approximate temporal and spatial derivatives at point (xi,j,tk) yields that

    δ2tuki,ja2(δ2xuki,j+δ2yuki,j)=f(uki,j,ukn1i,j,xi,j,tk)+(R1)ki,j,xi,jΩ,0kn, (3.1)

    where

    (R1)ki,j=τ2124ut4(xi,j,tk)a2h2x124ux4(xi,j,tk)a2h2y124uy4(xi,j,tk)+O(τ4+h4x+h4y). (3.2)

    Omitting the small term (R1)ki,j, and then replacing uki,j with Uki,j in Eq (3.1), we get the first EFDM as follows

    δ2tUki,ja2(δ2xUki,j+δ2yUki,j)=f(Uki,j,Ukn1i,j,xi,j,tk),xi,jΩh,0kn, (3.3a)
    Uki,j=ϕ(xi,j,tk),xi,j¯Ωh,n1k0, (3.3b)
    Uki,j=ψ(xi,j,tk),xi,jΩh,0kn. (3.3c)

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

    Besides, we give an assumption for f(μ, υ, x,t) as follow.

    Assumption A. Let u(x,t) be the exact solution to the problem (1.2a)–(1.2c). Then, there are positive constants c1 and ε0, and εl (l=1,2) satisfying |εl|<ε0, (l=1,2), it holds that

    f(u(x,t)+ε1,u(x,ts)+ε2,x,t)f(u(x,t),u(x,ts),x,t)∣≤c1(ε1+ε2).

    Moreover, denote eki,j=uki,jUki,j, 0im1, 0jm2, n1kn and h=max{hx,hy}, and give several lemmas used later.

    Lemma 3.1 [45] Let VS0h, we have

    [6(b2b1)2+6(d2d1)2]V2|V|21,|V|21V21[1+16(b2b1)2+6(d2d1)2]|V|21,h2xδxV24V2,h2yδyV24V2.

    Lemma 3.2 [26] Let both A and B be positive constants, {Fk|k0} be a non-negative sequence and satisfy the following inequality Fk+1A+BτKk=0Fk. Then we can obtain max0kK+1FkAexp(B(K+1)τ). Moreover, if Fk+1A+BτK+1k=0Fk, we have max0kK+1FkAexp(2B(K+1)τ) as long as τ is sufficiently small such that Bτ12.

    Lemma 3.3 [46] Let rχ=(|a|τ)/hχ (χ=x or y), Ek=||δtek+12||2+a2ek,ek+11,x+a2ek,ek+11,y, c2=1r2xr2y. Then when r2x+r2y<1, the following inequalities

    δtek+122Ekc2,|ek+12|21Eka2, (3.4a)
    a2|ek+1|212c2Ek, (3.4b)
    ek216(b2b1)2+6(d2d1)2|ek|21(b2b1)2(d2d1)26[(b2b1)2+(d2d1)2]2c2a2Ek1 (3.4c)

    hold.

    Proof. By using the equality αβ=[(α+β)2(αβ)2]/4 and simple computations, we have

    Ek=δtek+122+a2|ek+12|21a2[||δxek+12||2ek,ek+11,x]a2[δyek+122ek,ek+11,y]=δtek+122+a2|ek+12|21a24|ek+1ek|21=c2δtek+122+a2|ek+12|21+r2xδtek+122r2x4hxhym11i=0m21j=1(δtek+12i+1,jδtek+12i,j)2+r2yδtek+122r2y4hxhym11i=1m21j=0(δtek+12i,j+1δtek+12i,j)2c2δtek+122+a2|ek+12|21,

    which is used along with r2x+r2y<1 to yield Eq (3.4a).

    According to δxek+1i+12,j=τ2δxδtek+12i+12,j+δxek+12i+12,j, by simple computation, it is easy to find that

    (δxek+1i+12,j)2=δxek+12i+12,jδxek+1i+12,j+rx2|a|(δtek+12i+1,jδtek+12i,j)δxek+1i+12,j(δxek+12i+12,j)2+(δxek+1i+12,j)22+[rx|a|(δtek+12i+1,jδtek+12i,j)]2+(δxek+1i+12,j)2434(δxek+1i+12,j)2+12(δxek+12i+12,j)2+r2x4a2(δtek+12i+1,jδtek+12i,j)2. (3.5)

    Multiplying a2hxhy to both sides of Eq (3.5) and summing i from 1 to m11 and j from 1 to m21, we obtain

    a2δxek+123a24δxek+12+a22δxek+122+r2x2δtek+122,

    which shows that

    a2δxek+122a2δxek+122+2r2xδtek+122. (3.6)

    Also, noting δyek+1i,j+12=τ2δyδtek+12i,j+12+δyek+12i,j+12 and utilizing the analytical methods similar to Eq (3.6) get that

    a2δyek+122a2δyek+122+2r2yδtek+122. (3.7)

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

    a2|ek+1|212Ek+2(r2x+r2y)c2Ek=2c2Ek,

    which is used along with Lemma 3.1 to find that

    ek216(b2b1)2+6(d2d1)2|ek|21(b2b1)2(d2d1)26[(b2b1)2+(d2d1)2]2c2a2Ek1. (3.8)

    The proof is finished.

    Let exact solution u(x,t)C4,4(Ω×[s,T]). Then there exists positive constant c3 such that

    |(R1)ki,j|c3(τ2+h2x+h2y). (3.9)

    Theorem 3.1 Let u(x,t) C4,4(Ω×[s,T]) be the exact solution of (1.2a)–(1.2c), and Uki,j be the solution of the scheme (3.3a)–(3.3c) at the time level k. Assuming that there exist positive constants μ, σ and ϵ, such that μhhx h, μh hyh and τσh12+ϵ, we can derive that

    ek1c4(τ2+h2x+h2y) (3.10)

    under the conditions that

    r2x+r2y<1,hmin((ε0μ2σ2c4)12ϵ,ε0μ4c4),τ{4+8c23(b2b1)2(d2d1)23c22a2[(b2b1)2+(d2d1)2]}1,

    and Assumption A are valid. Here,

    c4=(1+(b2b1)2(d2d1)26[(d2d1)2+(b2b1)2])2c5a2c2,
    c5=Tc23(b2b1)(d2d1)c2exp({4+8c21(b2b1)2(d2d1)23c22a2[(b2b1)2+(d2d1)2]}T).

    Proof. For convenience, denote f(uki,j)=f(uki,j,ukn1i,j,xi,yj,tk),f(Uki,j)=f(Uki,j,Ukn1i,j,xi,yj,tk) without confusion. Subtracting Eq (3.3a) from Eq (3.1), the error equations are derived as follows

    δ2teki,ja2(δ2xeki,j+δ2yeki,j)=f(uki,j)f(Uki,j)+(R1)ki,j,xi,jˉΩh,0kn, (3.11a)
    eki,j=0,xi,jˉΩh,n1k0, (3.11b)
    eki,j=0,xi,jΩh,0kn. (3.11c)

    As eki,j=0 (n1k0), it is obvious that Eq (3.10) holds for n1k0. Assuming that Eq (3.10) is valid for n1kl, we will prove that it is also true for k=l+1.

    As μhhxh and μhhyh, it is easy to find that

    1hx1μh,1hy1μh. (3.12)

    Thus, using τσh12+ϵ and Eq (3.12), we obtain

    τ2hxhy+h32xhy+h32yhxσ2h1+2ϵμhμh+h32μh+h32μhσ2μh2ϵ+2hμ. (3.13)

    Besides, applying hmin{((ε0μ2σ2c4) 12ϵ, ε0μ4c4)}, it is easy to find that

    σ2μh2ϵε02c4,σ2μh2ϵε02c4. (3.14)

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

    |eki,j|h12xh12yekh12xh12yek1h12xh12yc4(τ2+h2x+h2y)=c4(τ2hxhy+h32xhy+h32yhx)c4(τ2μh+h32μh+h32μh)c4(σ2h2ϵμ+2hμ)<ε0,n1kl,

    which follows from Assumption A that

    |f(uki,j)f(Uki,j)|c1(|eki,j|+|ekn1i,j|),n1kl. (3.15)

    Taking inner product with Eq (3.11a) by 2δˆtek gives that

    δ2tek,2δˆteka2δ2xek,2δˆteka2δ2yek,2δˆtek=f(uk)f(Uk),2δˆtek+(R1)k,2δˆtek. (3.16)

    In what follows, we estimate each terms of Eq (3.16) step by step. Using the equality α2β2 =(αβ)(α+β) gives that

    δ2tek,2δˆtek=1τ(δtek+122δtek122). (3.17)

    Applying the discrete Green formula yields that

    a2δ2xek,2δˆtek=a2τhxhym11i=0m2j=0(δxeki+12,j)(δxek+1i+12,jδxek1i+12,j)=a2τ[ek,ek+11,xek,ek11,x]. (3.18)

    Similar to Eq (3.18), we also arrive at

    a2δ2yek,2δˆtek=a2τ[ek,ek+11,yek,ek11,y]. (3.19)

    Next, we estimate each terms at the right of Eq (3.16). Applying Eq (3.15), the inequalities 2αβ1εα2+εβ2 and (α+β)22(α2+β2), and Lemma 3.3, we have

    f(Uk)f(uk),2δˆtek2c21c2(ek2+ekn12)+c22(δtek+122+δtek122)2c21c2(ek||2+ekn12)+12(Ek+Ek1). (3.20)

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

    (R1)k,2δˆtekc23(b2b1)(d2d1)c2(τ2+h2x+h2y)2+12(Ek+Ek1). (3.21)

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

    EkEk1τ(Ek+Ek1)+2c21τc2(||ek||2+||ekn1||2)+c22(b2b1)(d2d1)τc2(τ2+h2x+h2y)2. (3.22)

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

    Ekτkγ=1Eγ+τk1γ=0Eγ+τ(2c1c2a)2k1γ=0(b2b1)2(d2d1)26[(b2b1)2+(d2d1)2]Eγ+τ(2c1c2a)2kn11γ=n1(b2b1)2(d2d1)26[(b2b1)2+(d2d1)2]Eγ+τkγ=0c23(b2b1)(d2d1)c2(τ2+h2x+h2y)2τ{2+8c21c22a2(b2b1)2(d2d1)26[(b2b1)2+(d2d1)2]}kγ=0Eγ+c23(b2b1)(d2d1)kτc2(τ2+h2x+h2y)2τ{2+4c21(b2b1)2(d2d1)23c22a2[(b2b1)2+(d2d1)2]}kγ=0Eγ+c23(b2b1)(d2d1)Tc2(τ2+h2x+h2y)2. (3.23)

    Hence, applying Lemma 3.2 to Eq (3.23), it follows from τ{4+8c23(b2b1)2(d2d1)23c22a2[(b2b1)2+(d2d1)2]}1 that

    Ekc5(τ2+h2x+h2y)2,0kl. (3.24)

    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.

    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 τ, hx and hy.

    Theorem 3.2 Suppose that u(x,t) C6,4(Ω×[s,T]) is the exact solution of Eqs (1.2a)–(1.2c), and Uki,j(τ,hx,hy) is the numerical solution of the 1st EFDM (3.3a)–(3.3c) at (xi,j,tk) using meshsizes hx, hy and τ. Define a REM as follows

    (UE)ki,j={43Uk2i,2j(τ,hx2,hy2)13Uki,j(τ,hx,hy),xi,jˉΩh,1kN,ψ(xi,j,tk),xi,jˉΩh,mk0. (3.25)

    Then, as 4a2[(τhx)2+(τhy)2]<1, and grid parameters τ, hx and hy are small enough, we obtain

    uk(UE)k1C(τ2+h4x+h4y). (3.26)

    Proof. Let r1(x,t)=a2124ux4(x,t) and r2(x,t)=a2124uy4(x,t). Then a combination of Taylor expansion formula with (3.2) gives that

    (R1)ki,j=τ2124ut4(xi,j,tk)h2xr1(xi,j,tk)h2yr2(xi,j,tk)+O(τ4+h4x+h4y). (3.27)

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

    δ2teki,ja2Δheki,j=f(uki,j,ukn1i,j,xi,j,tk)f(Uki,j,Ukn1i,j,xi,j,tk)+τ2124ut4(xi,j,tk)h2x(r1)ki,jh2y(r2)ki,j+O(τ4+h4x+h4y),xi,jˉΩh,0kn; (3.28a)
    eki,j=0,xi,jˉΩh,n1k0; (3.28b)
    eki,j=0,xi,jΩh,0kn. (3.28c)

    Furthermore, two auxiliary problems are introduced to derive the asymptotic expansion of the numerical solutions. Namely, let v1(x,t) and v2(x,t) be the exact solutions of the following initial-boundary-value problems (IBVPs)

    (v1)tta2Δv1=r1(x,t)+fμ(u(x,t),u(x,ts),x,t)v1(x,t)+fυ(u(x,t),u(x,ts),x,t)v1(x,ts),(x,t)Ω×[0,T]; (3.29a)
    v1(x,t)=0,(x,t)Ω×[s,0]; (3.29b)
    v1(x,t)=0,(v1)t(x,t)=0,xΩ,t×[0,T]; (3.29c)
    (v2)tta2Δv2=r2(x,t)+fμ(u(x,t),u(x,ts),x,t)v2(x,t)+fυ(u(x,t),u(x,ts),x,t)v2(x,ts)(x,t)Ω×[0,T]; (3.30a)
    v2(x,t)=0,(x,t)Ω×[s,0]; (3.30b)
    v2(x,t)=0,(v2)t(x,t)=0,xΩ,t×[0,T], (3.30c)

    respectively.

    Denote (v1)ki,j =v1(xi,j,tk) and (v2)ki,j =v2(xi,j,tk). 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

    δ2t(v1)ki,ja2Δh(v1)ki,j=r1(xi,j,tk)+fμ(uki,j,ukn1i,j,xi,j,tk)(v1)ki,j+fυ(uki,j,ukn1i,j,xi,j,tk)(v1)kn1i,j+O(τ2+h2x+h2y),xi,jˉΩh,0kn; (3.31a)
    (v1)ki,j=0,xi,jˉΩh,n1k0; (3.31b)
    (v1)ki,j=0,xi,jΩh,0kn, (3.31c)
    δ2t(v2)ki,ja2Δh(v2)ki,j=r2(xi,j,tk)+fμ(Uki,j,Ukn1i,j,xi,j,tk)(v2)ki,j+fυ(Uki,j,Ukn1i,j,xi,j,tk)(v2)kn1i,j+O(τ2+h2x+h2y),xi,jˉΩh,0kn; (3.32a)
    (v2)ki,j=0,xi,jˉΩh,n1k0; (3.32b)
    (v2)ki,j=0,xi,jΩh,0kn, (3.32c)

    respectively.

    Denote pki,j=eki,j+h2x(v1)ki,j+h2y(v2)ki,j, 0im1, 0jm2, n1kn. Then, multiplying Eqs (3.31b)–(3.31d) by h2x and Eqs (3.32b)–(3.32d) by h2y, respectively, and then adding the obtained results to Eqs (3.28b)–(3.28d), we have

    δ2tpki,ja2Δhpki,j=(˜H1)ki,j+(~R1)ki,j,xi,jˉΩh,0kn; (3.33a)
    pki,j=0,xi,jˉΩh,n1k0; (3.33b)
    pki,j=0,xi,jΩh,0kn, (3.33c)

    where

    (˜H1)ki,j=f(uki,j+h2x(v1)ki,j+h2y(v2)ki,j,ukn1i,j+h2x(v1)kn1i,j+h2y(v2)kn1i,j,xi,j,tk)f(Uki,j,ukn1i,j,xi,j,tk),

    and

    (~R1)ki,j=τ2124ut4(xi,yj,tk)+O(τ4+h4x+h4y). (3.34)

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

    f(uki,j,ukn1i,j,xi,j,tk)+h2x[fμ(uki,j,ukn1i,j,xi,j,tk)(v1)ki,j+fυ(uki,j,Ukn1i,j,xi,j,tk)(v1)kn1i,j]+h2y[fμ(uki,j,ukn1i,j,xi,j,tk)(v2)ki,j+fυ(uki,j,ukn1i,j,xi,j,tk)(v2)kn1i,j]=f(uki,j+h2x(v1)ki,j+h2y(v2)ki,j,ukn1i,j+h2x(v1)kn1i,j+h2y(v2)kn1i,j,xi,j,tk)+O(h4x+h4y+h2xh2y)

    is used.

    Then, as a2[(τhx)2+(τhy)2]<1, using the same techniques as those proposed in the proof of Theorem 3.1, we find that

    pk1C(τ2+h4x+h4y),0kn, (3.35)

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

    uk[Uk(τ,hx,hy)h2x(v1)kh2y(v2)k]1C(τ2+h4x+h4y),0kn. (3.36)

    On the other hand, from Richardson extrapolation solutions

    (UE)ki,j=43Uk2i,2j(τ,hx2,hy2)13Uki,j(τ,hx,hy),xi,jΩh,0kn, (3.37)

    it follows that

    uki,j(UE)ki,j=43{uki,j[Uk2i,2j(τ,hx2,hy2)(hx2)2(v1)ki,j(h2y2)2(v2)ki,j]}13{uki,j[Uki,j(τ,hx,hy)h2x(v1)ki,jh2y(v2)ki,j]}. (3.38)

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

    uk(UE)k143uk[Uk(τ,hx2,hy2)(hx2)2(v1)k(h2y2)2(v2)k]1+13uk[Uk(τ,hx,hy)h2x(v1)kh2y(v2)k]14C3[τ2+(hx2)4+(hy2)4]+C3(τ2+h4x+h4y)5C3(τ2+h4x+h4y), (3.39)

    holds as long as 4a2[(τhx)2+(τhy)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(τ2+h2x+h2y) as r2x+r2y<1. Besides, Theorem 4.1 shows that REM (3.25) is also conditionally convergent with an order of O(τ2+h24+h4y) as 4(r2x+r2y)<1 (see Eq (3.39)). Besides, as hx=hy=h and τ=h2, Richardson extrapolation solutions have a convergent rate of O(h4). In Section 6, to obtain numerical solutions with optimal convergent rates, we take temporal and spatial grids according to these conditions.

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

    Using the following difference formulas

    uxx(xi,j,tk)=δ2xuki,jτ2h2xδ2tuki,j+h2x124ux4(xi,j,tk)+τ2h2x2ut2(xi,j,tk)+O(τ4h2x+h4x), (4.1a)
    uyy(xi,j,tk)=δ2yuki,jτ2h2yδ2tuki,j+h2y124uy4(xi,j,tk)+τ2h2y2ut2(xi,j,tk)+O(τ4h2y+h4y) (4.1b)

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

    (1+r2x+r2y)δ2tuki,ja2Δhuki,j=f(uki,j,ukn1i,j,xi,j,tk)+(R2)ki,j,xi,jΩh,0kn, (4.2)

    where

    (R2)ki,j=τ2124ut4(xi,j,tk)a2h2x124ux4(xi,j,tk)a2h2y124uy4(xi,j,tk)+a2(τ2h2x+τ2h2y)2ut2(xi,j,tk)+O(τ4+h4x+h4y+τ4h2x+τ4h2y). (4.3)

    Omitting the small term (R2)ki,j and replacing uki,j by its approximation Uki,j in Eq (4.2), we get the second EFDM as follows

    (1+r2x+r2y)δ2tUki,ja2ΔhUki,j=f(Uki,j,Ukn1i,j,xi,j,tk),xi,jΩh,0kn, (4.4a)
    Uki,j=ϕ(xi,j,tk),xi,j¯Ωh,n1k0, (4.4b)
    Uki,j=ψ(xi,j,tk),xi,jΩh,0kn. (4.4c)

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

    Lemma 4.1 Let Gk=(1+r2x+r2y)δtek+122+a2ek,ek+11,x+a2ek,ek+11,y. Then the following inequalities

    δtek+122Gk,|ek+12|21Gka2, (4.5a)
    |ek+1|212(1+r2x+r2y)a2Gk, (4.5b)
    ek+121{1+(b2b1)2(d2d1)26[(b2b1)2+(d2d1)2]}2(1+r2x+r2y)a2Gk (4.5c)

    are valid.

    Proof. Using Lemma 3.1, we obtain

    r2xδtek+122a2τ24δxδtek+122=r2xδtek+122a2τ24h2xh2xδxδtek+120, (4.6a)
    r2yδtek+122a2τ24δyδtek+122=r2yδtek+122a2τ24h2yh2yδxδtek+120. (4.6b)

    Applying Eqs (4.6a), (4.6b) and the equality αβ= (α+β)2/4 (αβ)2/4, we have that

    Gk=(1+r2x+r2y)δtek+122+(a2δxek+122a2τ24h2xδxδtek+122)+(a2δyek+122a2τ24h2yδyδtek+122)δtek+122+a2|ek+12|21,

    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

    |ek+1|212|ek+12|21+2(r2x+r2y)a2||δtek+12||22(1+r2x+r2y)a2Gk, (4.7a)
    ek+12(b2b1)2(d2d1)26[(b2b1)2+(d2d1)2]|ek+1|21(b2b1)2(d2d1)2(1+r2x+r2y)3a2[(b2b1)2+(d2d1)2]Gk. (4.7b)

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

    ek+121{1+(b2b1)2(d2d1)26[(b2b1)2+(d2d1)2]}2(1+r2x+r2y)a2Gk.

    The proof is completed.

    Assuming u(x,t) C4,4(Ω×[s,T]), there exists positive constant c6 such that

    |(R2)ki,j|c6(τ2+h2x+h2y+τ2h2x+τ2h2y). (4.8)

    Theorem 4.1 Let u(x,t) C4,4(Ω×[s,T]) be exact solution of Eqs (1.2a)–(1.2c), Uki,j be the solution of the scheme Eqs (4.4a)–(4.4c) at node (xi,j,tk). Assume that there are constants μ, σ and ϵ such that μhhxh, μhhyh, τσh32+ϵ. It holds that

    ek1c7(τ2+h2x+h2y+τ2h2x+τ2h2y) (4.9)

    under the conditions that

    hmin{(με03c7σ2)12(1+ϵ),με06c7,(μ3ε06c7σ2)12ϵ},τ{4+4c21(b2b1)2(d2d1)2(1+r2x+r2y)3[(b2b1)2+(d2d1)2]a2}1

    and Assumption A hold. Here,

    c7=(1+(b2b1)2(d2d1)26[(d2d1)2+(b2b1)2])2(1+r2x+r2y)c8a2,
    c8=Tc26(b2b1)(d2d1)exp((4+4c21(b2b1)2(d2d1)2(1+r2x+r2y)3[(b2b1)2+(d2d1)2]a2)T).

    Proof. Also, let f(uki,j)=f(uki,j,ukn1i,j,xi,j,tk), f(Uki,j)=f(Uki,j,Ukn1i,j,xi,j,tk) for convenience. Subtracting Eq (4.4a) from Eq (4.2), the error equations are derived as follows

    (1+r2x+r2y)δ2teki,ja2Δheki,j=f(uki,j)f(Uki,j)+(R2)ki,j,xi,jΩh,0kn, (4.10a)
    eki,j=0,xi,jΩh,n1k0, (4.10b)
    eki,j=0,xi,jΩh,0kn. (4.10c)

    As eki,j=0 (n1k0), it is obvious that Eq (4.9) holds for n1k0. Assuming that Eq (4.9) is valid for n1kl, we will prove that it is also true for k=l+1.

    As μhhxh and μhhyh, it is easy to find that

    1hx1μh,1hy1μh. (4.11)

    Thus, using τσh32+ϵ and Eq (4.11), we obtain

    (τ2hxhy+h32xhy+h32yhx+τ2h5xhy+τ2hxh5y)σ2h3+2ϵμhμh+h32μh+h32μh+σ2h3+2ϵ(μh)5μh+σ2h3+2ϵμh(μh)5σ2μh2+2ϵ+2hμ+2σ2h2ϵμ3. (4.12)

    Besides, applying hmin((με03c7σ2)12(1+ϵ),με06c7,(μ3ε06c7σ2)12ϵ), it is easy to find that

    σ2μh2(1+ϵ)ε03c7,2hμε03c7,2σ2μ3h2ϵε03c7. (4.13)

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

    |eki,j| h12xh12yek h12xh12yek1h12xh12yc7(τ2+h2x+h2y+τ2h2x+τ2h2y)=c7(τ2hxhy+h32xhy+h32yhx+τ2h5xhy+τ2hxh5y)c7(σ2h2+2ϵμ+2hμ+2σ2h2ϵμ3)<ε0,n1kl. (4.14)

    Hence, applying Assumption A directly infers that

    |f(Uki,j)f(uki,j)|c1(|eki,j|+|ekn1i,j|),n1kl. (4.15)

    Acting the inner product with Eq (4.10a) by 2δˆtek yields that

    (1+r2x+r2y)δ2tek,2δˆteka2δ2xek,2δˆteka2δ2yek,2δˆtek=f(uk)f(Uk),2δˆtek+(R2)k,2δˆtek. (4.16)

    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

    (1+r2x+r2y)(δ2tek,2δˆtek)a2(δ2xek,2δˆtek)a2(δ2yek,2δˆtek)=1+r2x+r2yτ(δtek+122δtek122)+a2τ[ek+1,ek1,xek,ek11,x]+a2τ[ek+1,ek1,yek,ek11,y]=1τ(GkGk1). (4.17)

    In what follows, we estimate each terms at the right of Eq (4.16). Applying the inequalities 2αβ1εα2+εβ2 and (α+β)22(α2+β2), Eq (4.15) and Lemma 4.1, we can get that

    f(uk)f(Uk),2δˆtek2c21(ek2+ekn12)+12(||δtek+122+||δtek122)2c21(ek2+||ekn12)+12(Gk+Gk1). (4.18)

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

    (R2)k,2δˆtekc26(b2b1)(d2d1)(τ2+h2x+h2y+τ2h2x+τh2y)2+12(Gk+Gk1). (4.19)

    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 γ, summing γ from 0 to k, we can arrive at

    Gkτ(kγ=0Gγ+k1γ=1Gγ)+τc21(b2b1)2(d2d1)2(1+r2x+r2y)3[(b2b1)2+(d2d1)2]a2k1γ=1Gγ+τc21(b2b1)2(d2d1)2(1+r2x+r2y)3[(b2b1)2+(d2d1)2]a2kn11γ=n11Gγ+τkγ=0(R2)γ2τ{2+2c21(b2b1)2(d2d1)2(1+r2x+r2y)3[(b2b1)2+(d2d1)2]a2}kγ=0Gγ+Tc26(b2b1)(d2d1)(τ2+h2x+h2y+τ2h2x+τ2h2y)2. (4.20)

    As τ{4+4c21(b2b1)2(d2d1)2(1+r2x+r2y)3[(b2b1)2+(d2d1)2]a2}1, applying Lemma 3.2 to Eq (4.20) derives that

    Gkc8(τ2+h2x+h2y+τ2h2x+τ2h2y)2,0kl. (4.21)

    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(τ2+h2x+h2y+τ2h2x+τ2h2y) in H1-norm without any restriction no rx and ry. From truncation error Eq (4.3), it follows that the second EFDM (4.4a)–(4.4c) is conditionally consistent, and (R2)ki,j tends to zero as τ=o(hx), τ=o(hy), and hx and hy tend to zero. What's more, although the condition τσh32+ϵ (σ and ϵ are both positive constant.) is necessary in Theorem 4.1, we would like to adopt the grid τ=h2 and hx=hy=h to get numerical solutions with an optimal convergent order of O(h2) in practical computations.

    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(x,t) C4,4(Ω×[s,T]) is the exact solution of Eqs (1.2a)–(1.2c), and Uki,j(τ,hx,hy) is the numerical solution of the second EFDM (4.4a)–(4.4c) at (xi,j,tk) using meshsizes hx, hy and τ. Define the following REM

    (UE)ki,j={19Uki,j(τ,hx,hy)49U2ki,j(τ2,hx,hy)49U2k2i,2j(τ2,hx2,hy2)+169U4k2i,2j(τ4,hx2,hy2),xi,jˉΩh,1nN,ψ(xi,j,tk),xi,jˉΩh.mk0, (4.22)

    Then under the conditions of Theorem 4.1, hx=O(hy), τhx and τhy 0 as grid parameters τ, hx, hy 0, we arrive at

    uk(UE)k1C[τ2+τ4+h4x+h4y+τ4h2x+τ4h2y+(τhx)4+(τhy)4], (4.23)

    as long as parameter grids τ hx and hy are small enough.

    Proof. Let w1(x,t)=1124ut4(x,t), w2(x,t)=a2124ux4(x,t), w3(x,t)=a2124uy4(x,t) and w4(x,t)=a22ut2(x,t). Then, we define the grid functions (wl)ki,j =wl(xi,j,tk), (l=1,2,3,4), xi,jˉΩh, n1kn. Thus, it follows from Eq (4.3) that

    (R2)ki,j=τ2(w1)ki,jh2x(w2)ki,jh2y(w3)ki,j+(τ2h2x+τ2h2y)(w4)ki,j+O(τ4+h4x+h4y+τ4h2x+τ4h2y). (4.24)

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

    (1+r2x+r2y)δ2teki,ja2(δ2xeki,j+δ2yeki,j)=f(uki,j)f(Uki,j)+τ2(w1)ki,jh2x(w2)ki,jh2y(w3)ki,j+(τ2h2x+τ2h2y)(w4)ki,j+O(τ4+h4x+h4y+τ4h2x+τ4h2y),xi,jΩh,0kn, (4.25a)
    eki,j=0,xi,jΩh,n1k0, (4.25b)
    ek0,j=0,ekm1,j=0,0jm2,0kn, (4.25c)
    eki,0=0,eki,m2=0,0im1,0kn. (4.25d)

    Suppose that functions v1(x,t), v2(x,t), v3(x,t) and v4(x,t) are exact solutions of the following IBVPs

    2v1t2a2(2v1x2+2v1y2)=w1(x,t)+h1(x,t),(x,t)Ω×[0,T], (4.26a)
    v1(x,t)=0,(x,t)Ω×[s,0], (4.26b)
    v1(x,t)=0,(x,t)Ω×[0,T]; (4.26c)
    2v2t2a2(2v2x2+2v2y2)=w2(x,t)+h2(x,t),(x,t)Ω×[0,T], (4.27a)
     v2(x,t)=0,(x,t)Ω×[s,0], (4.27b)
    v2(x,t)=0,(x,t)Ω×[0,T]; (4.27c)
    2v3t2a2(2v3x2+2v3y2)=w3(x,t)+h3(x,t),(x,t)Ω×[0,T], (4.28a)
    v3(x,t)=0,(x,t)Ω×[s,0], (4.28b)
    v3(x,t)=0,(x,t)Ω×[0,T]; (4.28c)
    2v4t2a2(2v4x2+2v4y2)=w4(x,t)+h4(x,t),(x,t)Ω×[0,T], (4.29a)
    v4(x,t)=0,(x,t)Ω×[s,0], (4.29b)
    v4(x,t)=0,(x,t)Ω×[0,T]; (4.29c)

    respectively, where

    h1(x,t)=fμ(u(x,t),u(x,ts),x,t)v1(x,t)+fυ(u(x,t),u(x,ts),x,t)v1(x,ts),
    h2(x,t)=fμ(u(x,t),u(x,ts),x,t)υ2(x,t)+fυ(u(x,t),u(x,ts),x,y,t)v2(x,ts), (4.30b)
    h3(x,t)=fμ(u(x,t),u(x,ts),x,t)v3(x,t)+fυ(u(x,t),u(x,ts),x,t)v3(x,ts), (4.30c)
    h4(x,t)=fμ(u(x,t),u(x,ts),x,t)v4(x,t)+fυ(u(x,t),u(x,ts),x,t)v4(x,ts). (4.30d)

    Denote (vl)ki,j=vl(xi,j,tk), 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

    (1+r2x+r2y)δ2t(v1)ki,ja2[δ2x(v1)ki,j+δ2y(v1)ki,j]=(w1)ki,j+h1(xi,j,tk)+O(τ2+h2x+h2y+τ2h2x+τ2h2y),xi,jΩh,0kn, (4.31a)
    (v1)ki,j=0,xi,jΩh,n1k0, (4.31b)
    (v1)ki,j=0,xi,jΩh,0kn, (4.31c)
    (1+r2x+r2y)δ2t(v2)ki,ja2[δ2x(v2)ki,j+δ2y(v2)ki,j]=(w2)ki,j+h2(xi,j,tk)+O(τ2+h2x+h2y+τ2h2x+τ2h2y),xi,jΩh,0kn, (4.32a)
    (v2)ki,j=0,xi,jΩh,n1k0, (4.32b)
    (v2)ki,j=0,xi,jΩh,0kn, (4.32c)
    (1+r2x+r2y)δ2t(v3)ki,ja2[δ2x(v3)ki,j+δ2y(v3)ki,j]=(w3)ki,j+h3(xi,j,tk)+O(τ2+h2x+h2y+τ2h2x+τ2h2y),xi,jΩh,0kn, (4.33a)
    (v3)ki,j=0,xi,jΩh,n1k0, (4.33b)
    (v3)ki,j=0,xi,jΩh,0kn, (4.33c)
    (1+r2x+r2y)δ2t(v4)ki,ja2[δ2x(v4)ki,j+δ2y(v4)ki,j]=(w4)ki,j+h4(xi,j,tk)+O(τ2+h2x+h2y+τ2h2x+τ2h2y),xi,jΩh,0kn, (4.34a)
    (v4)ki,j=0,xi,jΩh,n1k0, (4.34b)
    (v4)ki,j=0,xi,jΩh,0kn. (4.34c)

    Denote qki,j=eki,j+h2x(v1)ki,j+h2y(v2)ki,j+(τ2h2x+τ2h2y)(v3)ki,j+τ2(v4)ki,j, 0im1, 0jm2, n1kn. Then, multiplying h2x, h2y, (τ2h2x+τ2h2y) and τ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

    δ2tqki,ja2(δ2xqki,j+δ2yqki,j)=˜hki,j+(~R2)ki,j,xi,jˉΩh,0kn, (4.35a)
    qki,j=0,i,jˉΩh,n1k0, (4.35b)
    qki,j=0,xi,jΩh,0kn, (4.35c)

    where

    (˜h1)ki,j=f(uki,j+h2x(v1)ki,j+h2y(v2)ki,j+(τ2h2x+τ2h2y)(v3)ki,j+τ2(v4)ki,j, (4.36a)
    ukn1i,j+h2x(v1)kn1i,j+h2y(v2)kn1i,j+(τ2h2x+τ2h2y)(v3)kn1i,j+τ2(v4)kn1i,j,xi,j,tk)f(Uki,j,Ukn1i,j,xi,j,tk), (4.36b)
    (~R2)ki,j=O(τ4+h4x+h4y+τ2h2x+τ2h2y+τ4h2x+τ4h2y+((hxhy)2+(hxhy)2)τ2 (4.36c)
    +(τhx)4+(τhx)4+(τ2hxhy)2). (4.36d)

    It is worth mentioning that the following Taylor expansion formula

    f(uki,j,ukn1i,j,xi,j,tk)+h2x[fμ(uki,j,ukn1i,j,xi,j,tk)(v1)ki,j+fυ(uki,j,ukn1i,j,xi,j,tk)(v1)kn1i,j]+h2y[fμ(uki,j,ukn1i,j,xi,j,tk)(v2)ki,j+fυ(uki,j,ukn1i,j,xi,j,tk)(v2)kn1i,j]+[(τhx)2+(τhy)2][fμ(uki,j,ukn1i,j,xi,j,tk)(v3)ki,j+fυ(uki,j,ukn1i,j,xi,j,tk)(v3)kn1i,j]+τ2[fμ(uki,j,ukn1i,j,xi,j,tk)(v4)ki,j+fυ(uki,j,ukn1i,j,xi,j,tk)(v4)kn1i,j]+O(τ4+τ4h2x+τ4h2x+τ2h2x+τ2h2y+(h2xh2y+h2yh2x)τ2+h4x+h2xh2y+h4y+(τhx)4+(τ2hxhy)2+(τhy)4) (4.37)
    =f(uki,j+h2x(v1)ki,j+h2y(v2)ki,j+(τ2h2x+τ2h2y)(v3)ki,j+τ2(v4)ki,j,ukn1i,j+h2x(v1)kn1i,j+h2y(v2)kn1i,j+(τ2h2x+τ2h2y)(v3)kn1i,j+τ2(v4)kn1i,j,xi,j,tk)+O(τ4+τ4h2x+τ4h2y+τ2h2x+τ2h2y+(h2xh2y+h2yh2x)τ2+h4x+h2xh2y+h4y+(τhx)4+(τ2hxhy)2+(τhy)4)

    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 αβ(α2+β2)/2 and max((hyhx)2,(hxhy)2)μ2, we have

    qk1C(τ2+h4x+h4y+τ4+τ4h2x+τ4h2y+(τhx)4+(τhy)4),0kn, (4.38)

    where C is a positive constant and independent of grid parameters τ, hx and hy. Finally, similar to Eq (3.39), using the triangle inequality to the following equality

    uki,j(UE)ki,j=19{(uki,jUki,j(hx,hy,τ)+h2x(v1)ki,j+h2y(v2)ki,j+(τ2h2x+τ2h2y)(v3)ki,j+τ2(v4)ki,j}49{(uki,jU2ki,j(hx,hy,τ2)+h2x(v1)ki,j+h2y(v2)ki,j+((τ2)2h2x+(τ2)2h2y)(v3)ki,j+(τ2)2(v4)ki,j}49{(uki,jU2k2i,2j(hx2,hy2,τ2)+(hx2)2(v1)ki,j+(hy2)2(v2)ki,j+[(τ2)2(hx2)2+(τ2)2(hy2)2](v3)ki,j+(τ2)2(v4)ki,j}+169{(uki,jU4k2i,2j(hx2,hy2,τ4)+(hx2)2(v1)ki,j+(hy2)2(v2)ki,j+[(τ4)2(hx2)2+(τ4)2(hy2)2](v3)ki,j+(τ4)2(v4)ki,j}

    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(τ2+τ4+h4x+h4y+τ4h2x+τ4h2y+(τhx)4+(τhy)4) in H1-norm. According to the convergence rate O(τ2+τ4+h4x+h4y+τ4h2x+τ4h2y+(τhx)4+(τhy)4), the optimal estimation uk(UE)k1=O(h4), 0kn can be obtained as an optimally temporal and spatial relationship hx=hy=h and τ=h2 is used. Thus this relationship hx=hy=h and τ=h2 is applied in practical computations.

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

    ut=2ux2+υ(g(u(x,tτ1)))ux+c[f(u(x,tτ2))u(x,t)],

    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 τ1 and τ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

    {utta2(uxx+uyy)=f(u(x,t),u(x,ts1(t)),u(x,ts2(t)),,u(x,tsL(t)),x,t),(x,t)Ω×[0,T],u(x,t)=ψ(x,t),xˉΩ,t[s,0]u(x,t)=ϕ(x,t),(x,t)Ω×[0,T], (5.1)

    where sl(t)>0 (l=1,2,,L) and s=max1lLsl(t).

    As sκ(t) =sκ (κ=1,2,,L) are all constant delays, i.e. independent of temporal variable t, constrained temporal grid τ=s1/n1=s2/n2==sL/nL,(nlZ+) is adopted to make t=Ksl (KZ+) 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(Uki,j, Ukn1i,j, xi,j,tk) with f(Uki,j, Ukn1i,j, Ukn2i,j, , UknLi,j, xi,j, tk). Also, corresponding theoretical results are derived by using similar analytical methods as long as f(u(x,t), u(x,ts1), u(x,ts2), , u(x,tsL), x, t) satisfies ε0-continuous with respect to its first, second, ,L+1-th variables.

    As sl(t) (l=1,2,,L) are variable delays, we should use non-constrained temporal grid, i.e. τ=T/n. The linear Lagrange interpolation is used to approximate the nonlinear delayed terms. Let sl(tk) =Plτ+θlτ, (PlZ+ 0θl1). Then delay term u(x,tsl(tk)) can be approximated by applying the linear Lagrange interpolation as follows u(xi,j,tsl(tk)) =(1θl)ukPli,j +θlukPl1i,j +O(τ2). Write ˆUk,li,j= =(1θl)UkPli,j +θlUkPl1i,j, which is viewed as a second-order approximations to u(xi,j,tsl(tk)). 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(Uki,j,Ukn1i,j,xi,j,tk) with f(Uki,j,ˆUk,1i,j,ˆUk,2i,j,,ˆUk,Li,j,xi,j,tk). 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.

    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., hx=hy=h). Besides, as EFDM-I and REM-I are used, we should carefully choose temporal and spatial meshsizes to satisfy r2x+r2y<1 for EFDM-I and r2x+r2y<14 for REM-I. By Theorem 3.2, Theorem 4.1, Theorem 4.2, and Remarks I–III, hx=hy=h and τ=h2 are often applied for REM-I, EFDM-II and REM-II to obtain numerical solutions with optimally convergent rates.

    The errors ME=max0knek, LE=max0kn||ek||2 and HE=max0knek21 and corresponding orders denoted by r, rL2 and rH1, 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 Ω=[0,1]×[0,1]. On Ω×[s,T], consider the numerical solutions of the following delayed wave equation

    2ut2=2ux2+2uy2+u2(x,t)u2(x,ts)+f(x,t),xΩ,0tT

    with Dirichlet initial-boundary conditions by using our numerical methods. Here, f(x,t)=(sinx+cosy)2(sint)2(sints)2. Initial-boundary-values conditions (IBVCs) are determined by its exact solution u(x,t)=(sinx+cosy)sint.

    In Table 1, τ=h2 is adopted to satisfy the condition r2x+r2y<1. Table 1 shows that EFDM-I is second-order accurate in time.

    Table 1.  Numerical results obtained by EFDM-I for Example 6.1, s=0.01, T=1, (τ=0.5h).
    h ME r LE rL2 HE rH1 CPU
    1/50 3.3184×106 1.7790×106 8.2345×106 0.082
    1/100 8.3000×107 1.999 4.4483×107 2.000 2.0594×106 2.000 0.514
    1/200 2.0752×107 2.000 1.1121×107 2.000 5.1490×107 2.000 3.931
    1/400 5.1880×108 2.000 2.7803×108 2.000 1.2873×107 2.000 31.872

     | Show Table
    DownLoad: CSV

    In Table 2 and Table 3, timestep τ=1×104 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.

    Table 2.  Numerical results obtained by EFDM-I for Example 6.1, s=0.01, T=1, (τ=1×104).
    h ME r LE rL2 HE rH1 CPU
    1/2 2.2247×103 1.1124×103 4.586×103 0.045
    1/4 6.5754×104 1.759 3.5141×104 1.663 1.5731×103 1.544 0.081
    1/8 1.6883×104 1.962 9.1542×105 1.941 4.1991×104 1.905 0.215
    1/16 4.2969×105 1.974 2.3100×105 1.987 1.0669×104 1.977 0.697

     | Show Table
    DownLoad: CSV
    Table 3.  Numerical results obtained by REM-I for Example 6.1, s=0.01, T=1, (τ=1×104).
    h ME r LE rL2 HE rH1 CPU
    1/2 1.3515×104 6.7575×105 2.7862×104 0.127
    1/4 6.6217×106 4.351 4.4421×106 3.927 2.0899×105 3.737 0.282
    1/8 4.5813×107 3.853 2.9362×107 3.919 1.6356×106 3.676 0.890
    1/16 2.8925×108 3.985 1.6823×108 4.126 1.2135×107 3.753 3.120

     | Show Table
    DownLoad: CSV

    In Table 4, τ=h2 is taken. By Theorem 3.1, we have maxn1kn{uk(UE)k1} =O(τ2+h4x+h4y) =O(h4) as τ=h2. It is clearly observed from Table 4 that numerical solutions obtained by REM-I is convergent with an order of O(h4) in L2-, H1- and L-norms.

    Table 4.  Numerical results obtained by REM-I for Example 6.1, s=0.01, T=1, (τ=h2).
    h ME r LE rL2 HE rH1 CPU
    1/10 9.4174×107 - 4.8164×107 - 2.228×106 - 0.051
    1/20 5.9300×108 3.989 3.0126×108 3.999 1.4220×107 3.970 0.230
    1/30 1.1705×108 4.002 5.9652×109 3.994 2.8672×108 3.949 1.025
    1/40 3.5988×109 4.100 1.8360×109 4.096 9.0027×109 4.027 3.097

     | Show Table
    DownLoad: CSV

    In Table 5 and Table 6, τ=h2 is used. By Remak II and Theorem 4.1, we can infer maxn1knek1=O(h2) provided by the EFDM-II with τ=h2. By Remak III and Theorem 4.2, we can derive maxn1kn{uk(UE)k1}=O(h4) yielded by the REM-II with τ=h2. Table 5 and Table 6 clearly show that numerical solutions obtained by EFDM-II and REM-II possess the convergent orders of O(h2) and O(h4) in L2-, H1- and L-norms, respectively.

    Table 5.  Numerical results obtained by EFDM-II for Example 6.1, s=0.01, T=1, (τ=h2).
    h ME r LE rL2 HE rH1 CPU
    1/10 2.7125×103 1.467×103 6.751×103 0.0311
    1/20 6.8927×104 1.977 3.6976×104 1.988 1.7092×103 1.982 0.061
    1/30 3.0704×104 1.994 1.6457×104 1.997 7.6142×104 1.994 0.232
    1/40 1.7280×104 1.998 9.2619×105 1.998 4.2865×104 1.997 0.695

     | Show Table
    DownLoad: CSV
    Table 6.  Numerical results obtained by REM-II for Example 6.1, s=0.01, T=1, (τ=h2).
    h ME r LE rL2 HE rH1 CPU
    1/10 1.5748×106 1.0746×106 5.559×106 0.1500
    1/20 9.6861×108 4.023 6.4358×108 4.062 3.4267×107 4.020 1.292
    1/30 1.8432×108 4.092 1.2358×108 4.070 6.6824×108 4.032 6.175
    1/40 3.6389×109 5.640 2.6048×109 5.412 1.6076×108 4.953 19.494

     | Show Table
    DownLoad: CSV

    From Table 7, we can find that to obtain much more accurate solutions, REMs often cost much lower time. For example, to achieve ME1.0×108, 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.

    Table 7.  Numerical results obtained by REM-II for Example 6.1, s=0.01, T=1, (τ=h2).
    Methods EFDM-I REM-I EFDM-II REM-II
    (τ,h) (1800,1400) (1900,130) (11402,1140) (1900,130)
    ME 5.1880×108 1.1705×108 1.9212×105 1.8432×108
    CPU 31.872 1.025 48.804 6.175

     | Show Table
    DownLoad: CSV

    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.

    Figure 1.  Example 6.1: Numerical solution at T=1 obtained by EFDM-I with (τ,h) =(1.0×1004,116) (left), numerical solution at T=1 obtained by EFDM-II with (τ,h) =(1.0×1004,116) (middle) and exact solution at T=1 (right).

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

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

    uttΔu(x,t)=g(u(x,t),u(x,ts1),u(x,ts2),u(x,ts3))+f(x,t), (6.1)

    where s1=0.1875, s2=0.3125, s3=0.5625, s=max(s1,s2,s3),

    g(u(x,t),u(x,ts1),u(x,ts2),u(x,ts3))=u2(x,t)+[u(x,ts1)1]u(x,ts2)[u(x,ts3)+1],f(x,t)=[exp(x+y)sin(t)3]exp(x+y)sin(t)[exp(x+y)sin(ts1)1]exp(x+y)sin(ts2)[exp(x+y)sin(ts3)+1].

    IBVCs are determined by its exact solution u(x,t)=exp(x+y)sin(t).

    In Table 8, τ=h2 is used to satisfy the condition r2x+r2y<1. Table 8 shows that the generalized EFDM-I is second-order temporally accurate.

    Table 8.  Numerical results obtained by EFDM-I for Example 6.2, T=1, (τ=h2).
    h ME r LE rL2 HE rH1 CPU
    1/8 7.1666×104 3.8867×104 1.809×103 0.0192
    1/16 1.8363×104 1.965 9.8195×105 1.985 4.6060×104 1.974 0.024
    1/32 4.6021×105 1.997 2.4611×105 1.996 1.1567×104 1.994 0.055
    1/64 1.1531×105 1.997 6.1565×106 1.999 2.8950×105 1.998 0.319
    1/512 1.8021×107 2.000 9.6215×108 2.0000 4.5251×107 2.0000 256.052

     | Show Table
    DownLoad: CSV

    In Table 9 and Table 10, τ=1.0×104 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.

    Table 9.  Numerical results obtained by EFDM-I for Example 6.2, T=1, (τ=1×104).
    1/4 3.0678×103 1.6821×103 7.6089×103 0.1363
    1/8 8.1830×104 1.907 4.4351×104 1.923 2.0644×103 1.882 0.406
    1/16 2.0984×104 1.963 1.1218×104 1.983 5.2619×104 1.972 1.417
    1/32 5.2587×105 1.997 2.8120×105 1.996 1.3216×104 1.993 5.499

     | Show Table
    DownLoad: CSV
    Table 10.  Numerical results obtained by REM-I for Example 6.2, T=1, (τ=1×104).
    h ME r LE rL2 HE rH1 CPU
    1/4 4.194×105 2.7443×105 1.2633×104 0.520
    1/8 2.773×106 3.919 1.7571×106 3.965 9.2531×106 3.771 1.862
    1/16 1.707×107 4.022 1.0624×107 4.048 6.8458×107 3.757 7.240
    1/32 1.645×108 3.376 9.8264×109 3.435 6.3410×108 3.432 30.986

     | Show Table
    DownLoad: CSV

    By Remark I and Theorem 3.2, we can find that the generalized REM-I possess the convergent rate of O(h4) in H1-norm as hx=hy=h and τ=h2. Table 11 exactly confirms that numerical solutions are convergent with an order of O(h4) in L2-, L- and H1-norms as τ=h2 which satisfes r2x+r2y<1/4.

    Table 11.  Numerical solutions provided by REM-I for Example 6.2, T=1, (τ=h2).
    h ME r LE rL2 HE rH1 CPU
    1/4 6.2221×105 2.0943×105 1.5250×104 0.038
    1/8 4.0593×106 3.938 1.3489×106 3.957 1.0813×105 3.818 0.048
    1/16 2.3929×107 4.084 8.5942×108 3.972 7.4360×107 3.862 0.252
    1/32 1.5077×108 3.988 5.4254×109 3.986 5.6416×108 3.720 4.051

     | Show Table
    DownLoad: CSV

    By Remark II and Theorem 4.1, we can derive that the generalized EFDM-II possesses the convergent rate of O(h2) in H1-norm as hx=hy=h and τ=h2. By Remark III and Theorem 4.2, we can deduce that the generalized REM-II owns the convergent order of O(h4) in H1-norm as hx=hy=h and τ=h2. Table 12 and Table 13 exactly confirm that the generalized EFDM-II and the generalized REM-II have the convergent rates of O(h2) and O(h4) with respect to L2-, L- and H1-norms, respectively.

    Table 12.  Numerical results provided by EFDM-II for Example 6.2, T=1, (τ=h2).
    h ME r LE rL2 HE rH1 CPU
    1/4 3.8157×102 - 2.0972×102 - 9.4933×102 - 0.024
    1/8 1.0545×102 1.855 5.7124×103 1.876 2.6592×102 1.836 0.021
    1/16 2.7222×103 1.954 1.4551×103 1.973 6.8253×103 1.962 0.065
    1/32 6.8337×104 1.994 3.6541×104 1.994 1.7174×103 1.991 0.734

     | Show Table
    DownLoad: CSV
    Table 13.  Numerical results obtained by REM-II for Example 6.2, T=1, (τ=h2).
    h ME r LE rL2 HE rH1 CPU
    1/4 4.8076×104 1.4869×104 1.0973×103 0.080
    1/8 2.9946×105 4.005 8.8740×106 4.067 7.1303×105 3.944 0.156
    1/16 1.7963×106 4.059 5.4497×107 4.025 4.5193×106 3.980 1.147
    1/32 1.1331×107 3.987 3.3910×108 4.006 2.8498×107 3.987 17.646

     | Show Table
    DownLoad: CSV

    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.

    Table 14.  Comparisons of numerical results obtained by the current algorithms for Example 6.2, T=1.
    Methods EFDM-I REM-I EFDM-II REM-II
    (τ,h) (11024,1512) (11024,132) (11402,1140) (11024,132)
    ME 1.8021×107 2.3929×107 3.5806×105 1.1331×107
    CPU 256.05 0.252 252.273 17.646

     | Show Table
    DownLoad: CSV

    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.

    Figure 2.  Example 6.2: Numerical solutions at T=1 obtained by REM-I with (τ,h) =(1322,132) (left), numerical solution at T=1 obtained by REM-II with (τ,h) =(1322,132) (middle) and exact solution at T=1 (right).

    In a word, Tables 814 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 Ω=(0,1)×(0,1). To further compare EFDM-I with EFDM-II, they are used to solve wave equations with two variable delays as follows

    utta2Δu(x,t)=g(u(x,t),u(x,tν1(t)),u(x,tν2(t))+f(x,t),(x,t)Ω×(0,T]

    where

    g(u(x,t),u(x,tν1(t)),u(x,tν2(t))=βθdu(x,t)θd+ud(x,t)+bu(x,tν1(t))1+u2(x,tν2(t),
    f(x,t)=(12a2)exp(x+y+t)exp(x+y+t12|sin(t)|)1+exp(2x+2y+2t12(1cos(t)))92exp(x+y+t)9+exp(2x+2y+2t).

    Numerical results are exhibited in Table 15 and Table 16, respectively. The exact solution to this problem is u(x,t)=exp(x+y+t), by which we can get IBVCs. Here, we take the values of the parameters as follows: β=0.5, θ=3, d=2, b=1, ν1(t)=12|sin(t)|, and ν2(t)=14[1cos(t)]. From Table 15 and Table 16, we can obtain several conclusions as follows. (1) As stable condition r2x+r2y<1 is satisfied, numerical solutions obtained by the generalized EFDM-I with linear Lagrangian interpolation are convergent with an order of O(h2) in L2-, H1- and L-norms in the case of a=1 and τ=h2. Besides, as τ=h2, numerical solutions obtained by the generalized EFDM-II with linear Lagrangian interpolation converge to exact solution with an order of O(h2) 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 rx and ry, in contrast, r2x+r2y<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 r2x+r2y<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.

    Table 15.  Numerical results at T=1 provided by EFDM-I for Example 6.3, (τ=h2).
    a=1
    h ME r LE rL2 HE rH1 CPU
    1/4 6.4324×103 3.3163×103 1.5167×102 0.022
    1/8 1.7193×103 1.904 9.0800×104 1.869 4.2917×103 1.821 0.024
    1/16 4.2600×104 2.013 2.3171×104 1.970 1.1086×103 1.953 0.064
    1/32 1.0658×104 1.999 5.8214×105 1.993 2.7952×104 1.988 0.681
    a=10
    h ME r LE rL2 HE rH1 CPU
    1/4 7.7627×1020 3.8921×1020 1.7310×1021 0.021
    1/8 3.0460×1063 1.5232×1063 6.8955e×1063 0.024
    1/16 4.8082×104 - 1.4334×104 1.2615×103 0.064
    1/32 1.2264×104 1.971 3.5381×105 2.018 3.1772×104 1.989 0.677
    1/64 3.0881×105 1.990 8.8345×106 2.002 7.9602×105 1.997 10.822
    a=50
    h ME r LE rL2 HE rH1 CPU
    1/4 3.7278×1043 1.8725×1043 8.3312×1043 0.021
    1/8 1.1456×10151 5.7339×10150 2.5974×10151 0.024
    1/16 inf nan inf 0.063
    1/32 inf nan inf 0.713
    1/64 9.8788×10321 nan inf 10.812
    1/72 2.6889×105 1.2652×105 6.4906×105 17.489
    1/80 2.1797×105 1.993 1.0233×105 2.014 5.2583×105 1.998 26.321

     | Show Table
    DownLoad: CSV
    Table 16.  Numerical results at T=1 provided by EFDM-II, (τ=h2).
    a=1
    h ME r LE rL2 HE rH1 CPU
    1/2 2.6289×101 1.3144×101 5.4195×101 0.021
    1/4 7.5543×102 1.799 4.0166×102 1.710 1.7985×101 1.591 0.023
    1/8 1.9199×102 1.976 1.0250×102 1.970 4.8114×102 1.902 0.026
    1/16 4.7029×103 2.029 2.5659×103 1.998 1.2252×102 1.973 0.068
    1/32 1.1735×103 2.003 6.4144×104 2.000 3.0783×103 1.993 0.757
    a=10
    h ME r LE rL2 HE rH1 CPU
    1/2 7.7084×102 2.1323×102 1.7722×101 0.023
    1/4 2.2280×102 1.791 1.2451×102 0.776 5.8129×102 1.608 0.026
    1/8 5.3874×103 2.048 2.5637×103 2.28 1.2847×102 2.178 0.069
    1/16 1.3888×103 1.956 7.8055×104 1.716 3.6596×103 1.812 0.757
    1/64 3.4522×104 2.008 1.0268×104 2.926 8.8984×104 2.040 11.981
    a=50
    h ME r LE rL2 HE rH1 CPU
    1/2 7.5423×102 2.1693×102 1.7635×101 0.022
    1/4 2.0843×102 1.855 5.8594×103 1.888 5.3050×102 1.733 0.026
    1/8 5.6388×103 1.886 2.0030×103 1.549 1.4149×102 1.907 0.070
    1/16 1.4379×103 1.971 7.5796×104 1.402 3.6415×103 1.958 0.756
    1/64 3.5599×104 2.014 1.5314×104 2.307 9.1584×104 1.991 11.962
    1/72 2.8676×104 1.836 1.0786×104 2.976 7.0903×104 2.173 19.357

     | Show Table
    DownLoad: CSV

    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(τ2+h2x+h2y) in H1-norm as r2x+r2y<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(τ2+h4x+h4y) as 4(r2x+r2y)<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 (R2)ki,j (see Eq (4.3)), as τ=o(hx), τ=o(hy), and hx and hy tend to zero, the truncation error (R2)ki,j tends to zero. Although numerical solutions obtained by the second EFDM (4.4a)–(4.4c) are convergent with an order of O(τ2+h2x+h2y +(τhx)2+(τhy)2) in H1-norm under the conditions of τ=o(h32) (h=max(hx,hy)) and other grids conditions, it is suggested that hx=hy=h and τ=h2 are applied to obtain numerical solutions with an optimal convergent rate of O(h2) H1-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 rx and ry.

    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.

    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. 20202B 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.



    [1] S. A. Messaoudi, A. Fareh, N. Doudi, Well posedness and exponential stability in a wave equation with a strong damping and a strong delay, J. Math. Phys., 57 (2016), 111501. https://doi.org/10.1063/1.4966551 doi: 10.1063/1.4966551
    [2] G. Liu, H. Yue, H. Zhang, Long time behavior for a wave equation with time delay, Taiwan. J. Math., 21 (2017), 107–129. https://doi.org/10.1353/jqr.2017.0005 doi: 10.1353/jqr.2017.0005
    [3] M. Kafini, S. Messaoudi, Local existence and blow up of solutions to a logarithmic nonlinear wave equation with delay, Appl. Anal., 99 (2020), 530–547. https://doi.org/10.1080/00036811.2018.1504029 doi: 10.1080/00036811.2018.1504029
    [4] K. Zhu, Y. Xie, F. Zhou, Pullback attractors for a damped semilinear wave equation with delays, Acta Math. Sin., 34 (2018), 1131–1150. https://doi.org/10.1007/s10114-018-7420-3 doi: 10.1007/s10114-018-7420-3
    [5] Y. Wang, Pullback attractors for a damped wave equation with delays, Stoch. Dyn., 15 (2015), 1550003. https://doi.org/10.1142/S0219493715500033 doi: 10.1142/S0219493715500033
    [6] M. Jornet, Exact solution to a multidimensional wave equation with delay, Appl. Math. Comput., 409 (2021), 126421. https://doi.org/10.1016/j.amc.2021.126421 doi: 10.1016/j.amc.2021.126421
    [7] F. Rodríguez, M. Roales, A. Martín, Exact solutions and numerical approximations of mixed problems for the wave equation with delay, Appl. Math. Comput., 219 (2012), 3178–3186. https://doi.org/10.1016/j.amc.2012.09.050 doi: 10.1016/j.amc.2012.09.050
    [8] J. Z. Lobo, Y. S. Valaulikar, Group analysis of the one dimensional wave equation with delay, Appl. Math. Comput., 378 (2020), 125193. https://doi.org/10.1016/j.amc.2020.125193 doi: 10.1016/j.amc.2020.125193
    [9] J. K. Hale, Theory of Functional Differential Equations, New York: Spring-Verlag, 1977.
    [10] J. Wu, Theory and applications of partial functional differential equations, Berlin: Springer Science & Business Media, 1996.
    [11] M. C. Mackey, Unified hypothesis for the origin of aplastic anemia and periodic hematopoiesis, Blood, 51 (1978), 941–956. https://doi.org/10.1182/blood.V51.5.941.941 doi: 10.1182/blood.V51.5.941.941
    [12] Z. Ling, Z. Lin, Traveling wavefront in a Hematopoiesis model with time delay, Appl Math. Lett., 23 (2010), 426–431. https://doi.org/10.1016/j.aml.2009.11.011 doi: 10.1016/j.aml.2009.11.011
    [13] L. Berezansky, E. Braverman, Mackey-Glass equation with variable coefficients, Comput. Math. Appl., 51 (2006), 1–16.
    [14] W. S. C. Gurney, S. P. Blythe, R. M. Nisbet, Nicholson's blowflies revisited, Nature, 287 (1980), 17–21. https://doi.org/10.1038/287017a0 doi: 10.1038/287017a0
    [15] J. W. H. So, X. Zou., Traveling waves for the diffusive Nicholson's blowflies equation, Appl. Math. Comput., 122 (2001), 385–392. https://doi.org/10.1016/S0096-3003(00)00055-2 doi: 10.1016/S0096-3003(00)00055-2
    [16] M. R. S. Kulenovic, G. Ladas, Y. G. Sficas, Global attractivity in Nicholson's blowflies, Appl. Anal., 43 (1992), 109–124. https://doi.org/10.1016/0045-8732(92)90107-W doi: 10.1016/0045-8732(92)90107-W
    [17] G. E. Hutchinson, Circular causal systems in ecology, Ann. NY Acad. Sci, 50 (1948), 221–246. https://doi.org/10.1111/j.1749-6632.1948.tb39854.x doi: 10.1111/j.1749-6632.1948.tb39854.x
    [18] A. Yu. Kolesov, E. F. Mishchenkob, N. Kh. Rozov, A modification of Hutchinson's equation, Comp. Math. Math. Phys., 50 (2010), 1990–2002. https://doi.org/10.1134/S0965542510120031 doi: 10.1134/S0965542510120031
    [19] Y. N. Kyrychko, S. J. Hogan, On the use of delay equations in engineering applications, J VIB CONTROL, 16 (2010), 943–960. https://doi.org/10.1134/S0965542510120031 doi: 10.1134/S0965542510120031
    [20] Q. He, L. Kang, D. J. Evans, Convergence and stability of the finite difference scheme for nonlinear parabolic systems with time delay, Numer. Algor., 16 (1997), 129–153. https://doi.org/10.1023/A:1019130928606 doi: 10.1023/A:1019130928606
    [21] Z. Sun, Z. Zhang, A linearized compact difference scheme for a class of nonlinear delay partial differential equations, Appl. Math. Model., 37 (2013), 742–752. https://doi.org/10.1016/j.apm.2012.02.036 doi: 10.1016/j.apm.2012.02.036
    [22] D. Li, C. Zhang, J. Wen, A note on compact finite difference method for reaction–diffusion equations with delay, Appl. Math. Model., 39 (2015), 1749–1754. https://doi.org/10.1016/j.apm.2014.09.028 doi: 10.1016/j.apm.2014.09.028
    [23] C. Tang, C. Zhang, A fully discrete θ-method for solving semi-linear reaction–diffusion equations with time-variable delay, Math. Comput. Simulat., 179 (2021), 48–56. https://doi.org/10.1016/j.matcom.2020.07.019 doi: 10.1016/j.matcom.2020.07.019
    [24] A. V. Lekomtsev, V. G. Pimenov, Convergence of the alternating direction method for the numerical solution of a heat conduction equation with delay, Proc. Steklov Inst. Math., 272 (2011), 101–118.
    [25] Q. Zhang, C. Zhang, L. Wang, The compact and Crank–Nicolson ADI schemes for two-dimensional semilinear multidelay parabolic equations, J. Comput. Appl. Math., 306 (2016), 217–230.
    [26] D. Deng, The study of a fourth-order multistep ADI method applied to nonlinear delay reaction-diffusion equations, Appl. Numer. Math., 96 (2015), 118–133. https://doi.org/10.1016/j.apnum.2015.05.007 doi: 10.1016/j.apnum.2015.05.007
    [27] Q. Zhang, D. Li, C. Zhang, D. Xu, Multistep finite difference schemes for the variable coefficient delay parabolic equations, J. Differ. Equ. Appl., 22 (2016), 745–765. https://doi.org/10.1080/10236198.2016.1142539 doi: 10.1080/10236198.2016.1142539
    [28] J. Xie, Z. Zhang, The high-order multistep ADI solver for two-dimensional nonlinear delayed reaction–diffusion equations with variable coefficients, Comput. Math. Appl., 75 (2018), 3558–3570.
    [29] D. Li, C. Zhang, H. Qin, LDG method for reaction–diffusion dynamical systems with time delay, Appl. Math. Comput., 217 (2011), 9173–9181. https://doi.org/10.1016/j.amc.2011.03.153 doi: 10.1016/j.amc.2011.03.153
    [30] Z. Jackiewicz, B. Zubik-Kowal, Spectral collocation and waveform relaxation methods for nonlinear delay partial differential equations, Appl. Numer. Math., 56 (2006), 433–443. https://doi.org/10.1016/j.apnum.2005.04.021 doi: 10.1016/j.apnum.2005.04.021
    [31] R. Burger, R. Ruiz-Baier, C. Tian, Stability analysis and finite volume element discretization for delay-driven spatio-temporal patterns in a predator–prey model, Math. Comput. Simulat., 132 (2017), 28–52. https://doi.org/10.1016/S0737-0806(17)30375-1 doi: 10.1016/S0737-0806(17)30375-1
    [32] D. Deng, J. Xie, Y. Jiang and D. Liang, A second-order box solver for nonlinear delayed convection-diffusion equations with Neumann boundary conditions, Int. J. Comput Math., 96 (2019), 1879–1898.
    [33] Q. Zhang, C. Zhang, A new linearized compact multisplitting scheme for the nonlinear convection-reaction-diffusion equations with delay, Commun. Nonlinear Sci. Numer. Simulat. 18 (2013), 3278–3288. https://doi.org/10.1016/j.cnsns.2013.05.018
    [34] B. Liu, C. Zhang, A spectral Galerkin method for nonlinear delay convection-diffusion-reaction equations, Comput. Math. Appl., 69 (2015), 709–724. https://doi.org/10.1016/j.camwa.2015.02.027 doi: 10.1016/j.camwa.2015.02.027
    [35] M. A. Castro, F. Rodríguez, J. Cabrera, J. A. Martín, Difference schemes for time-dependent heat conduction models with delay, Int. J. Comput. Math., 91 (2014), 53–61.
    [36] Q. Zhang, C. Zhang, A compact difference scheme combined with extrapolation techniques for solving a class of neutral delay parabolic differential equations, Appl. Math. Lett., 26 (2013), 306–312.
    [37] C. Zhang, W. Wang, B. Liu, T. Qin, Multi-domain Legendre spectral collocation method for nonlinear neutral equations with piecewise continuous argument, Int. J. Comput. Math., 95 (2018), 2419–2432.
    [38] H. Liang, Convergence and asymptotic stability of Galerkin methods for linear parabolic equations with delays, Appl. Math. Comput., 264 (2015), 160–178. https://doi.org/10.1016/j.amc.2015.04.104 doi: 10.1016/j.amc.2015.04.104
    [39] C. Zhang, Z. Tan, Linearized compact difference methods combined with Richardson extrapolation for nonlinear delay Sobolev equations, Commun. Nonlinear Sci. Numer. Simulat., 91 (2020), 105461. https://doi.org/10.1016/j.cnsns.2020.105461 doi: 10.1016/j.cnsns.2020.105461
    [40] A. B. Chiyaneh, H. Duru, On adaptive mesh for the initial boundary value singularly perturbed delay Sobolev problems, Numer. Meth. Part. D. E., 36 (2019), 228–248. https://doi.org/10.1002/num.22417 doi: 10.1002/num.22417
    [41] A. B. Chiyaneh, H. Duru, Uniform difference method for singularly pertubated delay Sobolev problems, Quaest. Math., 43 (2020), 1713–1736. https://doi.org/10.2989/16073606.2019.1653395 doi: 10.2989/16073606.2019.1653395
    [42] Q. Zhang, C. Zhang, D. Deng, A compact difference scheme and Richardson extrapolation algorithm for solving a class of the nonlinear delay hyperbolic partial differential equations (in Chinese), J. Numer. Meth. Comput. Appl., 34 (2013), 167–176.
    [43] Q. Zhang C. Zhang, D. Deng, Compact alternating direction implicit method to solve two-dimensional nonlinear delay hyperbolic differential equations, Int. J. Comput. Math., 91 (2014), 964–982. https://doi.org/10.1080/00207160.2013.810216 doi: 10.1080/00207160.2013.810216
    [44] E. C. Du Fort, S. P. Frankel, Conditions in the numerical treatment of parabolic differential equations, Math Tables Other Aids Comput., 7 (1953), 135–152. https://doi.org/10.2307/2002754 doi: 10.2307/2002754
    [45] Z. Sun, Numerical methods for partial differential equations (In Chinese), Beijing: Science Press, 2012.
    [46] D. Deng, D. Liang, The time fourth-order compact ADI methods for solving two-dimensional nonlinear wave equations, Appl. Math. Comput., 329 (2018), 188–209. https://doi.org/10.1016/j.amc.2018.02.010 doi: 10.1016/j.amc.2018.02.010
    [47] D. Deng, Q. Wang, A class of weighted energy-preserving Du Fort-Frankel difference schemes for solving sine-Gordon-type equations, Commun. Nonlinear Sci. Numer. Simulat., 117 (2023), 106916. https://doi.org/10.1016/j.cnsns.2022.106916 doi: 10.1016/j.cnsns.2022.106916
    [48] D. Deng, Z. Li, High-order structure-preserving Du Fort-Frankel schemes and their analyses for the nonlinear Schrödinger equation with wave operator, J. Comput. Appl. Math., 417 (2023), 114616. https://doi.org/10.1016/j.cam.2022.114616 doi: 10.1016/j.cam.2022.114616
  • This article has been cited by:

    1. Dingwen Deng, Zhu-an Wang, Zilin Zhao, The maximum norm error estimate and Richardson extrapolation methods of a second-order box scheme for a hyperbolic-difference equation with shifts, 2024, 1607-3606, 1, 10.2989/16073606.2024.2385424
    2. V. G. Pimenov, A. B. Lozhnikov, Richardson Method for a Diffusion Equation with Functional Delay, 2023, 321, 0081-5438, S204, 10.1134/S0081543823030173
  • Reader Comments
  • © 2023 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(1526) PDF downloads(85) Cited by(2)

Figures and Tables

Figures(2)  /  Tables(16)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog