1.
Introduction
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]
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,t−s),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
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≥0. Parameters α, β, θ, γ and d∈(1,∞) are positives constants. Besides, as
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
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.
2.
Partition and Notations
To begin with, we divide the spatial domain Ω= [b1,b2] ×[d1,d2] ⊂ R2, where b1, b2, d1 and d2 are constant numbers. Let hx=(b2−b1)/m1, hy=(d2−d1)/m2, (m1,m2∈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 τ (τ=s/n1=T/n) (n1,n∈Z+) is used. If the restricted grid is not used, Lagrange interpolation formula will be applied to approximate the delay term u(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 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
Besides, we introduce
On Ωhτ, further define
Then, on Ωh×Ωτ, we further introduce the following difference notations
Finally, ∀ U, V ∈S0h, 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 uki,j=u(xi,j,tk), xi,j∈Ωh, −n1≤k≤n. 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
where
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
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(μ, υ, 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
Moreover, denote eki,j=uki,j−Uki,j, 0≤i≤m1, 0≤j≤m2, −n1≤k≤n and h=max{hx,hy}, and give several lemmas used later.
Lemma 3.1 [45] Let V∈S0h, we have
Lemma 3.2 [26] Let both A and B be positive constants, {Fk|k≥0} be a non-negative sequence and satisfy the following inequality Fk+1≤A+BτK∑k=0Fk. Then we can obtain max0≤k≤K+1Fk≤Aexp(B(K+1)τ). Moreover, if Fk+1≤A+BτK+1∑k=0Fk, we have max0≤k≤K+1Fk≤Aexp(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+a2⟨ek,ek+1⟩1,x+a2⟨ek,ek+1⟩1,y, c2=1−r2x−r2y. Then when r2x+r2y<1, the following inequalities
hold.
Proof. By using the equality αβ=[(α+β)2−(α−β)2]/4 and simple computations, we have
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
Multiplying a2hxhy to both sides of Eq (3.5) and summing i from 1 to m1−1 and j from 1 to m2−1, we obtain
which shows that
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
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(x,t)∈C4,4(Ω×[−s,T]). Then there exists positive constant c3 such that
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 μh≤hx ≤h, μh≤ hy≤h and τ≤σh12+ϵ, we can derive that
under the conditions that
and Assumption A are valid. Here,
Proof. For convenience, denote f(uki,j)=f(uki,j,uk−n1i,j,xi,yj,tk),f(Uki,j)=f(Uki,j,Uk−n1i,j,xi,yj,tk) without confusion. Subtracting Eq (3.3a) from Eq (3.1), the error equations are derived as follows
As eki,j=0 (−n1≤k≤0), it is obvious that Eq (3.10) holds for −n1≤k≤0. Assuming that Eq (3.10) is valid for −n1≤k≤l, we will prove that it is also true for k=l+1.
As μh≤hx≤h and μh≤hy≤h, it is easy to find that
Thus, using τ≤σh12+ϵ and Eq (3.12), we obtain
Besides, applying h≤min{((ε0μ2σ2c4) 12ϵ, ε0√μ4c4)}, 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δˆtek gives that
In what follows, we estimate each terms of Eq (3.16) step by step. Using the equality α2−β2 =(α−β)(α+β) 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αβ≤1εα2+εβ2 and (α+β)2≤2(α2+β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 γ and summing γ from 0 to k, we have
Hence, applying Lemma 3.2 to Eq (3.23), it follows from τ{4+8c23(b2−b1)2(d2−d1)23c22a2[(b2−b1)2+(d2−d1)2]}≤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 τ, 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
Then, as 4a2[(τhx)2+(τhy)2]<1, and grid parameters τ, hx and hy are small enough, we obtain
Proof. Let r1(x,t)=a212∂4u∂x4(x,t) and r2(x,t)=a212∂4u∂y4(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 v1(x,t) and v2(x,t) be the exact solutions of the following initial-boundary-value problems (IBVPs)
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
respectively.
Denote pki,j=eki,j+h2x(v1)ki,j+h2y(v2)ki,j, 0≤i≤m1, 0≤j≤m2, −n1≤k≤n. 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
where
and
In Eq (3.33a), the following Taylor expansion formula
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
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 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.
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 (R2)ki,j and replacing uki,j by its approximation Uki,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 Gk=(1+r2x+r2y)‖δtek+12‖2+a2⟨ek,ek+1⟩1,x+a2⟨ek,ek+1⟩1,y. Then the following inequalities
are valid.
Proof. Using Lemma 3.1, we obtain
Applying Eqs (4.6a), (4.6b) and the equality αβ= (α+β)2/4 −(α−β)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(x,t) ∈C4,4(Ω×[−s,T]), there exists positive constant c6 such that
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 μh≤hx≤h, μh≤hy≤h, τ≤σh32+ϵ. It holds that
under the conditions that
and Assumption A hold. Here,
Proof. Also, let f(uki,j)=f(uki,j,uk−n1i,j,xi,j,tk), f(Uki,j)=f(Uki,j,Uk−n1i,j,xi,j,tk) for convenience. Subtracting Eq (4.4a) from Eq (4.2), the error equations are derived as follows
As eki,j=0 (−n1≤k≤0), it is obvious that Eq (4.9) holds for −n1≤k≤0. Assuming that Eq (4.9) is valid for −n1≤k≤l, we will prove that it is also true for k=l+1.
As μh≤hx≤h and μh≤hy≤h, it is easy to find that
Thus, using τ≤σh32+ϵ and Eq (4.11), we obtain
Besides, applying h≤min((με03c7σ2)12(1+ϵ),√με06c7,(μ3ε06c7σ2)12ϵ), 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δˆtek 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αβ≤1εα2+εβ2 and (α+β)2≤2(α2+β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 γ, summing γ from 0 to k, we can arrive at
As τ{4+4c21(b2−b1)2(d2−d1)2(1+r2x+r2y)3[(b2−b1)2+(d2−d1)2]a2}≤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(τ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.
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(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
Then under the conditions of Theorem 4.1, hx=O(hy), τhx and τhy ⟶0 as grid parameters τ, hx, hy ⟶0, we arrive at
as long as parameter grids τ hx and hy are small enough.
Proof. Let w1(x,t)=112∂4u∂t4(x,t), w2(x,t)=a212∂4u∂x4(x,t), w3(x,t)=a212∂4u∂y4(x,t) and w4(x,t)=a2∂2u∂t2(x,t). Then, we define the grid functions (wl)ki,j =wl(xi,j,tk), (l=1,2,3,4), xi,j∈ˉΩh, −n1≤k≤n. Thus, it follows from Eq (4.3) that
Furthermore, the error Eqs (4.10a)–(4.10c) can be rewritten as
Suppose that functions v1(x,t), v2(x,t), v3(x,t) and v4(x,t) are exact solutions of the following IBVPs
respectively, where
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
Denote qki,j=eki,j+h2x(v1)ki,j+h2y(v2)ki,j+(τ2h2x+τ2h2y)(v3)ki,j+τ2(v4)ki,j, 0≤i≤m1, 0≤j≤m2, −n1≤k≤n. 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
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 αβ≤(α2+β2)/2 and max((hyhx)2,(hxhy)2)≤μ−2, we have
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
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)k‖1=O(h4), 0≤k≤n 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.
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 τ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
where sl(t)>0 (l=1,2,…,L) and s=max1≤l≤Lsl(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,(nl∈Z+) is adopted to make t=Ksl (K∈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(Uki,j, Uk−n1i,j, xi,j,tk) with f(Uki,j, Uk−n1i,j, Uk−n2i,j, …, Uk−nLi,j, xi,j, tk). Also, corresponding theoretical results are derived by using similar analytical methods as long as f(u(x,t), u(x,t−s1), u(x,t−s2), …, u(x,t−sL), 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τ, (Pl∈Z+ 0≤θl≤1). Then delay term u(x,t−sl(tk)) can be approximated by applying the linear Lagrange interpolation as follows u(xi,j,t−sl(tk)) =(1−θl)uk−Pli,j +θluk−Pl−1i,j +O(τ2). Write ˆUk,li,j= =(1−θl)Uk−Pli,j +θlUk−Pl−1i,j, which is viewed as a second-order approximations to u(xi,j,t−sl(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,Uk−n1i,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.
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., 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=max0≤k≤n‖ek‖∞, LE=max0≤k≤n||ek||2 and HE=max0≤k≤n‖ek‖21 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
with Dirichlet initial-boundary conditions by using our numerical methods. Here, f(x,t)=−(sinx+cosy)2(sint)2(sint−s)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.
In Table 2 and Table 3, timestep τ=1×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, τ=h2 is taken. By Theorem 3.1, we have max−n1≤k≤n{‖uk−(UE)k‖1} =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.
In Table 5 and Table 6, τ=h2 is used. By Remak II and Theorem 4.1, we can infer max−n1≤k≤n‖ek‖1=O(h2) provided by the EFDM-II with τ=h2. By Remak III and Theorem 4.2, we can derive max−n1≤k≤n{‖uk−(UE)k‖1}=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.
From Table 7, we can find that to obtain much more accurate solutions, REMs often cost much lower time. For example, to achieve ME≈1.0×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 Ω=(0,1)×(0,1). On Ω×[−s,T], we consider numerical solutions of the following wave equations with three constant delays
where s1=0.1875, s2=0.3125, s3=0.5625, s=max(s1,s2,s3),
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.
In Table 9 and Table 10, τ=1.0×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(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.
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 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 Ω=(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
where
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[1−cos(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.
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(τ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.
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. 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.