Optimal Rate of Convergence for a Nonstandard Finite Difference Galerkin Method Applied to Wave Equation Problems

Abstract

The optimal rate of convergence of the wave equation in both the energy and the 𝐿²-norms using continuous Galerkin method is well known. We exploit this technique and design a fully discrete scheme consisting of coupling the nonstandard finite difference method in the time and the continuous Galerkin method in the space variables. We show that, for sufficiently smooth solution, the maximal error in the 𝐿²-norm possesses the optimal rate of convergence 𝑂(ℎ²+(Δ𝑡)²) where ℎ is the mesh size and Δ𝑡 is the time step size. Furthermore, we show that this scheme replicates the properties of the exact solution of the wave equation. Some numerical experiments should be performed to support our theoretical analysis.

1. Introduction

Most physical phenomena such as the acoustics, electromagnetic, and elastic problems are modeled by the wave equation. The qualitative solution of the model in the spacetime domain is always a very delicate but a fundamental issue that needs careful study. Our point of departure of this paper is to consider the following model of the wave equation: find 𝑢(𝑥,𝑡) such that 𝜕²𝑢𝜕𝑡²()(),𝑢()−Δ𝑢=𝑓in𝑥,𝑡∈Ω×0,𝑇𝑥,0=𝑢₀(𝑥),𝜕𝑢(𝑥,0)𝜕𝑡=𝑢₁(𝑥)𝑥∈Ω, (1)(2)(3) where Ω is a smooth bounded domain in ℝ² with smooth boundary 𝜕Ω. Problem (1)–(3) consists of constant coefficients and 𝑓(𝑥,𝑡) the source term, 𝑢₀(𝑥) and 𝑢₁(𝑥) are prescribed as the initial data. Furthermore, (0,𝑇) is taken to be a finite time interval and the boundary conditions satisfied by 𝑢(𝑥,𝑡) are given by 𝑢(𝑥,𝑡)=0,(𝑥,𝑡)∈𝜕Ω×(0,𝑇). (4) The methods which have been heavily used for the study of the wave equation (1)–(4) in physical life are the continuous as well as the discontinuous Galerkin methods; see [1, 2] for more details. The reason for these methods may be due to the way they deal with heterogeneous media and arbitrary shaped geometric objects, represented by unstructured grids.

The advantages of the continuous Galerkin method are enormous. Firstly, the convergence theory of this method is based on lower regularity (differentiability) requirements than the finite difference and the spectral methods. Secondly, the method retains the important energy conservation properties provided by the discrete version of the initial/boundary valued problem such as the one under consideration. Thirdly, the computation and the analysis from the Galerkin method could be extended in the approximation of the nonlinear wave equation. Furthermore, the method could be applicable to problems of any desired order of accuracy. For more details on these advantages, see in [3–5].

The apriori error estimate for continuous Galerkin approximation of the wave equation (1)–(4) was first derived by Dupont [6] and later improved by Baker [7], both for continuous and discrete time schemes. Gekeler [8] analyzed general multistep methods for the time discretization of second-order hyperbolic equation, when a Galerkin procedure is used in space. The nonclassical finite element treatment of the wave equation can be seen in Johnson [9] and Richter [10].

In this paper, we exploit and present a reliable technique consisting of coupling the nonstandard finite difference (NSFD) method in time and the continuous Galerkin (CG) method in the space variables. A similar approach was done for the first time using parabolic problems more specifically the diffusion equations in the nonsmooth domain as seen in [11]. The NSFD method was initiated by Mickens in [12] and major contributions to the foundation of the NSFD method could be seen in [13, 14]. Since its initiation, the NSFD method has been extensively applied to many concrete problems in engineering and science; see [12, 15, 16] for an overview. This paper compliments the technique used in [11]. The technique is geared toward obtaining a sufficiently smooth solution, the maximal error in the 𝐿²-norm, and to show that the error across the entire interval convergences optimally as 𝑂(ℎ²+Δ𝑡²) where ℎ is the mesh size and Δ𝑡 is the time step size. The reliability of this technique comes from the fact that the NSFD-CG method preserves both the energy features and the hyperbolicity of the exact solution of the wave equation (1)–(4).

The organization of the paper is as follows. Under Section 2, we review some of the useful spaces and their notations needed in the paper. In Section 3, we gather essential tools necessary to prove the main result of the paper. We present, in Section 4, a reliable scheme NSFD-CG and show that the numerical solution obtained from this scheme attains the optimal convergence rate in the energy as well as in the 𝐿²-norms. Furthermore, we show that the scheme under consideration replicates the properties of the exact solution. Section 5 is devoted to some numerical experiments using a numerical example which confirms the optimal rate of convergence of the solution proved analytically in Section 4. The concluding remarks are given in Section 6 and these underline how the work fits in the existing literature and also how it can be extended for further work.

2. Notations

In this section, we will review some of the spaces which we will be using in the paper together with their notations and possibly properties. For 𝑠≥0, 𝐻^𝑠(Ω) will denote the Sobolev space of real-valued functions on Ω, and the norm on 𝐻^𝑠(Ω) will be denoted by ‖⋅‖_𝑠. See [17] for the definitions and the relevant properties of these spaces. In a particular case, where 𝑠=0 the space 𝐻⁰(Ω)=𝐿²(Ω) and its inner product together with the norm will be denoted, respectively, by (𝑢,𝑣)=_Ω𝑢𝑣𝑑𝑥,𝑢,𝑣∈𝐿²(Ω),‖𝑢‖_𝐿²_(Ω)={(𝑢,𝑢)}^1/2,𝑢∈𝐿²(Ω). (5) In addition, 𝐶^∞₀(Ω) will denote the space of infinitely differentiable functions with support compactly contained in Ω. The space 𝐻¹₀(Ω) will denote the subspace of 𝐻¹(Ω) obtained by completing 𝐶^∞₀(Ω) with respect to the norm ‖⋅‖₁. Following [17], for 𝑋 a Hilbert space, we will more generally use the Sobolev space 𝐻^𝑠[(0,𝑇);𝑋], where 𝑠≥0 and in the case where 𝑠=0 we will have 𝐻⁰[(0,𝑇);𝑋]≡𝐿²[(0,𝑇);𝑋] with norm ‖𝑣‖_𝐿²_{[(0,𝑇);𝑋]}=^𝑇₀‖𝑣(⋅,𝑡)‖_𝑋𝑑𝑡^1/2. (6) In practice, 𝑋 will be the Sobolev space 𝐻^𝑚(Ω) or 𝐻^𝑚₀(Ω). Associated with (1) is the bilinear form 𝑎(𝑢,𝑣)=_Ω∇𝑢∇𝑣𝑑𝑥,𝑢,𝑣∈𝐻¹(Ω). (7)𝑎(⋅,⋅) will be symmetric and positive definite; that is, 𝑎(𝑢,𝑣)=𝑎(𝑣,𝑢),𝑎(𝑢,𝑢)≥0. (8)

3. The Continuous Galerkin Method

Having the previously mentioned notations in place, we proceed under this section to gather essential tools necessary to prove the main result of our paper. We begin by stating the following weak problem of (1)–(4).

Find 𝑢∈𝐿²[(0,𝑇);𝐻¹₀(Ω)] given such that 𝜕²𝑢(⋅,𝑡)𝜕𝑡²,𝑣−𝑎(𝑢(⋅,𝑡),𝑣)=𝑓(⋅,𝑡),𝑣∀𝑣∈𝐻¹₀(Ω)(𝑢())=𝑢,𝑡≥0,⋅,0,𝑣₀,(),𝑣𝜕𝑢⋅,0=𝑢𝜕𝑡,𝑣₁,𝑣. (9) See [8] for the existence and the uniqueness of a solution 𝑢 of (9). Hence forth, in appropriate places to follow, additional conditions on the regularity of 𝑢 which guarantee the convergence results will be imposed. We continue next by providing the framework for stating the discrete version of (9). To this end, we let 𝒯_ℎ be a regular family of triangulations of Ω consisting of compatible triangles 𝑇 of diameter ℎ_𝑇≤ℎ; see [17] for more. For each mesh size 𝒯_ℎ, we associate the finite element space 𝑉_ℎ of continuous piecewise linear function that are zero on the boundary 𝑉_ℎ∶=𝑣_ℎ∈𝐶⁰Ω;𝑣_ℎ|_𝜕Ω=0,𝑣_ℎ|_𝑇∈𝑃₁,∀𝑇∈𝒯_ℎ, (10) where 𝑃₁ is the space of polynomials of degree less than or equal to 1 and 𝑉_ℎ is a finite dimensional subspace of the Sobolev space 𝐻¹₀(Ω). It is well known that 𝑉_ℎ parametrized by ℎ∈(0,1) possesses the following approximation properties: there exists a constant 𝐶 such that if 𝑣∈𝐻¹₀(Ω)∩𝐻²(Ω) we have inf_{𝜒∈𝑉ℎ}𝑣−𝜒+ℎ𝑣−𝜒₁≤𝐶ℎ²‖𝑣‖₂. (11) By the use of the energy method and Gronwall’s Lemma, there exists a discrete Galerkin solution 𝑢_ℎ∈𝑉_ℎ such that 𝜕²𝑢_ℎ𝜕𝑡²,𝑣_ℎ−𝑎𝑢_ℎ,𝑣_ℎ=𝑓,𝑣_ℎ,∀𝑣_ℎ∈𝑉_ℎ[],𝑢,𝑡∈0,𝑇_ℎ,𝑣_ℎ=𝑃_ℎ𝑢₀,𝑣_ℎ,𝜕𝑢_ℎ𝜕𝑡,𝑣_ℎ=𝑃_ℎ𝑢₁,𝑣_ℎ, (12) where 𝑃_ℎ denote the 𝐿²-projection onto 𝑉_ℎ. Furthermore, we let 𝑢∈𝐻²(Ω) and we define the Galerkin projection 𝜋_ℎ𝑢 of 𝑢 in 𝑉_ℎ by requiring that 𝑎𝜋_ℎ𝑢,𝑣=𝑎(𝑢,𝑣),∀𝑣∈𝑉_ℎ,𝜕^𝑘𝜋_ℎ𝑢𝜕𝑡^𝑘=𝜋_ℎ𝜕^𝑘𝑢𝜕𝑡^𝑘,𝑘=0,1,2,𝑡∈[0,𝑇]. (13)(14) The previous essential tools lead to the immediate consequence of the approximation properties (11).

Lemma 1. Let 𝑢 be the solution of (9). Then, there exists a unique mapping 𝜋_ℎ𝑢∈𝐿²[(0,𝑇);𝑉_ℎ] which satisfies (13). Furthermore, if for some integer 𝑘≥0, 𝜕^𝑘𝑢/𝜕𝑡^𝑘∈𝐿^𝑝[(0,𝑇);𝐻²(Ω)], then 𝜕^𝑘𝜋_ℎ𝑢𝜕𝑡^𝑘∈𝐿^𝑝(0,𝑇);𝑉_ℎ,𝜕𝜕𝑡^𝑘𝑢−𝜋_ℎ𝑢_{𝐿𝑝[(0,𝑇);𝐿2(Ω)]}≤𝐶ℎ²𝜕𝜕𝑡^𝑘𝑢_{𝐿𝑝[(0,𝑇);𝐻2(Ω)]} (15) for some constant 𝐶 independent of the mesh size.

4. Coupled Nonstandard Finite Difference and Continuous Galerkin Methods

Instead of the continuous Galerkin method summarized previously, we present in this section a reliable scheme NSFD-CG consisting of the nonstandard finite difference method in time and the continuous Galerkin method in the space variable. We show that the numerical solution obtained from the scheme NSFD-CG attained the optimal convergence in both the energy and 𝐿²-norms. We proceed, in this regard, by letting the step size 𝑡_𝑛=𝑛Δ𝑡 for 𝑛=0,1,2,…,𝑁. For a sufficiently smooth function 𝑣(𝑥,𝑡), we set 𝜕𝜕𝑡^𝑘𝑣^𝑛=𝜕𝜕𝑡^𝑘𝑣⋅,𝑡_𝑛,𝑣^𝑛=𝑣⋅,𝑡_𝑛,𝑘≥0. (16) With this, we proceed to find the NSFD-CG approximation {𝑈^𝑛_ℎ} such that 𝑈^𝑛_ℎ≈𝑢^𝑛_ℎ at discrete time 𝑡_𝑛. That is, find a sequence {𝑈^𝑛_ℎ}^𝑁_𝑛=0 in 𝑉_ℎ such that 𝛿²𝑈^𝑛_ℎ,𝑣_ℎ+𝑎𝑈^𝑛_ℎ,𝑣_ℎ=𝑓^𝑛,𝑣_ℎ,∀𝑣_ℎ∈𝑉_ℎ,𝑛=1,2,…,𝑁−1, (17) where 𝛿²𝑈^𝑛_ℎ=𝑈_ℎ^𝑛+1−2𝑈^𝑛_ℎ+𝑈_ℎ^𝑛−1𝜓(Δ𝑡)²,𝑛=1,2,…,𝑁−1 (18) and 𝜓(Δ𝑡)=2sin(Δ𝑡/2) restricted between 0<𝜓(Δ𝑡)<1. The initial conditions 𝑈⁰_ℎ∈𝑉_ℎ and 𝑈¹_ℎ∈𝑉_ℎ are given by 𝑈⁰_ℎ,𝑣_ℎ=𝑃_ℎ𝑢₀,𝑣_ℎ𝑈_1ℎ,,𝑣_ℎ=𝑈⁰_ℎ+𝜓(Δ𝑡)𝑃_ℎ𝑢₁+𝜓(Δ𝑡)²2𝑈⁰_ℎ,𝑣_ℎ,∀𝑣_ℎ∈𝑉_ℎ, (19)(20) where 𝑈⁰_ℎ∈𝑉_ℎ is defined by 𝑈⁰_ℎ,𝑣_ℎ=𝑓⁰,𝑣_ℎ−𝑎𝑢₀,𝑣_ℎ,∀𝑣_ℎ∈𝑉_ℎ. (21) If 𝑓=0 in (1), we will have in view of (17) an exact scheme 𝑈_ℎ^𝑛+1−2𝑈^𝑛_ℎ+𝑈_ℎ^𝑛−14sin²(Δ𝑡/2),𝑣_ℎ+∇𝑈^𝑛_ℎ,∇𝑣_ℎ=0, (22) which according to Mickens [12] replicates both the energy preserving features and the properties of the exact solution (1)–(4). In order to state and prove the main result, we need a framework on which this result is based. To this end, we proceed by decomposing the error denoted by 𝑒^𝑛=𝑢^𝑛−𝑈^𝑛_ℎ into the following error equation: 𝑒^𝑛=𝑢^𝑛−𝑤^𝑛_ℎ+𝑤^𝑛_ℎ−𝑈^𝑛_ℎ=𝜂^𝑛−𝜙^𝑛, (23) where 𝑤^𝑛_ℎ=𝜋_ℎ𝑢^𝑛∈𝑉_ℎ is the Galerkin projector of 𝑢^𝑛. Due to the regularity assumption mentioned earlier, the exact solution 𝑢 of (1)–(4) satisfies 𝜕²𝑢^𝑛𝜕𝑡²,𝑣_ℎ+𝑎𝑢^𝑛,𝑣_ℎ=𝑓^𝑛,𝑣_ℎ,∀𝑣_ℎ∈𝑉_ℎ,𝑛=1,2,…,𝑁. (24) Subtracting (17) from (24) and using some properties of the Galerkin projection in the space we have 𝛿²𝑤^𝑛_ℎ−𝛿²𝑈^𝑛_ℎ,𝑣_ℎ+𝑎𝑤^𝑛_ℎ−𝑢^𝑛,𝑣_ℎ=𝛿²𝑤^𝑛_ℎ−𝜕²𝑢^𝑛𝜕𝑡²,𝑣_ℎ, (25) from where we obtain in view of (23) 𝛿²𝜙^𝑛,𝑣_ℎ+𝜙^𝑛,𝑣_ℎ=𝑟^𝑛,𝑣_ℎ,∀𝑣_ℎ∈𝑉_ℎ,𝑛=1,2,…,𝑁−1, (26) where 𝑟^𝑛=𝛿²𝑤^𝑛_ℎ−𝜕²𝑢^𝑛/𝜕𝑡². In view of the necessity for error bound, we set 𝑟^𝑛=𝛿²𝑤^𝑛_ℎ−𝜕²𝑢^𝑛𝜕𝑡²𝜙,for𝑛=1,2,…,𝑁−1₁−𝜙⁰𝜓(Δ𝑡)²for𝑛=0, (27) from where we define 𝑅^𝑛=𝜓(Δ𝑡)^𝑛_𝑚=0𝑟^𝑚. (28) The previously mentioned framework leads to the following main result.

Theorem 2. Let 𝑢 be the solution of (9) and {𝑈^𝑛_ℎ}^𝑁_𝑛=0 in 𝑉_ℎ a sequence defined by (17)–(21). Suppose that 𝑢∈𝐿²[(0,𝑇);𝐻²(Ω)], 𝜕𝑢/𝜕𝑡∈𝐿²[(0,𝑇);𝐻²(Ω)], and 𝜕^𝑘𝑢/𝜕𝑡^𝑘∈𝐿²[(0,𝑇);𝐿²(Ω)] for 𝑘=3,4. Then, there exists a constant 𝐶>0 independent of Δ𝑡 and the mesh refinement size ℎ: max_{0≤𝑛≤𝑁}𝑢(⋅,Δ𝑡)−𝑈^𝑛_ℎ0≤𝐶ℎ²+(Δ𝑡)². (29) Furthermore, the discrete solution replicates all the properties of solution of the hyperbolic equation in the limiting case of the space independent equation.

We prove the previous Theorem 2 thanks to the following series of results.

Proposition 3. The following result holds for a constant 𝐶>0 that is independent of Δ𝑡 and the mesh refinement size ℎ: max_{𝑛=0,…,𝑁}𝑒^𝑛₀≤𝐶𝑒⁰₀+max_{𝑛=0,…,𝑁}𝜂^𝑛₀+𝜓(Δ𝑡)max_{𝑛=0,…,𝑁−1}𝑅^𝑛₀. (30)

Proof. In view of (23), we have the following relation using triangular inequality: max_{𝑛=0,…,𝑁}𝑒^𝑛₀≤max_{𝑛=0,…,𝑁}𝜙^𝑛₀+max_{𝑛=0,…,𝑁}𝜂^𝑛₀. (31) We bound 𝜙^𝑛 in (31) by first taking partial sums of the first term of (26), seconded by adding the remaining terms from 𝑛=1 to 𝑚 for 1≤𝑚≤𝑁−1, and multiplying both sides by 𝜓(Δ𝑡) yields 𝜙^𝑚+1−𝜙^𝑚𝜓(Δ𝑡),𝑣_ℎ−𝜙¹−𝜙⁰𝜓(Δ𝑡),𝑣_ℎ+𝜓(Δ𝑡)^𝑚_𝑛=1𝑎𝜙^𝑚,𝑣_ℎ=𝜓(Δ𝑡)^𝑚_𝑛=1𝑟^𝑛,𝑣_ℎ. (32) In view of (28), we can define Φ^𝑚=𝜓(Δ𝑡)^𝑚_𝑛=1𝜙^𝑚,Φ⁰=0 (33) and using this in (32) we have 𝜙^𝑚+1−𝜙^𝑚𝜓(Δ𝑡),𝑣_ℎ+𝑎Φ^𝑚,𝑣_ℎ=𝑅^𝑚,𝑣_ℎ,∀𝑣_ℎ∈𝑉_ℎ. (34) If we choose 𝑣_ℎ=𝜙^𝑚+1+𝜙^𝑚∈𝑉_ℎ in (34) and multiply the result by 𝜓(Δ𝑡), this yields 𝜙^𝑚+12₀−𝜙^𝑚2₀+𝜓(Δ𝑡)𝑎Φ^𝑚,𝜙^𝑚+1+𝜙^𝑚=𝜓(Δ𝑡)^𝑛−1_𝑚=1𝑅^𝑚,𝜙^𝑚+1+𝜙^𝑚,for0≤𝑚≤𝑁−1. (35) Summing this from 𝑚=0 to 𝑚=𝑛−1 for 1≤𝑛≤𝑁−1 gives 𝜙^𝑛2₀−𝜙⁰²₀+𝜓(Δ𝑡)^𝑛−1_𝑚=0𝑎Φ^𝑚,𝜙^𝑚+1+𝜙^𝑚=𝜓(Δ𝑡)^𝑛−1_𝑚=0𝑅^𝑚,𝜙^𝑚+1+𝜙^𝑚. (36) Since 𝑎 is symmetric, Φ⁰=0, and Φ^𝑚+1−Φ^𝑚−1=𝜓(Δ𝑡)𝜙^𝑚+𝜙^𝑚+1,𝑚=1,2,…,𝑁−1, (37) then, we deduce in view of the third term of (36) 𝜓(Δ𝑡)^𝑛−1_𝑚=0𝑎Φ^𝑚,𝜙^𝑚+1+𝜙^𝑚=^𝑛−1_𝑚=0𝑎Φ^𝑚,Φ^𝑚+1+Φ^𝑚−1=^𝑛−1_𝑚=0𝑎Φ^𝑚,Φ^𝑚+1−^𝑛−1_𝑚=0𝑎Φ^𝑚+1,Φ^𝑚−1=𝑎Φ^𝑛−1,Φ^𝑛. (38) By symmetry and coercivity properties of 𝑎 and the fact that Φ^𝑛−Φ^𝑛−1=𝜓(Δ𝑡)𝜙^𝑛for𝑛=1,2,…,𝑁, (39) we have in view of (38) 𝜓(Δ𝑡)²2𝑎𝜙^𝑛,𝜙^𝑛=12𝑎Φ^𝑛−Φ^𝑛−𝑖,Φ^𝑛−Φ^𝑛−1=12𝑎Φ^𝑛,Φ^𝑛−𝑎Φ^𝑛,Φ^𝑛−1+12𝑎Φ^𝑛−1,Φ^𝑛−1 (40) and so 𝜓(Δ𝑡)²2𝑎𝜙^𝑛,𝜙^𝑛≥−𝑎Φ^𝑛,Φ^𝑛−1 (41) and this together with (36) yields 𝜙^𝑛2₀−𝜓(Δ𝑡)²2𝑎𝜙^𝑛,𝜙^𝑛≤𝜙⁰²₀+𝜓(Δ𝑡)^𝑛−1_𝑚=0𝑅^𝑛,𝜙^𝑚+𝜙^𝑚+1,for1≤𝑛≤𝑁. (42) By Poincare inequality together with the continuity and coercivity of 𝑎, we have the following inequality: 𝜙^𝑛2₀−𝜓(Δ𝑡)²2𝜙^𝑛2₀≤𝜙⁰²₀+𝜓(Δ𝑡)^𝑛−1_𝑚=0𝑅^𝑛,𝜙^𝑚+𝜙^𝑚+1. (43) We suppose at this stage that ℎ and Δ𝑡 satisfy the CFL condition such that 𝜓(Δ𝑡)/ℎ<1. With this condition, the previous inequality becomes 𝐶_Δ𝑡𝜙^𝑛2₀≤𝜙⁰²₀+𝜓(Δ𝑡)^𝑛−1_𝑚=0𝑅^𝑛,𝜙^𝑚+𝜙^𝑚+1,1≤𝑛≤𝑁 (44) from where we have 1/2<𝐶_Δ𝑡=1−(𝜓(Δ𝑡))²/2ℎ²<1 and 1/2<𝐶_Δ𝑡<1. By the use of Cauchy-Schwarz inequality on (44) and some algebraic manipulations, we have 𝐶_Δ𝑡𝜙^𝑛2₀≤𝜙⁰²₀+𝜓(Δ𝑡)^𝑛−1_𝑚=0𝑅^𝑚,𝜙^𝑚+𝜙^𝑚+1≤𝜙⁰²₀+𝜓(Δ𝑡)^𝑛−1_𝑚=0𝑅^𝑚₀𝜙^𝑚₀+𝜙^𝑚+1₀≤𝜙⁰²₀+2𝜓(Δ𝑡)^𝑛−1_𝑚=0𝑅^𝑚₀max_{𝑚=0,…,𝑛}𝜙^𝑚₀≤𝜙⁰²₀+𝐶_Δ𝑡2max_{𝑚=1,…,𝑛}𝜙^𝑚2₀+2𝐶_Δ𝑡𝜓(Δ𝑡)^𝑛−1_𝑚=0𝑅^𝑛₀²for0≤𝑚≤𝑛 (45) and since the right-hand side is independent of 𝑛, then 𝐶_Δ𝑡−𝐶_Δ𝑡2𝜙^𝑚2₀≤𝜙⁰²₀+2𝐶_Δ𝑡𝜓(Δ𝑡)^𝑁−1_𝑛=0𝑅^𝑛₀² (46) which then implies that 𝐶_Δ𝑡2𝜙^𝑚2₀≤𝜙⁰²₀+2𝐶_Δ𝑡𝜓(Δ𝑡)^𝑁−1_𝑛=0𝑅^𝑛₀². (47) Furthermore, we have 𝜙^𝑛2₀≤2𝐶_Δ𝑡𝜙^𝑛2₀+4𝜓(Δ𝑡)²𝐶_Δ𝑡^𝑁−1_𝑛=0𝑅^𝑛₀² (48) and hence the following result: max_{𝑛=0,…,𝑁}𝜙^𝑛₀≤2𝐶_Δ𝑡𝜙⁰₀+2𝜓(Δ𝑡)𝐶_Δ𝑡^𝑁−1_𝑛=1𝑅^𝑛₀. (49) With (49) and the fact that 𝜙⁰₀≤𝑒⁰₀+𝜂⁰₀ (50) we have in view of (31) the desired estimate of the proposition.

We now proceed to bound the term containing 𝑅^𝑛 on the right-hand side of Proposition 3. We achieve this by estimating in the 𝐿²-norm the function 𝑟^𝑛 for the cases 𝑛=0 and follow by the case when 𝑛≥1.

Lemma 4. There holds 𝑟⁰₀≤𝐶𝜓(Δ𝑡)⁻¹ℎ²𝜕𝑢𝜕𝑡_𝐿²_{[(0,𝑇);𝐻}²_(Ω)](𝜕+𝜓Δ𝑡)²𝑢𝜕𝑡²_{𝐿2[(0,𝑇);𝐿2(Ω)]}, (51) with a constant 𝐶>0 that is independent of ℎ and the mesh size.

Proof. In view of (27), we have 𝑟⁰=(𝜓(Δ𝑡))⁻²(𝜙¹−𝜙⁰). We estimate ‖𝜙¹−𝜙⁰‖₀ by taking 𝑣_ℎ∈𝑉_ℎ arbitrary as follows: 𝜙¹−𝜙⁰,𝑣_ℎ=𝑤¹_ℎ−𝑈¹_ℎ,𝑣_ℎ−𝑤⁰_ℎ−𝑈⁰_ℎ,𝑣_ℎ=𝜋_ℎ−𝐼𝑢₁−𝑢₀,𝑣_ℎ+𝑢₁−𝑈¹,𝑣_ℎ, (52) where we have used (𝑢₀−𝑈⁰_ℎ,𝑣_ℎ)=(𝑢₀−𝑃_ℎ𝑢₀,𝑣)=0 in view of (19). It now follows from (52) that 𝜋_ℎ−𝐼𝑢₁−𝑢₀,𝑣_ℎ≤^𝑡1₀𝜕𝜋𝜕𝑡^ℎ−𝐼𝑢(⋅,𝑠),𝑣_ℎ≤𝑑𝑠ByLemma1_𝑡10𝜋_ℎ()−𝐼𝜕𝑢⋅,𝑠𝜕𝑡,𝑣_ℎ()()ℎ𝑑𝑠By18≤𝐶𝜓Δ𝑡²𝜕𝑢𝜕𝑡_𝐿²_{[()0,𝑇;𝐻}²_(Ω)]𝑣_ℎ0. (53) We now proceed to estimate the term (𝑢₁−𝑈¹_ℎ,𝑣) in (52) as follows: by Taylor’s formula and the fact 𝑢(⋅,0)=𝑢₀ and 𝜕𝑢(⋅,0)/𝜕𝑡=𝑢₁ in (1)–(3), we have 𝑢₁=𝑢₀+𝜓(Δ𝑡)𝑢₁+𝜓(Δ𝑡)²2𝜕²𝑢(⋅,0)𝜕𝑡²+12^𝑡1₀𝜓(Δ𝑡)−𝑠²𝜕³𝑢(⋅,𝑠)𝜕𝑡³𝑑𝑠. (54) In view of the definition of 𝑈¹_ℎ in (20) and the fact that 𝑢₀−𝑃_ℎ𝑢₀,𝑣=0,𝑢₁−𝑃_ℎ𝑢₁,𝑣=0 (55) we have 𝑈¹_ℎ=𝑈⁰_ℎ+𝜓(Δ𝑡)𝑃_ℎ𝑢₁+𝜓(Δ)²2𝑈⁰_ℎ. (56) We then deduce from (54) and (56) that 𝑢₁−𝑈¹_ℎ,𝑣_ℎ=𝜓(Δ𝑡)²2𝜕²𝑢(⋅,0)𝜕𝑡²−𝑈⁰_ℎ,𝑣_ℎ+12^𝑡1₀𝜓(Δ𝑡)−𝑠²𝜕³𝑢(⋅,𝑠)𝜕𝑡³,𝑣_ℎ𝑑𝑠. (57) But by the definition of 𝑈⁰_ℎ in (20) we have in view of (21) that 𝜕²𝑢(⋅,𝑠)𝜕𝑡²−𝑈⁰_ℎ,𝑣_ℎ=𝑓⁰,𝑣_ℎ−𝑎𝑢₀,𝑣_ℎ−𝑓⁰,𝑣_ℎ+𝑎𝑢₀,𝑣_ℎ=0 (58) implying that 𝑢₁−𝑈¹_ℎ≤1,𝑣2^𝑡1₀𝜓(Δ𝑡)−𝑠²𝜕³𝑢(⋅,𝑠)𝜕𝑡³𝜓(),𝑣𝑑𝑠≤𝐶Δ𝑡³𝜕³𝑢𝜕𝑡³_{𝐿2[(0,𝑇);𝐿2(Ω)]}‖𝑣‖₀. (59) Since 𝜙¹−𝜙⁰∈𝑉_ℎ, then (52) together with (53) and (59) yields 𝜙¹−𝜙⁰₀≤𝐶𝜓(Δ𝑡)ℎ²𝜕𝑢𝜕𝑡_𝐿²_{[()0,𝑇;𝐻}²_(Ω)]+𝜓(Δ𝑡)³𝜕³𝑢𝜕𝑡³_{𝐿2[(0,𝑇);𝐿2(Ω)]} (60) and dividing throughout by (𝜓(Δ𝑡))² yields 𝑟⁰₀≤𝐶𝜓(Δ𝑡)⁻¹ℎ²𝜕𝑢𝜕𝑡_𝐿²_{[()0,𝑇;𝐻}²_(Ω)](𝜕+𝜓Δ𝑡)³𝑢𝜕𝑡³_{𝐿2[()0,𝑇;𝐿2(Ω)]} (61) as required.

Lemma 5. For 1≤𝑛≤𝑁−1, there holds 𝑟^𝑛₀ℎ≤𝐶²𝜓(Δ𝑡)^𝑡𝑛+1_𝑡𝑛−1𝜕²𝑢(⋅,𝑠)𝜕𝑡²₂()𝑑𝑠+𝜓Δ𝑡^𝑡𝑛+1_𝑡𝑛−1𝜕⁴𝑢(⋅,𝑠)𝜕𝑡⁴₀𝑑𝑠 (62) with a constant 𝐶>0 that is independent of ℎ and mesh size.

Proof. In view of (27), we have 𝑟^𝑛=𝛿²𝑤^𝑛_ℎ−𝜕²𝑢^𝑛𝜕𝑡²,for𝑛=1,2,…,𝑁−1. (63) By the triangular inequality, we have 𝑟^𝑛₀=𝛿²𝑤^𝑛_ℎ−𝜕²𝑢^𝑛𝜕𝑡²₀≤𝛿²(𝜋_ℎ−𝐼)𝑢^𝑛₀+𝛿²𝑤^𝑛_ℎ−𝜕²𝑢^𝑛𝜕𝑡²₀ (64) and using the following identity 𝑣⋅,𝑡_𝑛+1−2𝑣⋅,𝑡_𝑛+𝑣⋅,𝑡_𝑛−1=Δ𝑡^𝑡𝑛+1_𝑡𝑛−11−𝑠−𝑡_𝑛𝜕Δ𝑡²𝑣(⋅,𝑠)𝜕𝑡²𝑑𝑠 (65) on the first term of the right-hand side of (64) we have 𝛿²𝜋_ℎ−𝐼𝑢^𝑛₀≤1𝜓(Δ𝑡)^𝑡𝑛+1_𝑡𝑛−11−𝑠−𝑡_𝑛𝜕Δ𝑡²(𝜋_ℎ−𝐼)𝑢𝜕𝑡²₀ℎ𝑑𝑠≤𝐶²𝜓(Δ𝑡)^𝑡𝑛+1_𝑡𝑛−1𝜕²𝑢(⋅,𝑠)𝜕𝑡²₂𝑑𝑠 (66) after the use of Lemma 1 and (14).
The second term of (64) can be estimated by the use of the identity 𝛿²𝑤^𝑛_ℎ−𝜕²𝑢^𝑛𝜕𝑡²=16(Δ𝑡)^{2𝑡𝑛+1𝑡}_𝑛−1Δ𝑡−𝑠−𝑡_𝑛³𝜕⁴𝑢(⋅,𝑠)𝜕𝑡⁴𝑑𝑠 (67) which is obtained from Taylor’s formulae with integral remainder. This is deduced as follows: 𝛿²𝑤^𝑛_ℎ−𝜕²𝑢^𝑛𝜕𝑡²₀≤16(Δ𝑡)^{2𝑡𝑛+1𝑡}_𝑛−1Δ𝑡−𝑠−𝑡_𝑛³𝜕⁴𝑢(⋅,𝑠)𝜕𝑡⁴₀≤1𝑑𝑠6(Δ𝑡)^{2𝑡𝑛+1𝑡}_𝑛−1𝜓(Δ𝑡)³𝜕⁴𝑢(⋅,𝑠)𝜕𝑡⁴₀≤𝜓()𝑑𝑠Δ𝑡6^𝑡𝑛+1_𝑡𝑛−1𝜕⁴𝑢(⋅,𝑠)𝜕𝑡⁴₀𝑑𝑠 (68) using the relation that 𝜓(Δ𝑡)−|𝑠−𝑡_𝑛|≤𝜓(Δ𝑡) in (68).

In view of (64) using (66) and (68), we have the desired result for the bound of ‖𝑟^𝑛‖₀. We now assemble Lemmas 4 and 5 in the next proposition to obtain the bound 𝑅^𝑛 as follows.

Proposition 6. For 0≤𝑛≤𝑁−1, there holds 𝑅^𝑛₀≤𝐶𝜓(Δ𝑡)²𝜕³𝑢𝜕𝑡³_{𝐿2)(0,𝑇;𝐿2(Ω)}+𝜕⁴𝑢𝜕𝑡⁴_{𝐿2[()0,𝑇;𝐿2(Ω)]}+𝐶ℎ²𝜕𝑢𝜕𝑡_𝐿²_{(0,𝑇);𝐻}²_(Ω)+𝜕²𝑢𝜕𝑡²_{𝐿2[(0,𝑇);𝐻2(Ω)]}. (69)

Proof. Using the bounds on ‖𝑟^𝑛‖₀ in Lemmas 4 and 5, we have 𝑅^𝑛₀≤𝜓(Δ𝑡)𝑟⁰₀+𝜓(Δ𝑡)^𝑁−1_𝑚=1𝑟^𝑚≤𝐶𝜓(Δ𝑡)²𝜕³𝑢𝜕𝑡³_{𝐿2[()0,𝑇;𝐿2(Ω)]}+𝜕⁴𝑢𝜕𝑡⁴_{𝐿2[()0,𝑇;𝐿2(Ω)]}+𝐶ℎ²𝜕𝑢𝜕𝑡_𝐿²_{[(0,𝑇);𝐻}²_(Ω)]+𝜕²𝑢𝜕𝑡²_{𝐿2[(0,𝑇);𝐻2(Ω)]} (70) and the proof is completed.

With all these results, we are now in the position to prove the main result as follows.

Proof of Theorem 2. Using Proposition 3 and the fact that 𝜓(Δ𝑡)^𝑁−1_𝑛=0𝑅^𝑛₀≤𝑇max_{𝑛=0,…,𝑁}𝑅^𝑛₀ (71) we have max_{𝑛=0,…,𝑁}𝑒^𝑛≤𝐶𝑒⁰₀+max_{𝑛=0,…,𝑁}𝜂^𝑛₀+𝑇max_{𝑛=0,…,𝑁−1}𝑅^𝑛₀. (72) By the use of Lemma 1, we can bound the second term on the right-hand side as follows. max_{𝑛=0,…,𝑁}𝜂^𝑛₀≤𝐶ℎ²‖𝑢‖_{𝐿2[(0,𝑇);𝐻2(Ω)]}. (73) From the approximation property of the 𝐿²-projection, we have 𝑒⁰₀≤𝐶ℎ²𝑢₀₂≤𝐶ℎ²‖𝑢‖_{𝐿2[(0,𝑇);𝐻2(Ω)]}. (74) Finally, by the bound of max_{𝑛=0,…,𝑁−1}‖𝑅^𝑛‖₀ obtained via Proposition 6, we have the result max_{𝑛=0,…,𝑁}𝑒^𝑛₀≤𝐶ℎ²‖𝑢‖_{𝐿2[(0,𝑇);𝐻2(Ω)]}+𝐶𝜓(Δ𝑡)²𝜕³𝑢𝜕𝑡³_{𝐿2[()0,𝑇;𝐿2(Ω)]}+𝜕⁴𝑢𝜕𝑡⁴_{𝐿2[()0,𝑇;𝐿2(Ω)]}. (75) In view of the relationship 𝜓(Δ𝑡)=2sinΔ𝑡2≈Δ𝑡+0(Δ𝑡)² (76) as proposed for such schemes in Mickens [12], we have the required estimate as Δ𝑡→0. Furthermore, using the fact that its uniform convergence results imply pointwise convergence for 𝑥∈Ω^′ completes the proof.

5. Numerical Experiments

In this section, we present the numerical experiments on problem (1) using both the standard finite difference (SFD) and NSFD-CG methods. These experiments are performed in Ω=(0,1)²×(0,𝑇) where Ω was discretized using regular meshes of sizes ℎ=1/𝑀 in the space and Δ𝑡=𝑇/𝑁 in the time. The 𝑀 and 𝑁 in such a discretization denote the number of nodes and time respectively. The initial data was considered to be 𝑢₀(𝑥)=𝑥₁𝑥₂(1−𝑥₁)(1−𝑥₂) and 𝑢₁(𝑥)=0 where these data were deduced from the exact solution 𝑢(𝑥,𝑡)=𝑥₁𝑥₂1−𝑥₁1−𝑥₂cos(𝜋𝑡). (77) Using the previous exact solution we obtained the right-hand side 𝑓 of (1). In the computation, (1) together with (17)–(20) led to a system of equations 𝐀𝐗=𝐛. (78) In solving for 𝐗 in the above system, we took the following values of 𝑀=10, 15, 20, 25, 50, and 100. For a fix 𝑀=20, 𝑁=10, and 𝑇=1, we had Figures 1–3 illustrating various solutions corresponding to their respective schemes.

Figure 1 

The Exact solution.

Figure 2 

Approximate solution for SFD-CG scheme.

Figure 3 

Approximate solution for NSFD-CG scheme.

Figure 1 shows the exact solution, Figure 2 the solution from the SFD-CG, and Figure 3 the solution from the NSFD-CG schemes, followed by Tables 1 and 2 which demonstrate various optimal rates of convergence in both 𝐻¹ and 𝐿²-norms of these schemes.

Table 1 

NSFD method error in both 𝐿² and 𝐻¹-norms.

Table 2 

SFD method error in both 𝐿² and 𝐻¹-norms.

These optimal rates of convergence were calculated by using the formula 𝑒Rate=ln₂/𝑒₁lnℎ₂/ℎ₁, (79) where ℎ₁ and ℎ₂ together with 𝑒₁ and 𝑒₂ are successive triangle diameters and errors, respectively. These results are self-explanatory and we could conclude that the results as shown by these experiments exhibit the desired results as expected from our theoretical analysis.

6. Conclusion

We presented a reliable scheme of the wave equation consisting of the nonstandard finite difference method in time and the continuous Galerkin method in the space variable (NSFD-CG). We proved theoretically that the numerical solution obtained from this scheme attains the optimal rate of convergence in both the energy and the 𝐿²-norms. Furthermore, we showed that the scheme under investigation replicates the properties of the exact solution of the wave equation. We proceeded by the help of a numerical example and showed that the optimal rate of convergence as proved theoretically is guaranteed in both the energy and the 𝐿²-norms. This convergence results hold for any fully discrete NSFD-CG method where the scheme under consideration has a bilinear form which is symmetric, continuous, and coercive.

The method presented in this paper can be extended to the nonlinear hyperbolic or parabolic problems with either smooth or nonsmooth domain if at all these cases followed the procedure as proposed by Mickens [12]. We will also exploit another form of nonstandard finite differential method as proposed in [18, 19].

Conflict of Interests

The author declares that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

Big thanks goes to the Almighty God from whose strength I draw my inspiration to complete this project. The research contained in this paper has been supported by the University of Pretoria, South Africa. The authors would like to thank the reviewers who spent a lot of time to go through the paper. Thanks also go to the following colleagues Dr. Chapwanya, Dr. Appadu, Mr. Mbehou, and Mr. Aderogba who help checked the revised version of the work.