Abstract

we propose a new high-order approximation for the solution of two-space-dimensional quasilinear hyperbolic partial differential equation of the form 𝑢𝑡𝑡=𝐴(𝑥,𝑦,𝑡,𝑢)𝑢𝑥𝑥+𝐵(𝑥,𝑦,𝑡,𝑢)𝑢𝑦𝑦+𝑔(𝑥,𝑦,𝑡,𝑢,𝑢𝑥,𝑢𝑦,𝑢𝑡), 0<𝑥, 𝑦<1, 𝑡>0 subject to appropriate initial and Dirichlet boundary conditions , where 𝑘>0 and >0 are mesh sizes in time and space directions, respectively. We use only five evaluations of the function 𝑔 as compared to seven evaluations of the same function discussed by (Mohanty et al., 1996 and 2001). We describe the derivation procedure in details and also discuss how our formulation is able to handle the wave equation in polar coordinates. The proposed method when applied to a linear hyperbolic equation is also shown to be unconditionally stable. Some examples and their numerical results are provided to justify the usefulness of the proposed method.

1. Introduction

We consider the following two-space dimensional quasilinear hyperbolic partial differential equation:𝑢𝑡𝑡=𝐴(𝑥,𝑦,𝑡,𝑢)𝑢𝑥𝑥+𝐵(𝑥,𝑦,𝑡,𝑢)𝑢𝑦𝑦+𝑔𝑥,𝑦,𝑡,𝑢,𝑢𝑥,𝑢𝑦,𝑢𝑡,0<𝑥,𝑦<1,𝑡>0(1.1) subject to the initial conditions 𝑢(𝑥,𝑦,0)=𝜙(𝑥,𝑦),𝑢𝑡(𝑥,𝑦,0)=𝜓(𝑥,𝑦),0𝑥,𝑦1,(1.2) and the boundary conditions𝑢(0,𝑦,𝑡)=𝑎0(𝑦,𝑡),𝑢(1,𝑦,𝑡)=𝑎1(𝑦,𝑡),0𝑦1,𝑡0,(1.3a)𝑢(𝑥,0,𝑡)=𝑏0(𝑥,𝑡),𝑢(𝑥,1,𝑡)=𝑏1(𝑥,𝑡),0𝑥1,𝑡0,(1.3b)where (1.1) is assumed to satisfy the hyperbolicity condition 𝐴(𝑥,𝑦,𝑡,𝑢)>0 and 𝐵(𝑥,𝑦,𝑡,𝑢)>0 in the solution region Ω{(𝑥,𝑦,𝑡)0<𝑥,𝑦<1,𝑡>0}. Further we assume that 𝑢(𝑥,𝑦,𝑡)𝐶6, 𝐴(𝑥,𝑦,𝑡,𝑢),𝐵(𝑥,𝑦,𝑡,𝑢)𝐶4, and 𝜙(𝑥,𝑦) and 𝜓(𝑥,𝑦) are sufficiently differentiable function of as higher-order as possible.

Second-order quasilinear hyperbolic partial differential equation with appropriate initial and boundary conditions serves as models in many branches of physics and technology. Ciment and Leventhal [1, 2] have discussed operator compact implicit method to solve wave equation. Using finite volume technique, Schwartzkopff et al. [3] have studied high-order ADER schemes for first-order linear hyperbolic systems. Even much higher order ADER finite volume schemes have been presented and analyzed by Schwartzkopff et al. [4] on Cartesian meshes, and by Dumbser and Käser [5] on general two- and three-dimensional unstructured meshes. In 1996, Mohanty et al. [6] have developed a high-order numerical method for the solution of two-space dimensional second-order nonlinear hyperbolic equation. Later, Mohanty et al. [7] have extended their technique to solve second-order quasilinear hyperbolic equations. In both cases they have used seven evaluations of the function 𝑔 and 19-grid points. It has been shown that the linear schemes discussed in [6, 7] are conditionally stable. In the recent past, many researchers (see [816]) have developed unconditionally stable implicit finite difference methods for the solution of two-space dimensional linear hyperbolic equations with significant first derivative terms. Most recently, Mohanty and Singh [17] have derived a high accuracy numerical method based on Numerov type discretization for the solution of one space dimensional nonlinear hyperbolic equations, in which they have shown that the linear scheme is unconditionally stable.

In this paper, using nineteen grid-points, we derive a new compact three-level implicit numerical method of accuracy two in time and four in space for the solution of two-space dimensional quasilinear hyperbolic equation (1.1). In this method we require only five evaluations of the function 𝑔 as compared to seven evaluations of the same function discussed in [6, 7]. In the next section, we give formulation of the method. In Section 3, we give the complete derivation of the method. In Section 4, we discuss the application of the proposed method to two-space dimensional wave equation in polar coordinates and discuss the stability analysis. In this section, we modify our method in such a way that the solution retains its order and accuracy everywhere in the vicinity of the singularity. In Section 5, we discuss superstable method for two-space dimensional telegraphic equation. In Section 6, we examine our new method over a set of linear and nonlinear second-order hyperbolic equations whose exact solutions are known and compared the results with the results of other known methods. Concluding remarks are given in Section 7.

2. Formulation of the Numerical Method

In this section, we aim to discuss a numerical method for the solution of nonlinear wave equation: 𝑢𝑡𝑡=𝐴(𝑥,𝑦,𝑡)𝑢𝑥𝑥+𝐵(𝑥,𝑦,𝑡)𝑢𝑦𝑦+𝑔𝑥,𝑦,𝑡,𝑢,𝑢𝑥,𝑢𝑦,𝑢𝑡,0<𝑥,𝑦<1,𝑡>0.(2.1)

Let >0 and 𝑘>0 be the mesh spacing in the space and time directions, respectively. We replace the solution region Ω={(𝑥,𝑦,𝑡)0<𝑥,𝑦<1,𝑡>0} by a set of grid points (𝑥𝑙,𝑦𝑚,𝑡𝑗), where 𝑥𝑙=𝑙, 𝑙=0(1)𝑁+1, 𝑦𝑚=𝑚, 𝑚=0(1)𝑁+1, and 𝑡𝑗=𝑗𝑘, 0<𝑗<𝐽, 𝑁 and 𝐽 being positive integers and (𝑁+1)=1. Let 𝑝=𝑘/>0 be the mesh ratio parameter. The exact values of 𝑢(𝑥,𝑦,𝑡), 𝐴(𝑥,𝑦,𝑡), and 𝐵(𝑥,𝑦,𝑡) at the grid point (𝑥𝑙,𝑦𝑚,𝑡𝑗) are denoted by 𝑈𝑗𝑙,𝑚, 𝐴𝑗𝑙,𝑚, and 𝐵𝑗𝑙,𝑚, respectively. Similarly at the grid point (𝑥𝑙,𝑦𝑚,𝑡𝑗), we denote 𝐴𝑗𝑥𝑙,𝑚=𝜕𝐴𝑗𝑙,𝑚/𝜕𝑥,𝐴𝑗𝑦𝑙,𝑚=𝜕𝐴𝑗𝑙,𝑚/𝜕𝑦,𝐴𝑗𝑥𝑥𝑙,𝑚=𝜕2𝐴𝑗𝑙,𝑚/𝜕𝑥2, and so forth. Let 𝑢𝑗𝑙,𝑚 be the approximate value of 𝑢(𝑥,𝑦,𝑡) at the same grid point.

At the grid point (𝑥𝑙,𝑦𝑚,𝑡𝑗), we consider the following approximations:𝑈𝑗𝑡𝑙,𝑚=𝑈𝑗+1𝑙,𝑚𝑈𝑗1𝑙,𝑚,(2𝑘)(2.2a)𝑈𝑗𝑡𝑙±1,𝑚=𝑈𝑗+1𝑙±1,𝑚𝑈𝑗1𝑙±1,𝑚(2𝑘),(2.2b)𝑈𝑗𝑡𝑙,𝑚±1=𝑈𝑗+1𝑙,𝑚±1𝑈𝑗1𝑙,𝑚±1(2𝑘),(2.2c)𝑈𝑗𝑡𝑡𝑙,𝑚=𝑈𝑗+1𝑙,𝑚2𝑈𝑗𝑙,𝑚+𝑈𝑗1𝑙,𝑚𝑘2,(2.2d)𝑈𝑗𝑡𝑡𝑙±1,𝑚=𝑈𝑗+1𝑙±1,𝑚2𝑈𝑗𝑙±1,𝑚+𝑈𝑗1𝑙±1,𝑚𝑘2,(2.2e)𝑈𝑗𝑡𝑡𝑙,𝑚±1=𝑈𝑗+1𝑙,𝑚±12𝑈𝑗𝑙,𝑚±1+𝑈𝑗1𝑙,𝑚±1𝑘2,(2.2f)𝑈𝑗𝑥𝑙,𝑚=𝑈𝑗𝑙+1,𝑚𝑈𝑗𝑙1,𝑚,(2)(2.3a)𝑈𝑗𝑥𝑙±1,𝑚=±3𝑈𝑗𝑙±1,𝑚4𝑈𝑗𝑙,𝑚±𝑈𝑗𝑙1,𝑚(2),(2.3b)𝑈𝑗𝑥𝑙,𝑚±1=𝑈𝑗𝑙+1,𝑚±1𝑈𝑗𝑙1,𝑚±1(2),(2.3c)𝑈𝑗𝑥𝑥𝑙,𝑚=𝑈𝑗𝑙+1,𝑚2𝑈𝑗𝑙,𝑚+𝑈𝑗𝑙1,𝑚2,(2.3d)𝑈𝑗𝑥𝑥𝑙,𝑚±1=𝑈𝑗𝑙+1,𝑚±12𝑈𝑗𝑙,𝑚±1+𝑈𝑗𝑙1,𝑚±12,(2.3e)𝑈𝑗𝑦𝑙,𝑚=𝑈𝑗𝑙,𝑚+1𝑈𝑗𝑙,𝑚1,(2)(2.4a)𝑈𝑗𝑦𝑙±1,𝑚=𝑈𝑗𝑙±1,𝑚+1𝑈𝑗𝑙±1,𝑚1(2),(2.4b)𝑈𝑗𝑦𝑙,𝑚±1=±3𝑈𝑗𝑙,𝑚±14𝑈𝑗𝑙,𝑚±𝑈𝑗𝑙,𝑚1(2),(2.4c)𝑈𝑗𝑦𝑦𝑙,𝑚=𝑈𝑗𝑙,𝑚+12𝑈𝑗𝑙,𝑚+𝑈𝑗𝑙,𝑚12,(2.4d)𝑈𝑗𝑦𝑦𝑙±1,𝑚=𝑈𝑗𝑙±1,𝑚+12𝑈𝑗𝑙±1,𝑚+𝑈𝑗𝑙±1,𝑚12,(2.4e)𝐺𝑗𝑙±1,𝑚𝑥=𝑔𝑙±1,𝑦𝑚,𝑡𝑗,𝑈𝑗𝑙±1,𝑚,𝑈𝑗𝑥𝑙±1,𝑚,𝑈𝑗𝑦𝑙±1,𝑚,𝑈𝑗𝑡𝑙±1,𝑚,(2.5a)𝐺𝑗𝑙,𝑚±1𝑥=𝑔𝑙,𝑦𝑚±1,𝑡𝑗,𝑈𝑗𝑙,𝑚±1,𝑈𝑗𝑥𝑙,𝑚±1,𝑈𝑗𝑦𝑙,𝑚±1,𝑈𝑗𝑡𝑙,𝑚±1,(2.5b)𝑈𝑗𝑥𝑙,𝑚=𝑈𝑗𝑥𝑙,𝑚+16𝐴𝑗𝑙,𝑚𝐺𝑗𝑙+1,𝑚𝐺𝑗𝑙1,𝑚16𝐴𝑗𝑙,𝑚𝑈𝑗𝑡𝑡𝑙+1,𝑚𝑈𝑗𝑡𝑡𝑙1,𝑚+𝐵𝑗𝑖,𝑚16𝐴𝑗𝑙,𝑚𝑈𝑗𝑦𝑦𝑙+1,𝑚𝑈𝑗𝑦𝑦𝑙1,𝑚+28𝐴𝑗𝑙,𝑚𝐴𝑗𝑥𝑙,𝑚𝑈𝑗𝑥𝑥𝑙,𝑚+𝐵𝑗𝑥𝑙,𝑚𝑈𝑗𝑦𝑦𝑙,𝑚,(2.6a)𝑈𝑗𝑦𝑙,𝑚=𝑈𝑗𝑦𝑙,𝑚+16𝐵𝑗𝑙,𝑚𝐺𝑗𝑙,𝑚+1𝐺𝑗𝑙,𝑚116𝐵𝑗𝑙,𝑚𝑈𝑗𝑡𝑡𝑙,𝑚+1𝑈𝑗𝑡𝑡𝑙,𝑚1+𝐴𝑗𝑙,𝑚16𝐵𝑗𝑙,𝑚𝑈𝑗𝑥𝑥𝑙,𝑚+1𝑈𝑗𝑥𝑥𝑙,𝑚1+28𝐵𝑗𝑙,𝑚𝐴𝑗𝑦𝑙,𝑚𝑈𝑗𝑥𝑥𝑙,𝑚+𝐵𝑗𝑦𝑙,𝑚𝑈𝑗𝑦𝑦𝑙,𝑚,(2.6b)𝐺𝑗𝑙,𝑚𝑥=𝑔𝑙,𝑦𝑚,𝑡𝑗,𝑈𝑗𝑙,𝑚,𝑈𝑗𝑥𝑙,𝑚,𝑈𝑗𝑦𝑙,𝑚,𝑈𝑗𝑡𝑙,𝑚.(2.7)

Then at each grid point (𝑥𝑙,𝑦𝑚,𝑡𝑗), 𝑙,𝑚=1,2,,𝑁, 𝑗=2,,𝐽, an approximation with accuracy of 𝑂(𝑘2+𝑘22+4) for the solution of differential equation (2.1) may be written as𝐿𝑢𝑝2𝐴𝑗𝑙,𝑚26𝐴𝑗𝑥𝑙,𝑚𝐴𝑗𝑙,𝑚𝐴𝑗𝑥𝑙,𝑚26𝐵𝑗𝑦𝑙,𝑚𝐵𝑗𝑙,𝑚𝐴𝑗𝑦𝑙,𝑚+2𝐴12𝑗𝑥𝑥𝑙,𝑚+2𝐴12𝑗𝑦𝑦𝑙,𝑚𝛿2𝑥𝑈𝑗𝑙,𝑚+𝑝2𝐵𝑗𝑙,𝑚26𝐴𝑗𝑥𝑙,𝑚𝐴𝑗𝑙,𝑚𝐵𝑗𝑥𝑙,𝑚26𝐵𝑗𝑦𝑙,𝑚𝐵𝑗𝑙,𝑚𝐵𝑗𝑦𝑙,𝑚+2𝐵12𝑗𝑥𝑥𝑙,𝑚+2𝐵12𝑗𝑦𝑦𝑙,𝑚𝛿2𝑦𝑈𝑗𝑙,𝑚+𝑝2𝐴12𝑗𝑦𝑙,𝑚𝐵𝑗𝑦𝑙,𝑚𝐵𝑗𝑙,𝑚𝐴𝑗𝑙,𝑚𝛿2𝑥2𝜇𝑦𝛿𝑦𝑈𝑗𝑙,𝑚+𝑝2𝐵12𝑗𝑥𝑙,𝑚𝐴𝑗𝑥𝑙,𝑚𝐴𝑗𝑙,𝑚𝐵𝑗𝑙,𝑚𝛿2𝑦2𝜇𝑥𝛿𝑥𝑈𝑗𝑙,𝑚+𝑝2𝐴12𝑗𝑙,𝑚+𝐵𝑗𝑙,𝑚𝛿2𝑥𝛿2𝑦𝑈𝑗𝑙,𝑚=𝑘2121𝐴𝑗𝑥𝑙,𝑚𝐴𝑗𝑙,𝑚𝑈𝑗𝑡𝑡𝑙+1,𝑚+1+𝐴𝑗𝑥𝑙,𝑚𝐴𝑗𝑙,𝑚𝑈𝑗𝑡𝑡𝑙1,𝑚+1𝐵𝑗𝑦𝑙,𝑚𝐵𝑗𝑙,𝑚𝑈𝑗𝑡𝑡𝑙,𝑚+1+1+𝐵𝑗𝑦𝑙,𝑚𝐵𝑗𝑙,𝑚𝑈𝑗𝑡𝑡𝑙,𝑚1+8𝑈𝑗𝑡𝑡𝑙,𝑚𝑘2121𝐴𝑗𝑥𝑙,𝑚𝐴𝑗𝑙,𝑚𝐺𝑗𝑙+1,𝑚+1+𝐴𝑗𝑥𝑙,𝑚𝐴𝑗𝑙,𝑚𝐺𝑗𝑙1,𝑚+1𝐵𝑗𝑦𝑙,𝑚𝐵𝑗𝑙,𝑚𝐺𝑗𝑙,𝑚+1+1+𝐵𝑗𝑦𝑙,𝑚𝐵𝑗𝑙,𝑚𝐺𝑗𝑙,𝑚1+8𝐺𝑗𝑙,𝑚+𝑇𝑗𝑙,𝑚,𝑙=1,2,,𝑁;𝑚=1,2,,𝑁;𝑗=1,2,,𝐽,(2.8) where 𝛿𝑥𝑈𝑗𝑙,𝑚=(𝑈𝑗𝑙+12,𝑚𝑈𝑗𝑙12,𝑚), 𝜇𝑥𝑈𝑗𝑙,𝑚=(1/2)(𝑈𝑗𝑙+1/2,𝑚+𝑈𝑗𝑙(1/2),𝑚), and so forth and 𝑇𝑗𝑙,𝑚=𝑂(𝑘4+𝑘42+𝑘24).

3. Derivation Procedure of the Approximation

At the grid point (𝑥𝑙,𝑦𝑚,𝑡𝑗), let us denote 𝛼𝑗𝑙,𝑚=𝜕𝑔𝜕𝑈𝑥𝑗𝑙,𝑚,𝛽𝑗𝑙,𝑚=𝜕𝑔𝜕𝑈𝑦𝑗𝑙,𝑚.(3.1)

Further at the grid point (𝑥𝑙,𝑦𝑚,𝑡𝑗), the exact solution value 𝑈𝑗𝑙,𝑚 satisfies 𝑈𝑗𝑡𝑡𝑙,𝑚𝑥𝐴𝑙,𝑦𝑚,𝑡𝑗𝑈𝑗𝑥𝑥𝑙,𝑚𝑥𝐵𝑙,𝑦𝑚,𝑡𝑗𝑈𝑗𝑦𝑦𝑙,𝑚𝑥=𝑔𝑙,𝑦𝑚,𝑡𝑗,𝑈𝑗𝑙,𝑚,𝑈𝑗𝑥𝑙,𝑚,𝑈𝑗𝑦𝑙,𝑚,𝑈𝑗𝑡𝑙,𝑚𝐺𝑗𝑙,𝑚(say).(3.2) Using Taylor expansion about the grid point (𝑥𝑙,𝑦𝑚,𝑡𝑗), first we obtain𝐿𝑢=𝑘2121𝐴𝑗𝑥𝑙,𝑚𝐴𝑗𝑙,𝑚𝑈𝑗𝑡𝑡𝑙+1,𝑚+1+𝐴𝑗𝑥𝑙,𝑚𝐴𝑗𝑙,𝑚𝑈𝑗𝑡𝑡𝑙1,𝑚+1𝐵𝑗𝑦𝑙,𝑚𝐵𝑗𝑙,𝑚𝑈𝑗𝑡𝑡𝑙,𝑚+1+1+𝐵𝑗𝑦𝑙,𝑚𝐵𝑗𝑙,𝑚𝑈𝑗𝑡𝑡𝑙,𝑚1+8𝑈𝑗𝑡𝑡𝑙,𝑚𝑘2121𝐴𝑗𝑥𝑙,𝑚𝐴𝑗𝑙,𝑚𝐺𝑗𝑙+1,𝑚+1+𝐴𝑗𝑥𝑙,𝑚𝐴𝑗𝑙,𝑚𝐺𝑗𝑙1,𝑚+1𝐵𝑗𝑦𝑙,𝑚𝐵𝑗𝑙,𝑚𝐺𝑗𝑙,𝑚+1+1+𝐵𝑗𝑦𝑙,𝑚𝐵𝑗𝑙,𝑚𝐺𝑗𝑙,𝑚1+8𝐺𝑗𝑙,𝑚𝑘+𝑂4+𝑘42+𝑘24.(3.3)

With the help of the approximations (2.2a)–(2.2e), and from (2.5a) and (2.5b), we obtain𝐺𝑗𝑙±1,𝑚=𝐺𝑗𝑙±1,𝑚23𝜕3𝑈𝑗𝑙,𝑚𝜕𝑥3𝛼𝑗𝑙,𝑚+26𝜕3𝑈𝑗𝑙,𝑚𝜕𝑦3𝛽𝑗𝑙,𝑚𝑘+𝑂2±𝑘2±3+4,𝐺𝑗𝑙,𝑚±1=𝐺𝑗𝑙,𝑚±1+26𝜕3𝑈𝑗𝑙,𝑚𝜕𝑥3𝛼𝑗𝑙,𝑚23𝜕3𝑈𝑗𝑙,𝑚𝜕𝑦3𝛽𝑗𝑙,𝑚𝑘+𝑂2±𝑘2±3+4.(3.4) Now we define the following approximations: 𝑈𝑗𝑥𝑙,𝑚=𝑈𝑗𝑥𝑙,𝑚+𝑎11𝐺𝑗𝑙+1,𝑚𝐺𝑗𝑙1,𝑚+𝑎12𝑈𝑗𝑡𝑡𝑙+1,𝑚𝑈𝑗𝑡𝑡𝑙1,𝑚+𝑎13𝑈𝑗𝑦𝑦𝑙+1,𝑚𝑈𝑗𝑦𝑦𝑙1,𝑚+2𝑎14𝑈𝑗𝑥𝑥𝑙,𝑚+𝑎15𝑈𝑗𝑦𝑦𝑙,𝑚,𝑈𝑗𝑦𝑙,𝑚=𝑈𝑗𝑦𝑙,𝑚+𝑏11𝐺𝑗𝑙,𝑚+1𝐺𝑗𝑙,𝑚1+𝑏12𝑈𝑗𝑡𝑡𝑙,𝑚+1𝑈𝑗𝑡𝑡𝑙,𝑚1+𝑏13𝑈𝑗𝑥𝑥𝑙,𝑚+1𝑈𝑗𝑥𝑥𝑙,𝑚1+2𝑏14𝑈𝑗𝑥𝑥𝑙,𝑚+𝑏15𝑈𝑗𝑦𝑦𝑙,𝑚,(3.5) where 𝑎1𝑖 and 𝑏1𝑖, 𝑖=1,2,,5, are free parameters to be determined.

By using the approximations (2.2a)–(2.4e) and (3.4), from (3.5), we get𝑈𝑗𝑥𝑙,𝑚=𝑈𝑗𝑥𝑙,𝑚+26𝑇1𝑘+𝑂2+𝑘22+4,𝑈𝑗𝑦𝑙,𝑚=𝑈𝑗𝑦𝑙,𝑚+26𝑇2𝑘+𝑂2+𝑘22+4,(3.6) where 𝑇1=112𝑎11𝐴𝑗𝑙,𝑚𝑈𝑗𝑥𝑥𝑥𝑙,𝑚𝑎+1211+𝑎12𝑈𝑗𝑡𝑡𝑥𝑙,𝑚+12𝑎11𝐵𝑗𝑙,𝑚+𝑎13𝑈𝑗𝑥𝑦𝑦𝑙,𝑚+6𝑎1412𝑎11𝐴𝑗𝑥𝑙,𝑚𝑈𝑗𝑥𝑥𝑙,𝑚+6𝑎1512𝑎11𝐵𝑗𝑥𝑙,𝑚𝑈𝑗𝑦𝑦𝑙,𝑚,𝑇2=112𝑏11𝐵𝑗𝑙,𝑚𝑈𝑗𝑦𝑦𝑦𝑙,𝑚𝑏+1211+𝑏12𝑈𝑗𝑡𝑡𝑦𝑙,𝑚+12𝑏11𝐴𝑗𝑙,𝑚+𝑏13𝑈𝑗𝑥𝑥𝑦𝑙,𝑚+6𝑏1412𝑏11𝐴𝑗𝑦𝑙,𝑚𝑈𝑗𝑥𝑥𝑙,𝑚+6𝑏1512𝑏11𝐵𝑗𝑦𝑙,𝑚𝑈𝑗𝑦𝑦𝑙,𝑚.(3.7) Thus from (2.7), we obtain 𝐺𝑗𝑙,𝑚=𝐺𝑗𝑙,𝑚+26𝑇1𝛼𝑗𝑙,𝑚+26𝑇2𝛽𝑗𝑙,𝑚𝑘+𝑂2+𝑘22+4.(3.8) Finally by the help of the approximations (3.4) and (3.8), from (2.8) and (3.3), we get 𝑇𝑗𝑙,𝑚=𝑘22𝜕363𝑈𝑗𝑙,𝑚𝜕𝑥3𝛼𝑗𝑙,𝑚+𝜕3𝑈𝑗𝑙,𝑚𝜕𝑦3𝛽𝑗𝑙,𝑚𝑇41𝛼𝑗𝑙,𝑚+𝑇2𝛽𝑗𝑙,𝑚𝑘+𝑂4+𝑘42+𝑘24.(3.9)

In order to obtain the difference approximation of 𝑂(𝑘2+𝑘22+4) the coefficient of 𝑘22 in (3.9) must be zero, which gives the values of parameters:𝑎11=116𝐴𝑗𝑙,𝑚,𝑎12=116𝐴𝑗𝑙,𝑚,𝑎13=𝐵𝑗𝑙,𝑚16𝐴𝑗𝑙,𝑚,𝑎14=𝐴𝑗𝑥𝑙,𝑚8𝐴𝑗𝑙,𝑚,𝑎15=𝐵𝑗𝑥𝑙,𝑚8𝐴𝑗𝑙,𝑚,𝑏11=116𝐵𝑗𝑙,𝑚,𝑏12=116𝐵𝑗𝑙,𝑚,𝑏13=𝐴𝑗𝑙,𝑚16𝐵𝑗𝑙,𝑚,𝑏14=𝐴𝑗𝑦𝑙,𝑚8𝐵𝑗𝑙,𝑚,𝑏15=𝐵𝑗𝑦𝑙,𝑚8𝐵𝑗𝑙,𝑚.(3.10)

Thus we obtain the difference method of 𝑂(𝑘2+𝑘22+4) for the differential equation (2.1) and the local truncation error reduces to 𝑇𝑗𝑙,𝑚=𝑂(𝑘4+𝑘42+𝑘24).

Now we consider the numerical method of 𝑂(𝑘2+𝑘22+4) for the solution of the quasilinear hyperbolic equation (1.1). Whenever the coefficients are 𝐴=𝐴(𝑥,𝑦,𝑡,𝑢) and 𝐵=𝐵(𝑥,𝑦,𝑡,𝑢), the difference scheme (2.8) needs to be modified. For this purpose, we use the following central differences: 𝐴𝑗𝑥𝑙,𝑚=𝐴𝑗𝑙+1,𝑚𝐴𝑗𝑙1,𝑚(2)+𝑂2,𝐴𝑗𝑥𝑥𝑙,𝑚=𝐴𝑗𝑙+1,𝑚2𝐴𝑗𝑙,𝑚+𝐴𝑗𝑙1,𝑚2+𝑂2,𝐴𝑗𝑦𝑙,𝑚=𝐴𝑗𝑙,𝑚+1𝐴𝑗𝑙,𝑚1(2)+𝑂2,𝐴𝑗𝑦𝑦𝑙,𝑚=𝐴𝑗𝑙,𝑚+12𝐴𝑗𝑙,𝑚+𝐴𝑗𝑙,𝑚12+𝑂2,𝐵𝑗𝑥𝑙,𝑚=𝐵𝑗𝑙+1,𝑚𝐵𝑗𝑙1,𝑚(2)+𝑂2,𝐵𝑗𝑥𝑥𝑙,𝑚=𝐵𝑗𝑙+1,𝑚2𝐵𝑗𝑙,𝑚+𝐵𝑗𝑙1,𝑚2+𝑂2,𝐵𝑗𝑦𝑙,𝑚=𝐵𝑗𝑙,𝑚+1𝐵𝑗𝑙,𝑚1(2)+𝑂2,𝐵𝑗𝑦𝑦𝑙,𝑚=𝐵𝑗𝑙,𝑚+12𝐵𝑗𝑙,𝑚+𝐵𝑗𝑙,𝑚12+𝑂2,(3.11) where𝐴𝑗𝑙,𝑚𝑥=𝐴𝑙,𝑦𝑚,𝑡𝑗,𝑈𝑗𝑙,𝑚,𝐴𝑗𝑙±1,𝑚𝑥=𝐴𝑙±1,𝑦𝑚,𝑡𝑗,𝑈𝑗𝑙±1,𝑚,𝐴𝑗𝑙,𝑚±1𝑥=𝐴𝑙,𝑦𝑚±1,𝑡𝑗,𝑈𝑗𝑙,𝑚±1,𝐵𝑗𝑙,𝑚𝑥=𝐵𝑙,𝑦𝑚,𝑡𝑗,𝑈𝑗𝑙,𝑚,𝐵𝑗𝑙±1,𝑚𝑥=𝐵𝑙±1,𝑦𝑚,𝑡𝑗,𝑈𝑗𝑙±1,𝑚,𝐵𝑗𝑙,𝑚±1𝑥=𝐵𝑙,𝑦𝑚±1,𝑡𝑗,𝑈𝑗𝑙,𝑚±1.(3.12)

By the help of the approximations (3.11), it is easy to verify that 𝐴𝑗𝑙,𝑚26𝐴𝑗𝑥𝑙,𝑚𝐴𝑗𝑙,𝑚𝐴𝑗𝑥𝑙,𝑚26𝐵𝑗𝑦𝑙,𝑚𝐵𝑗𝑙,𝑚𝐴𝑗𝑦𝑙,𝑚+2𝐴12𝑗𝑥𝑥𝑙,𝑚+2𝐴12𝑗𝑦𝑦𝑙,𝑚=𝐴𝑗𝑙,𝑚124𝐴𝑗𝑙,𝑚𝐴𝑗𝑙+1,𝑚𝐴𝑗𝑙1,𝑚2124𝐵𝑗𝑙,𝑚𝐵𝑗𝑙,𝑚+1𝐵𝑗𝑙,𝑚1𝐴𝑗𝑙,𝑚+1𝐴𝑗𝑙,𝑚1+1𝐴12𝑗𝑙+1,𝑚2𝐴𝑗𝑙,𝑚+𝐴𝑗𝑙1,𝑚+1𝐴12𝑗𝑙,𝑚+12𝐴𝑗𝑙,𝑚+𝐴𝑗𝑙,𝑚1+𝑂4,𝐵𝑗𝑙,𝑚26𝐴𝑗𝑥𝑙,𝑚𝐴𝑗𝑙,𝑚𝐵𝑗𝑥𝑙,𝑚26𝐵𝑗𝑦𝑙,𝑚𝐵𝑗𝑙,𝑚𝐵𝑗𝑦𝑙,𝑚+2𝐵12𝑗𝑥𝑥𝑙,𝑚+2𝐵12𝑗𝑦𝑦𝑙,𝑚=𝐵𝑗𝑙,𝑚124𝐵𝑗𝑙,𝑚𝐵𝑗𝑙,𝑚+1𝐵𝑗𝑙,𝑚12124𝐴𝑗𝑙,𝑚𝐴𝑗𝑙+1,𝑚𝐴𝑗𝑙1,𝑚𝐵𝑗𝑙+1,𝑚𝐵𝑗𝑙1,𝑚+1𝐵12𝑗𝑙+1,𝑚2𝐵𝑗𝑙,𝑚+𝐵𝑗𝑙1,𝑚+1𝐵12𝑗𝑙,𝑚+12𝐵𝑗𝑙,𝑚+𝐵𝑗𝑙,𝑚1+𝑂4,,andsoforth.(3.13)

Thus substituting the central difference approximations (3.11) into (2.8), we obtain the required numerical method of 𝑂(𝑘2+𝑘22+4) for the solution of quasilinear hyperbolic partial differential equation (1.1).

Note that the initial and Dirichlet boundary conditions are given by (1.2), (1.3a) and (1.3b), respectively. Incorporating the initial and boundary conditions, we can write the method (2.8) in a block tridiagonal matrix form. If the differential equation (1.1) is linear, we can solve the linear system using operator splitting method or, alternating direction implicit method (see [18, 19]); in the nonlinear case, we can use Newton-Raphson iterative method to solve the non-linear system (see [2024]).

4. Stability Analysis

In this section, we aim to discuss a stable difference scheme for the two-space dimensional linear hyperbolic equation with singular coefficients and ensure that the numerical method developed here retains its order and accuracy.

Let us consider the equation of the following form:𝑢𝑡𝑡=𝑢𝑥𝑥+𝑢𝑦𝑦+𝐷(𝑥)𝑢𝑥+𝑓(𝑥,𝑦,𝑡),0<𝑥,𝑦<1,𝑡>0.(4.1) Applying the approximation (2.8) to the differential equation (4.1), we obtain𝑝2𝛿2𝑥+𝛿2𝑦+16𝛿2𝑥𝛿2𝑦𝑈𝑗𝑙,𝑚=𝑘212𝑈𝑗𝑡𝑡𝑙+1,𝑚+𝑈𝑗𝑡𝑡𝑙1,𝑚+𝑈𝑗𝑡𝑡𝑙,𝑚+1+𝑈𝑗𝑡𝑡𝑙,𝑚1+8𝑈𝑗𝑡𝑡𝑙,𝑚𝑘212𝐺𝑗𝑙+1,𝑚+𝐺𝑗𝑙1,𝑚+𝐺𝑗𝑙,𝑚+1+𝐺𝑗𝑙,𝑚1+8𝐺𝑗𝑙,𝑚𝑘+𝑂4+𝑘42+𝑘24𝑙=1,2,,𝑁;𝑚=1,2,,𝑁;𝑗=2,3,,𝐽,(4.2) where we denote𝐷𝑙𝑥=𝐷𝑙,𝑓𝑗𝑙,𝑚𝑥=𝑓𝑙,𝑦𝑚,𝑡𝑗,𝐺𝑗𝑙±1,𝑚=𝐷𝑙±1𝑈𝑗𝑥𝑙±1,𝑚+𝑓𝑗𝑙±1,𝑚,𝐺𝑗𝑙,𝑚±1=𝐷𝑙𝑈𝑗𝑥𝑙,𝑚±1+𝑓𝑗𝑙,𝑚±1,𝐺𝑗𝑙,𝑚=𝐷𝑙𝑈𝑗𝑥𝑙,𝑚+𝑓𝑗𝑙,𝑚.(4.3)

Note that the linear scheme (4.2) is of 𝑂(𝑘2+𝑘22+4) for the solution of differential equation (4.1). However, the scheme (4.2) fails to compute at 𝑙=1, if the coefficient 𝐷(𝑥) and/or the function 𝑓(𝑥,𝑦,𝑡) contains the singular terms like 1/𝑥,1/𝑥2,, and so forth. We overcome this difficulty by modifying the method (4.2) in such a way that the solution retains its order and accuracy everywhere in the vicinity of singularity 𝑥=0.

We need the following approximations: 𝐷𝑙±1=𝐷𝑙±𝐷𝑥𝑙+22𝐷𝑙𝑥𝑥±𝑂3,𝑓𝑗𝑙±1,𝑚=𝑓𝑗𝑙,𝑚±𝑓𝑗𝑥𝑙,𝑚+22𝑓𝑗𝑥𝑥𝑙,𝑚±𝑂3,𝑓𝑗𝑙,𝑚±1=𝑓𝑗𝑙,𝑚±𝑓𝑗𝑦𝑙,𝑚+22𝑓𝑗𝑦𝑦𝑙,𝑚±𝑂3.(4.4)

Using the approximations (4.4) in the scheme (4.2) and neglecting local truncation error, we get11+𝛿122𝑥+𝑅02𝜇𝑥𝛿𝑥+𝛿2𝑦𝛿122𝑡𝑢𝑗𝑙,𝑚=𝑝2𝑅1𝛿2𝑥+𝑅22𝜇𝑥𝛿𝑥+𝛿2𝑦+𝑅06𝛿2𝑦2𝜇𝑥𝛿𝑥+16𝛿2𝑥𝛿2𝑦𝑢𝑗𝑙,𝑚+𝑓,(4.5) where𝑅0=𝐷𝑙2,𝑅1=1+2122𝐷𝑥𝑙+𝐷𝑙2,𝑅2=𝑅0+3𝐷24𝑥𝑥𝑙+𝐷𝑙𝐷𝑥𝑙,𝑘𝑓=21212𝑓𝑗𝑙,𝑚+2𝑓𝑗𝑥𝑥𝑙,𝑚+𝑓𝑗𝑦𝑦𝑙,𝑚+𝐷𝑙𝑓𝑗𝑥𝑙,𝑚.(4.6) The previous scheme in product form may be written as11+𝛿122𝑥+𝑅02𝜇𝑥𝛿𝑥𝛿1+2𝑦𝛿122𝑡𝑢𝑗𝑙,𝑚=𝑝2𝑅1𝛿2𝑥+𝑅22𝜇𝑥𝛿𝑥+𝛿2𝑦+𝑅06𝛿2𝑦2𝜇𝑥𝛿𝑥+16𝛿2𝑥𝛿2𝑦𝑢𝑗𝑙,𝑚+𝑓.(4.7)

It is difficult to find the stability region for the scheme (4.7). In order to obtain a valid stability region, we may modify the scheme (4.7) (see [23]) as11+𝑅121𝛿2𝑥+𝑅22𝜇𝑥𝛿𝑥𝛿1+2𝑦𝛿122𝑡𝑢𝑗𝑙,𝑚=𝑝2𝑅1𝛿2𝑥+𝑅22𝜇𝑥𝛿𝑥+𝛿2𝑦+16𝑅1𝛿2𝑥+𝑅22𝜇𝑥𝛿𝑥𝛿2𝑦𝑢𝑗𝑙,𝑚+𝑓.(4.8)

The additional terms added in (4.7) and (4.8) are of high orders and do not affect the accuracy of the scheme.

For stability, we put 𝑢𝑗𝑙,𝑚=𝐴𝑙𝐵𝑚𝜉𝑗𝑒𝑖𝛽𝑙𝑒𝑖𝛾𝑚 (where 𝜉=𝑒𝑖𝜙 such that |𝜉|=1) in the homogeneous part of the scheme (4.8), and we get𝜉2+𝜉1=4sin2𝜙2=𝑝2𝑀1+𝑀2+(1/6)𝑀1𝑀21+(1/12)𝑀11+(1/12)𝑀2,(4.9) where𝑀1=𝑅1𝐴+𝐴1cos𝛽2+𝑖𝐴𝐴1sin𝛽+𝑅2𝐴𝐴1cos𝛽+𝑖𝐴+𝐴1,𝑀sin𝛽2=𝐵+𝐵1cos𝛾2+𝑖𝐵𝐵1sin𝛾,(4.10) where 𝐴 and 𝐵 are nonzero real parameters to be determined. Left-hand side of (4.9) is a real quantity. Thus the imaginary part of right-hand side of (4.9) must be zero. Thus we obtain 𝑅1+𝑅2𝑅𝐴1𝑅2𝐴1=0,𝐵𝐵1=0.(4.11) On solving, we get 𝐴=(𝑅1𝑅2)/(𝑅1+𝑅2), 𝐴1=(𝑅1+𝑅2)/(𝑅1𝑅2), and 𝐵=𝐵1=1, provided that 𝑅1±𝑅2>0,

Substituting the values of 𝐴, 𝐴1, 𝐵 and 𝐵1 in (4.10), we obtain𝑀1𝑅=21+𝑅21𝑅222sin2𝛽21,𝑀2=4sin2𝛾2.(4.12)

Let 𝑀1=2𝑀3 and 𝑀2=4𝑀4, where 𝑀3=𝑅1+𝑅21𝑅222sin2𝛽21,𝑀4=sin2𝛾2.(4.13)

Substituting the values of 𝑀1 and 𝑀2 in (4.9), we obtain sin2𝜙2=𝑝2𝑀3+2𝑀4(2/3)𝑀3𝑀421(1/6)𝑀31(1/3)𝑀4.(4.14)

Since sin2(𝜙/2)1, hence the method (4.8) is stable as long as𝑝max2𝑀3+2𝑀423𝑀3𝑀421min16𝑀3113𝑀4.(4.15) this implies 0<𝑝2𝑅42+𝑅22𝑅3232+𝑅2+𝑅22𝑅32,(4.16) which is the required stability interval for the scheme (4.8).

In order to facilitate the computation, we may rewrite (4.8) in operator split form (see [18, 19]) as 𝛿1+𝑦2𝑢12𝑙,𝑚=𝑝2𝑅1𝛿𝑥2+𝑅22𝜇𝑥𝛿𝑥+𝛿𝑦2+16𝑅1𝛿𝑥2+𝑅22𝜇𝑥𝛿𝑥𝛿𝑦2𝑢𝑗𝑙,𝑚+1𝑓,1+𝑅121𝛿𝑥2+𝑅22𝜇𝑥𝛿𝑥𝛿𝑡2𝑢𝑗𝑙,𝑚=𝑢𝑙,𝑚(4.17) or,11+𝑅121𝛿𝑥2+𝑅22𝜇𝑥𝛿𝑥𝑢𝑗+1𝑙,𝑚=𝑢𝑙,𝑚+11+𝑅121𝛿2𝑥+𝑅22𝜇𝑥𝛿𝑥2𝑢𝑗𝑙,𝑚𝑢𝑗1𝑙,𝑚,(4.18) where 𝑢𝑙,𝑚 is an intermediate value. The intermediate boundary values required for solving (4.17) can be obtained from (4.18).

The left-hand sides of (4.17) and (4.18) are factorizations into 𝑦- and 𝑥-differences, respectively, which allows us to solve by sweeping first (4.17) in the 𝑦- and then (4.18) in the 𝑥-direction. It will be seen that these sweeps require only the solution of tri-diagonal systems.

5. Super Stable Method for Two-Dimensional Linear Hyperbolic Equation

Let us consider the equation of the following form:𝑢𝑡𝑡+2𝛼𝑢𝑡+𝛽2𝑢=𝑢𝑥𝑥+𝑢𝑦𝑦+𝑓(𝑥,𝑦,𝑡),0<𝑥,𝑦<1,𝑡>0,𝛼>0,𝛽0(5.1) subject to appropriate initial and Dirichlet boundary conditions that are prescribed. The equation above represents two-space dimensional linear telegraphic hyperbolic equation. Applying the difference method (2.8) to the differential equation (5.1), we obtain𝛿2𝑡𝑢𝑗𝑙,𝑚+𝑎2𝜇𝑡𝛿𝑡𝑢𝑗𝑙,𝑚+𝑎𝛿122𝑥+𝛿2𝑦2𝜇𝑡𝛿𝑡𝑢𝑗𝑙,𝑚+𝑏𝑢𝑗𝑙,𝑚+𝑏12𝑝2𝛿2𝑥+𝛿2𝑦𝑢𝑗𝑙,𝑚𝑝26𝛿2𝑥𝛿2𝑦𝑈𝑗𝑙,𝑚+1𝛿122𝑥+𝛿2𝑦𝛿2𝑡𝑈𝑗𝑙,𝑚=𝐹𝐹,(5.2) where𝑎=𝛼2𝑘2,𝑏=𝛽2𝑘2𝑘,𝐹𝐹=2𝑓12𝑗𝑙+1,𝑚+𝑓𝑗𝑙1,𝑚+𝑓𝑗𝑙,𝑚+1+𝑓𝑗𝑙,𝑚1+8𝑓𝑗𝑙,𝑚.(5.3)

The previous scheme is conditionally stable (see [7]). In order to obtain a superstable method, we simply follow the ideas given by Chawla [24] and the scheme (5.2) may be re-written as1+𝜂𝑏2+112𝛾𝑝2𝛿2𝑥+𝛿2𝑦𝛿2𝑡𝑢𝑗𝑙,𝑚+𝑎𝛿1+2𝑥+𝛿122𝑦122𝜇𝑡𝛿𝑡𝑢𝑗𝑙,𝑚=𝑝2𝑏𝛿122𝑥+𝛿2𝑦+𝑝26𝛿2𝑥𝛿2𝑦𝑢𝑏𝑗𝑙,𝑚+𝐹𝐹,(5.4) where 𝜂 and 𝛾 are free parameters to be determined. The additional terms are of high orders and do not affect the accuracy of the scheme.

To study stability for the scheme (5.4), we substitute 𝑢𝑗𝑙,𝑚=𝜉𝑗𝑒𝑖𝜃𝑙𝑒𝑖𝜙𝑚 in the homogeneous part of the (5.4), and we obtain the characteristic equation:𝐴𝜉2+𝐵𝜉+𝐶=0,(5.5) where𝐴=1+𝜂𝑏2+𝑎𝑎3sin2𝜃2+sin2𝜙2+4𝛾𝑝213sin2𝜃2+sin2𝜙2,𝐵=22𝜂𝑏2𝑝+42𝑏12sin2𝜃2+sin2𝜙2+2+𝑏38𝛾𝑝2sin2𝜃2+sin2𝜙28𝑝23sin2𝜃2sin2𝜙2,𝐶=1+𝜂𝑏2𝑎+𝑎3sin2𝜃2+sin2𝜙2+4𝛾𝑝213sin2𝜃2+sin2𝜙2.(5.6)

Using the transformation 𝜉=(1+𝑧)/(1𝑧) in (5.5), we obtain the transformed characteristic equation:(𝐴𝐵+𝐶)𝑧2+2(𝐴𝐶)𝑧+(𝐴+𝐵+𝐶)=0.(5.7) For stability of the scheme (5.4), we must have the following conditions: 𝐴+𝐵+𝐶>0,𝐴𝐶>0,𝐴𝐵+𝐶>0.(5.8) It is easy to verify that the coefficient is as follows:𝑏𝐴+𝐵+𝐶=6sin2𝜃2+sin2𝜙2+𝑏2cos2𝜃2+cos2𝜙2+4𝑝2sin2𝜃2+sin2𝜙22sin2𝜃2sin2𝜙2+16𝑝23sin2𝜃2sin2𝜙2>0for𝛼>0,𝛽0,andall𝜃,𝜙,except𝜃=𝜙=0or2𝜋.(5.9) We can treat this case separately.

The coefficient is as follows:𝐴𝐶=𝑎13sin2𝜃2+sin2𝜙2+cos2𝜃2+cos2𝜙2>0for𝛼>0,𝛽0,andall𝜃,𝜙.(5.10) For stability, it is required that the third coefficient is as follows:𝐴𝐵+𝐶=4+4𝜂𝑏2𝑏𝑏+3sin2𝜃2+sin2𝜙2+4(4𝛾1)𝑝213sin2𝜃2+sin2𝜙2+8𝑝23sin2𝜃2sin2𝜙2>0.(5.11) Multiplying throughout of (5.11) by 16𝜂, we get (64𝜂1)+(8𝜂𝑏1)2+16𝑏𝜂3sin2𝜃2+sin2𝜙2+64𝜂(4𝛾1)𝑝213sin2𝜃2+sin2𝜙2+128𝑝23sin2𝜃2sin2𝜙2>0.(5.12) Thus the scheme is stable if1𝜂64,𝛾1+3𝑝212𝑝2.(5.13) For 𝜃=𝜙=0 or 2𝜋 and 𝛽=0, we have the characteristic equation:1+𝑎𝜉22𝜉+1𝑎=0(5.14) whose roots are 𝜉1,2=1, (1𝑎)/(1+𝑎)=(1𝛼𝑘)/(1+𝛼𝑘).

In this case also |𝜉|1 and the scheme (5.4) is stable.

Hence for 𝛼>0, 𝛽0, 𝜂1/64, 𝛾(1+3𝑝2)/12𝑝2 (which are independent of and 𝑘), the scheme (5.4) is superstable.

Further, the scheme (5.4) in product form may be written as11+12𝛾𝑝2𝛿2𝑦1+𝜂𝑏2+112𝛾𝑝2𝛿2𝑥𝛿2𝑡𝑢𝑗𝑙,𝑚+𝑎𝛿1+2𝑥122𝜇𝑡𝛿𝑡𝑢𝑗𝑙,𝑚=𝑝2𝑏𝛿122𝑥+𝛿2𝑦+𝑝26𝛿2𝑥𝛿2𝑦𝑢𝑏𝑗𝑙,𝑚+𝐹𝐹.(5.15)

In this case also additional terms are of high orders, which do not affect the accuracy of the scheme. In order to facilitate the computation, we may rewrite the scheme (5.15) in two-step operator split form:11+12𝛾𝑝2𝛿2𝑦𝑢𝑙,𝑚=𝑝2𝑏𝛿122𝑥+𝛿2𝑦+𝑝26𝛿2𝑥𝛿2𝑦𝑢𝑏𝑗𝑙,𝑚+𝐹𝐹,(5.16a)1+𝜂𝑏2+112𝛾𝑝2𝛿2𝑥𝛿2𝑡𝑢𝑗𝑙,𝑚+𝑎𝛿1+2𝑥122𝜇𝑡𝛿𝑡𝑢𝑗𝑙,𝑚=𝑢𝑙,𝑚,(5.16b)where 𝑢𝑙,𝑚 is any intermediate value, and the intermediate boundary conditions required for the solution of 𝑢𝑙,𝑚 may be obtained from (5.16b). The left-hand side matrices represented by (5.16a) and (5.16b) are tridiagonal, thus very easily solved in the region 0<𝑥, 𝑦<1, 𝑡>0 using a tridiagonal solver.

6. Numerical Illustrations

In this section, we have solved some benchmark problems using the method described by (2.8) and compared our results with the results of the fourth-order numerical method discussed in [6, 7]. The exact solutions are provided in each case. The initial and boundary conditions may be obtained using the exact solution as a test procedure. The linear difference equations have been solved using a tridiagonal solver, whereas non-linear difference equations have been solved using the Newton-Raphson method. While using the Newton-Raphson method, the iterations were stopped when absolute error tolerance ≤10−12 was achieved. In order to demonstrate the fourth-order convergence of the proposed method, for the computation of Example 6.1, we have chosen the fixed value of the parameter 𝜎=𝑘/2=3.2 for Example 6.1, and for other examples we have chosen the value of 𝜎=1.6. For this choice, our method behaves like a fourth-order method in space. All computations were carried out using double precision arithmetic.

Note that the proposed method (2.8) for second-order hyperbolic equations is a three-level scheme. The value of 𝑢 at 𝑡=0 is known from the initial condition. To start any computation, it is necessary to know the numerical value of 𝑢 of required accuracy at 𝑡=𝑘. In this section, we discuss an explicit scheme of 𝑂(𝑘2) for 𝑢 at first time level, that is, at 𝑡=𝑘 in order to solve the differential equation (1.1) using the method (2.8), which is applicable to problems in both Cartesian and polar coordinates.

Since the values of 𝑢 and 𝑢𝑡 are known explicitly at 𝑡=0, this implies that all their successive tangential derivatives are known at 𝑡=0; that is, the values of 𝑢,𝑢𝑥,𝑢𝑥𝑥,𝑢𝑦,𝑢𝑦𝑦,,𝑢𝑡,𝑢𝑡𝑥,𝑢𝑡𝑦,, and so forth are known at 𝑡=0.

An approximation for 𝑢 of 𝑂(𝑘2) at 𝑡=𝑘 may be written as𝑢1𝑙,𝑚=𝑢0𝑙,𝑚+𝑘𝑢0𝑡𝑙,𝑚+𝑘22𝑢𝑡𝑡0𝑙,𝑚𝑘+𝑂3.(6.1)

From (1.1), we have𝑢𝑡𝑡0𝑙,𝑚=𝐴(𝑥,𝑦,𝑡,𝑢)𝑢𝑥𝑥+𝐵(𝑥,𝑦,𝑡,𝑢)𝑢𝑦𝑦+𝐺𝑥,𝑦,𝑡,𝑢,𝑢𝑥,𝑢𝑦,𝑢𝑡0𝑙,𝑚.(6.2)

Thus using the initial values and their successive tangential derivative values, from (6.2) we can obtain the value of (𝑢𝑡𝑡)0𝑙,𝑚, and then ultimately, from (6.1) we can compute the value of 𝑢 at first time level, that is, at 𝑡=𝑘.

Example 6.1 (the telegraphic equation with forcing function). 𝑢𝑡𝑡+2𝛼𝑢𝑡+𝛽2𝑢=𝑢𝑥𝑥+𝑢𝑦𝑦+24𝛼+𝛽2𝑒2𝑡sinh𝑥sinh𝑦,0<𝑥,𝑦<1,𝑡>0,(6.3) where 𝛼>0 and 𝛽0 are real parameters. The exact solution is 𝑢=𝑒2𝑡sinh𝑥sinh𝑦. The maximum absolute errors and CPU time (in seconds) are tabulated in Table 1 at 𝑡=5 for 𝛼>0 and 𝛽0.

Example 6.2 (wave equation in cylindrical polar coordinates). 𝑢𝑡𝑡=𝑢𝑟𝑟+𝑢𝑧𝑧+1𝑟𝑢𝑟3cosh𝑟+sinh𝑟𝑟cosh𝑧sin𝑡,0<𝑟,𝑧<1,𝑡>0.(6.4) The aforementioned equation represents the two-space dimensional wave equation in cylindrical polar coordinates. The exact solution is 𝑢=cosh𝑟cosh𝑧sin𝑡. The maximum absolute errors and CPU time (in seconds) are tabulated in Table 2 at 𝑡=1 and 𝑡=2.

Example 6.3 (Van der Pol type nonlinear wave equation). 𝑢𝑡𝑡=𝑢𝑥𝑥+𝑢𝑦𝑦𝑢+𝛾2𝑢1𝑡+2𝜋2+𝛾2𝑒2𝛾𝑡sin2(𝜋𝑥)sin2𝑒(𝜋𝑦)𝛾𝑡sin(𝜋𝑥)sin(𝜋𝑦),0<𝑥,𝑦<1,𝑡>0,(6.5) with exact solution 𝑢=𝑒𝛾𝑡sin(𝜋𝑥)sin(𝜋𝑦). The maximum absolute errors and CPU time (in seconds) are tabulated in Table 3 at 𝑡=2 for 𝛾=1, 2, and 3.

Example 6.4 (dissipative nonlinear wave equation). 𝑢𝑡𝑡=𝑢𝑥𝑥+𝑢𝑦𝑦2𝑢𝑢𝑡+2sin(𝜋𝑥)sin(𝜋𝑦)cos𝑡+2𝜋21sin(𝜋𝑥)sin(𝜋𝑦)sin𝑡,0<𝑥,𝑦<1,𝑡>0,(6.6) with exact solution 𝑢=sin(𝜋𝑥)sin(𝜋𝑦)sin𝑡. The maximum absolute errors and CPU time (in seconds) are tabulated in Table 4 at 𝑡=1.

Example 6.5 (nonlinear wave equation with variable coefficients). 𝑢𝑡𝑡=1+𝑥2𝑢𝑥𝑥+1+𝑦2𝑢𝑦𝑦𝑢+𝛾𝑢𝑥+𝑢𝑦+𝑢𝑡+𝛾𝑒𝑡(cosh𝑥cosh𝑦sinh(𝑥+𝑦))1𝑥2𝑦2𝑒𝑡cosh𝑥cosh𝑦,0<𝑥,𝑦<1,𝑡>0,(6.7) with exact solution 𝑢=𝑒𝑡cosh𝑥cosh𝑦. The maximum absolute errors and CPU time (in seconds) are tabulated in Table 5 at 𝑡=1 for 𝛾=1, 2 and 5.

Example 6.6 (quasilinear hyperbolic equation). 𝑢𝑡𝑡=1+𝑢2𝑢𝑥𝑥+𝑢𝑦𝑦𝑢+𝛾𝑢𝑥+𝑢𝑦+𝑢𝑡+𝛾𝑒𝑡𝑒(sinh𝑥cosh𝑦cosh(𝑥+𝑦))1𝑡𝑒sinh𝑥cosh𝑦2𝑡sinh𝑥cosh𝑦3,0<𝑥,𝑦<1,𝑡>0,(6.8) with exact solution 𝑢=𝑒𝑡sinh𝑥cosh𝑦. The maximum absolute errors and CPU time (in seconds) are tabulated in Table 6 at 𝑡=1 for 𝛾=1, 2, and 5.
The order of convergence may be obtained by using the following formula: 𝑒log1𝑒log2log1log2,(6.9) where 𝑒1and 𝑒2 are maximum absolute errors for two uniform mesh widths 1 and 2, respectively. For computation of order of convergence of the proposed method, we have considered errors for last two values of , that is, 1=1/32, 2=1/64 for two linear problems, and 1=1/16, 2=1/32 for all non-linear and quasilinear problems, and results are reported in Table 7.

7. Concluding Remarks

Available numerical methods for the numerical solution of second-order quasilinear wave equations are of order four, which require 19-grid points. In this article, using the same number of grid points and five evaluations of the function 𝑔 (as compared to seven evaluations of the function 𝑔 discussed in [6, 7]), we have derived a new stable numerical method of 𝑂(𝑘2+𝑘22+4) accuracy for the solution of quasilinear wave equation (1.1). Ultimately, we use less algebra for computation, and for a fixed parameter value 𝜎=𝑘/2, the proposed method behaves like a fourth-order method, which is exhibited from the computed results. The proposed numerical method is applicable to wave equation in polar coordinates, and for the damped wave equation and telegraphic equation the method is shown to be unconditionally stable.

Acknowledgments

The authors thank the anonymous reviewer for his valuable suggestions, which substantially improved the standard of the paper. This research was supported by The University of Delhi under research Grant no. Dean (R)/R & D/2010/1311.