A New High-Order Approximation for the Solution of Two-Space-Dimensional Quasilinear Hyperbolic Equations

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,𝑦,𝑡)=𝑎₁(𝑦,𝑡),0≤𝑦≤1,𝑡≥0,(1.3a)(𝑥,𝑡),𝑢(𝑥,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 𝑢(𝑥,𝑦,𝑡)∈𝐶⁶, 𝐴(𝑥,𝑦,𝑡,𝑢),𝐵(𝑥,𝑦,𝑡,𝑢)∈𝐶⁴, 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 [8–16]) 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 𝐴^𝑗𝑥_𝑙,𝑚=𝜕𝐴^{𝑗𝑙,𝑚}/𝜕𝑥,𝐴^𝑗𝑦_𝑙,𝑚=𝜕𝐴^{𝑗𝑙,𝑚}/𝜕𝑦,𝐴^𝑗𝑥𝑥_𝑙,𝑚=𝜕²𝐴^{𝑗𝑙,𝑚}/𝜕𝑥²,… 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.2d)𝑈^𝑗𝑡𝑡_𝑙±1,𝑚=𝑈^{𝑗+1𝑙±1,𝑚}−2𝑈^{𝑗𝑙±1,𝑚}+𝑈^{𝑗−1𝑙±1,𝑚}𝑘²,(2.2e)𝑈^𝑗𝑡𝑡_𝑙,𝑚±1=𝑈^{𝑗+1𝑙,𝑚±1}−2𝑈^{𝑗𝑙,𝑚±1}+𝑈^{𝑗−1𝑙,𝑚±1}𝑘²,(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.3d)𝑈^𝑗𝑥𝑥_𝑙,𝑚±1=𝑈^{𝑗𝑙+1,𝑚±1}−2𝑈^{𝑗𝑙,𝑚±1}+𝑈^{𝑗𝑙−1,𝑚±1}ℎ²,(2.3e) 𝑈^𝑗𝑦_𝑙,𝑚=𝑈^{𝑗𝑙,𝑚+1}−𝑈^{𝑗𝑙,𝑚−1}(2ℎ),(2.4a)𝑈^𝑗𝑦_𝑙±1,𝑚=𝑈^{𝑗𝑙±1,𝑚+1}−𝑈^{𝑗𝑙±1,𝑚−1}(2ℎ),(2.4b)𝑈^𝑗𝑦_𝑙,𝑚±1=±3𝑈^{𝑗𝑙,𝑚±1}∓4𝑈^{𝑗𝑙,𝑚}±𝑈^{𝑗𝑙,𝑚∓1}(2ℎ),(2.4c)𝑈^𝑗𝑦𝑦_𝑙,𝑚=𝑈^{𝑗𝑙,𝑚+1}−2𝑈^{𝑗𝑙,𝑚}+𝑈^{𝑗𝑙,𝑚−1}ℎ²,(2.4d)𝑈^𝑗𝑦𝑦_𝑙±1,𝑚=𝑈^{𝑗𝑙±1,𝑚+1}−2𝑈^{𝑗𝑙±1,𝑚}+𝑈^{𝑗𝑙±1,𝑚−1}ℎ²,(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,𝑚}+ℎ²8𝐴^{𝑗𝑙,𝑚}𝐴^{𝑗𝑥𝑙,𝑚}𝑈^𝑗𝑥𝑥_𝑙,𝑚+𝐵^𝑗𝑥_𝑙,𝑚𝑈^𝑗𝑦𝑦_𝑙,𝑚,(2.6a)𝑈^𝑗𝑦_𝑙,𝑚=𝑈^𝑗𝑦_𝑙,𝑚+ℎ16𝐵^{𝑗𝑙,𝑚}𝐺^{𝑗𝑙,𝑚+1}−𝐺^{𝑗𝑙,𝑚−1}−ℎ16𝐵^{𝑗𝑙,𝑚}+ℎ𝐴_{𝑗𝑙,𝑚}𝑈^𝑗𝑡𝑡_𝑙,𝑚+1−𝑈^𝑗𝑡𝑡_{𝑙,𝑚−1}16𝐵^{𝑗𝑙,𝑚}𝑈^𝑗𝑥𝑥_𝑙,𝑚+1−𝑈^𝑗𝑥𝑥_{𝑙,𝑚−1}+ℎ²8𝐵^{𝑗𝑙,𝑚}𝐴^{𝑗𝑦𝑙,𝑚}𝑈^𝑗𝑥𝑥_𝑙,𝑚+𝐵^𝑗𝑦_𝑙,𝑚𝑈^𝑗𝑦𝑦_𝑙,𝑚,(2.6b) 𝐺^{𝑗𝑙,𝑚}=𝑔𝑥_𝑙,𝑦_𝑚,𝑡_𝑗,𝑈^{𝑗𝑙,𝑚},𝑈^𝑗𝑥_𝑙,𝑚,𝑈^𝑗𝑦_𝑙,𝑚,𝑈^𝑗𝑡_𝑙,𝑚.(2.7)

Then at each grid point (𝑥_𝑙,𝑦_𝑚,𝑡_𝑗), 𝑙,𝑚=1,2,…,𝑁, 𝑗=2,…,𝐽, an approximation with accuracy of 𝑂(𝑘²+𝑘²ℎ²+ℎ⁴) for the solution of differential equation (2.1) may be written as𝐿_𝑢≡𝑝²𝐴^{𝑗𝑙,𝑚}−ℎ²6𝐴𝐴^{𝑗𝑥𝑙,𝑚𝑗𝑙,𝑚}𝐴^{𝑗𝑥𝑙,𝑚}−ℎ²6𝐵𝐵^{𝑗𝑦𝑙,𝑚𝑗𝑙,𝑚}𝐴^{𝑗𝑦𝑙,𝑚}+ℎ²12𝐴^{𝑗𝑥𝑥𝑙,𝑚}+ℎ²12𝐴^{𝑗𝑦𝑦𝑙,𝑚}𝛿^2𝑥𝑈^{𝑗𝑙,𝑚}+𝑝²𝐵^{𝑗𝑙,𝑚}−ℎ²6𝐴𝐴^{𝑗𝑥𝑙,𝑚𝑗𝑙,𝑚}𝐵^{𝑗𝑥𝑙,𝑚}−ℎ²6𝐵𝐵^{𝑗𝑦𝑙,𝑚𝑗𝑙,𝑚}𝐵^{𝑗𝑦𝑙,𝑚}+ℎ²12𝐵^{𝑗𝑥𝑥𝑙,𝑚}+ℎ²12𝐵^{𝑗𝑦𝑦𝑙,𝑚}𝛿^2𝑦𝑈^{𝑗𝑙,𝑚}+ℎ𝑝²12𝐴^{𝑗𝑦𝑙,𝑚}−𝐵^{𝑗𝑦𝑙,𝑚}𝐵^{𝑗𝑙,𝑚}𝐴^{𝑗𝑙,𝑚}𝛿^2𝑥2𝜇_𝑦𝛿_𝑦𝑈^{𝑗𝑙,𝑚}+ℎ𝑝²12𝐵^{𝑗𝑥𝑙,𝑚}−𝐴^{𝑗𝑥𝑙,𝑚}𝐴^{𝑗𝑙,𝑚}𝐵^{𝑗𝑙,𝑚}𝛿^2𝑦2𝜇_𝑥𝛿_𝑥𝑈^{𝑗𝑙,𝑚}+𝑝²12𝐴^{𝑗𝑙,𝑚}+𝐵^{𝑗𝑙,𝑚}𝛿^2𝑥𝛿^2𝑦𝑈^{𝑗𝑙,𝑚}=𝑘²121−ℎ𝐴𝐴^{𝑗𝑥𝑙,𝑚𝑗𝑙,𝑚}𝑈^𝑗𝑡𝑡_𝑙+1,𝑚+1+ℎ𝐴^{𝑗𝑥𝑙,𝑚}𝐴^{𝑗𝑙,𝑚}𝑈^𝑗𝑡𝑡_{𝑙−1,𝑚}+1−ℎ𝐵^{𝑗𝑦𝑙,𝑚}𝐵^{𝑗𝑙,𝑚}𝑈^𝑗𝑡𝑡_𝑙,𝑚+1+1+ℎ𝐵^{𝑗𝑦𝑙,𝑚}𝐵^{𝑗𝑙,𝑚}𝑈^𝑗𝑡𝑡_{𝑙,𝑚−1}+8𝑈^𝑗𝑡𝑡_𝑙,𝑚−𝑘₂121−ℎ𝐴𝐴^{𝑗𝑥𝑙,𝑚𝑗𝑙,𝑚}𝐺^{𝑗𝑙+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 𝑇^{𝑗𝑙,𝑚}=𝑂(𝑘⁴+𝑘⁴ℎ²+𝑘²ℎ⁴).

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𝐿_𝑢=𝑘²121−ℎ𝐴𝐴^{𝑗𝑥𝑙,𝑚𝑗𝑙,𝑚}𝑈^𝑗𝑡𝑡_𝑙+1,𝑚+1+ℎ𝐴^{𝑗𝑥𝑙,𝑚}𝐴^{𝑗𝑙,𝑚}𝑈^𝑗𝑡𝑡_{𝑙−1,𝑚}+1−ℎ𝐵^{𝑗𝑦𝑙,𝑚}𝐵^{𝑗𝑙,𝑚}𝑈^𝑗𝑡𝑡_𝑙,𝑚+1+1+ℎ𝐵^{𝑗𝑦𝑙,𝑚}𝐵^{𝑗𝑙,𝑚}𝑈^𝑗𝑡𝑡_{𝑙,𝑚−1}+8𝑈^𝑗𝑡𝑡_𝑙,𝑚−𝑘₂121−ℎ𝐴𝐴^{𝑗𝑥𝑙,𝑚𝑗𝑙,𝑚}𝐺^{𝑗𝑙+1,𝑚}+1+ℎ𝐴^{𝑗𝑥𝑙,𝑚}𝐴^{𝑗𝑙,𝑚}𝐺^{𝑗𝑙−1,𝑚}+1−ℎ𝐵^{𝑗𝑦𝑙,𝑚}𝐵^{𝑗𝑙,𝑚}𝐺^{𝑗𝑙,𝑚+1}+1+ℎ𝐵^{𝑗𝑦𝑙,𝑚}𝐵^{𝑗𝑙,𝑚}𝐺^{𝑗𝑙,𝑚−1}+8𝐺^{𝑗𝑙,𝑚}+𝑂𝑘⁴+𝑘⁴ℎ²+𝑘²ℎ⁴.(3.3)

With the help of the approximations (2.2a)–(2.2e), and from (2.5a) and (2.5b), we obtain𝐺^{𝑗𝑙±1,𝑚}=𝐺^{𝑗𝑙±1,𝑚}−ℎ²3𝜕³𝑈^{𝑗𝑙,𝑚}𝜕𝑥³𝛼^{𝑗𝑙,𝑚}+ℎ²6𝜕³𝑈^{𝑗𝑙,𝑚}𝜕𝑦³𝛽^{𝑗𝑙,𝑚}+𝑂𝑘²±𝑘²ℎ±ℎ³+ℎ⁴,𝐺^{𝑗𝑙,𝑚±1}=𝐺^{𝑗𝑙,𝑚±1}+ℎ²6𝜕³𝑈^{𝑗𝑙,𝑚}𝜕𝑥³𝛼^{𝑗𝑙,𝑚}−ℎ²3𝜕³𝑈^{𝑗𝑙,𝑚}𝜕𝑦³𝛽^{𝑗𝑙,𝑚}+𝑂𝑘²±𝑘²ℎ±ℎ³+ℎ⁴.(3.4) Now we define the following approximations: 𝑈^𝑗𝑥_𝑙,𝑚=+𝑈𝑎^𝑗𝑥_𝑙,𝑚13+𝑎₁₁ℎ𝐺^{𝑗𝑙+1,𝑚}−𝐺^{𝑗𝑙−1,𝑚}+𝑎¹²ℎ𝑈^𝑗𝑡𝑡_𝑙+1,𝑚−𝑈^𝑗𝑡𝑡_{𝑙−1,𝑚}ℎ𝑈^𝑗𝑦𝑦_𝑙+1,𝑚−𝑈^𝑗𝑦𝑦_{𝑙−1,𝑚}+ℎ²𝑎¹⁴𝑈^𝑗𝑥𝑥_𝑙,𝑚+𝑎₁₅𝑈^𝑗𝑦𝑦_𝑙,𝑚,𝑈^𝑗𝑦_𝑙,𝑚=+𝑈𝑏^𝑗𝑦_𝑙,𝑚13+𝑏₁₁ℎ𝐺^{𝑗𝑙,𝑚+1}−𝐺^{𝑗𝑙,𝑚−1}+𝑏¹²ℎ𝑈^𝑗𝑡𝑡_𝑙,𝑚+1−𝑈^𝑗𝑡𝑡_{𝑙,𝑚−1}ℎ𝑈^𝑗𝑥𝑥_𝑙,𝑚+1−𝑈^𝑗𝑥𝑥_{𝑙,𝑚−1}+ℎ²𝑏¹⁴𝑈^𝑗𝑥𝑥_𝑙,𝑚+𝑏₁₅𝑈^𝑗𝑦𝑦_𝑙,𝑚,(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𝑈^𝑗𝑥_𝑙,𝑚=𝑈^𝑗𝑥_𝑙,𝑚+ℎ²6𝑇¹+𝑂𝑘²+𝑘²ℎ²+ℎ⁴,𝑈^𝑗𝑦_𝑙,𝑚=𝑈^𝑗𝑦_𝑙,𝑚+ℎ²6𝑇²+𝑂𝑘²+𝑘²ℎ²+ℎ⁴,(3.6) where 𝑇₁=1−12𝑎₁₁𝐴^{𝑗𝑙,𝑚}𝑈^{𝑗𝑥𝑥𝑥𝑙,𝑚}+12𝑎₁₁+𝑎₁₂𝑈^{𝑗𝑡𝑡𝑥}_𝑙,𝑚+12−𝑎₁₁𝐵^{𝑗𝑙,𝑚}+𝑎₁₃𝑈^{𝑗𝑥𝑦𝑦𝑙,𝑚}𝑇₂+6𝑎₁₄−12𝑎₁₁𝐴^𝑗𝑥_𝑙,𝑚𝑈^{𝑗𝑥𝑥𝑙,𝑚}+6𝑎₁₅−12𝑎₁₁𝐵^𝑗𝑥_𝑙,𝑚𝑈^{𝑗𝑦𝑦𝑙,𝑚},=1−12𝑏₁₁𝐵^{𝑗𝑙,𝑚}𝑈^{𝑗𝑦𝑦𝑦𝑙,𝑚}+12𝑏₁₁+𝑏₁₂𝑈^{𝑗𝑡𝑡𝑦}_𝑙,𝑚+12−𝑏₁₁𝐴^{𝑗𝑙,𝑚}+𝑏₁₃𝑈^{𝑗𝑥𝑥𝑦𝑙,𝑚}+6𝑏₁₄−12𝑏₁₁𝐴^𝑗𝑦_𝑙,𝑚𝑈^{𝑗𝑥𝑥𝑙,𝑚}+6𝑏₁₅−12𝑏₁₁𝐵^𝑗𝑦_𝑙,𝑚𝑈^{𝑗𝑦𝑦𝑙,𝑚}.(3.7) Thus from (2.7), we obtain 𝐺^{𝑗𝑙,𝑚}=𝐺^{𝑗𝑙,𝑚}+ℎ²6𝑇¹𝛼^{𝑗𝑙,𝑚}+ℎ²6𝑇²𝛽^{𝑗𝑙,𝑚}+𝑂𝑘²+𝑘²ℎ²+ℎ⁴.(3.8) Finally by the help of the approximations (3.4) and (3.8), from (2.8) and (3.3), we get 𝑇^{𝑗𝑙,𝑚}=−𝑘²ℎ²36𝜕³𝑈^{𝑗𝑙,𝑚}𝜕𝑥³𝛼^{𝑗𝑙,𝑚}+𝜕³𝑈^{𝑗𝑙,𝑚}𝜕𝑦³𝛽^{𝑗𝑙,𝑚}−4𝑇₁𝛼^{𝑗𝑙,𝑚}+𝑇₂𝛽^{𝑗𝑙,𝑚}+𝑂𝑘⁴+𝑘⁴ℎ²+𝑘²ℎ⁴.(3.9)

In order to obtain the difference approximation of 𝑂(𝑘²+𝑘²ℎ²+ℎ⁴) the coefficient of 𝑘²ℎ² in (3.9) must be zero, which gives the values of parameters:𝑎₁₁=116𝐴^{𝑗𝑙,𝑚},𝑎₁₂=−116𝐴^{𝑗𝑙,𝑚},𝑎₁₃=𝐵^{𝑗𝑙,𝑚}16𝐴^{𝑗𝑙,𝑚},𝑎₁₄=𝐴^{𝑗𝑥𝑙,𝑚}8𝐴^{𝑗𝑙,𝑚},𝑎₁₅=𝐵^{𝑗𝑥𝑙,𝑚}8𝐴^{𝑗𝑙,𝑚}𝑏₁₁,=116𝐵^{𝑗𝑙,𝑚},𝑏₁₂=−116𝐵^{𝑗𝑙,𝑚},𝑏₁₃=𝐴^{𝑗𝑙,𝑚}16𝐵^{𝑗𝑙,𝑚},𝑏₁₄=𝐴^{𝑗𝑦𝑙,𝑚}8𝐵^{𝑗𝑙,𝑚},𝑏₁₅=𝐵^{𝑗𝑦𝑙,𝑚}8𝐵^{𝑗𝑙,𝑚}.(3.10)

Thus we obtain the difference method of 𝑂(𝑘²+𝑘²ℎ²+ℎ⁴) for the differential equation (2.1) and the local truncation error reduces to 𝑇^{𝑗𝑙,𝑚}=𝑂(𝑘⁴+𝑘⁴ℎ²+𝑘²ℎ⁴).

Now we consider the numerical method of 𝑂(𝑘²+𝑘²ℎ²+ℎ⁴) 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ℎ)+𝑂ℎ²,=𝐴^{𝑗𝑙+1,𝑚}−2𝐴^{𝑗𝑙,𝑚}+𝐴^{𝑗𝑙−1,𝑚}ℎ²+𝑂ℎ²𝐴_{𝑗𝑦𝑙,𝑚},=𝐴^{𝑗𝑙,𝑚+1}−𝐴^{𝑗𝑙,𝑚−1}𝐴_{𝑗𝑦𝑦𝑙,𝑚}(2ℎ)+𝑂ℎ²,=𝐴^{𝑗𝑙,𝑚+1}−2𝐴^{𝑗𝑙,𝑚}+𝐴^{𝑗𝑙,𝑚−1}ℎ²+𝑂ℎ²𝐵_{𝑗𝑥𝑙,𝑚},=𝐵^{𝑗𝑙+1,𝑚}−𝐵^{𝑗𝑙−1,𝑚}𝐵_{𝑗𝑥𝑥𝑙,𝑚}(2ℎ)+𝑂ℎ²,=𝐵^{𝑗𝑙+1,𝑚}−2𝐵^{𝑗𝑙,𝑚}+𝐵^{𝑗𝑙−1,𝑚}ℎ²+𝑂ℎ²𝐵_{𝑗𝑦𝑙,𝑚},=𝐵^{𝑗𝑙,𝑚+1}−𝐵^{𝑗𝑙,𝑚−1}𝐵_{𝑗𝑦𝑦𝑙,𝑚}(2ℎ)+𝑂ℎ²,=𝐵^{𝑗𝑙,𝑚+1}−2𝐵^{𝑗𝑙,𝑚}+𝐵^{𝑗𝑙,𝑚−1}ℎ²+𝑂ℎ²,(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 𝐴^{𝑗𝑙,𝑚}−ℎ²6𝐴𝐴^{𝑗𝑥𝑙,𝑚𝑗𝑙,𝑚}𝐴^{𝑗𝑥𝑙,𝑚}−ℎ²6𝐵𝐵^{𝑗𝑦𝑙,𝑚𝑗𝑙,𝑚}𝐴^{𝑗𝑦𝑙,𝑚}+ℎ²12𝐴^{𝑗𝑥𝑥𝑙,𝑚}+ℎ²12𝐴^{𝑗𝑦𝑦𝑙,𝑚}=𝐴^{𝑗𝑙,𝑚}−124𝐴^{𝑗𝑙,𝑚}𝐴^{𝑗𝑙+1,𝑚}−𝐴^{𝑗𝑙−1,𝑚}²−124𝐵^{𝑗𝑙,𝑚}+1𝐵^{𝑗𝑙,𝑚+1}−𝐵^{𝑗𝑙,𝑚−1}𝐴^{𝑗𝑙,𝑚+1}−𝐴^{𝑗𝑙,𝑚−1}12𝐴^{𝑗𝑙+1,𝑚}−2𝐴^{𝑗𝑙,𝑚}+𝐴^{𝑗𝑙−1,𝑚}+112𝐴^{𝑗𝑙,𝑚+1}𝐵_{𝑗𝑙,𝑚}−2𝐴^{𝑗𝑙,𝑚}+𝐴^{𝑗𝑙,𝑚−1}+𝑂ℎ⁴,−ℎ²6𝐴𝐴^{𝑗𝑥𝑙,𝑚𝑗𝑙,𝑚}𝐵^{𝑗𝑥𝑙,𝑚}−ℎ²6𝐵𝐵^{𝑗𝑦𝑙,𝑚𝑗𝑙,𝑚}𝐵^{𝑗𝑦𝑙,𝑚}+ℎ²12𝐵^{𝑗𝑥𝑥𝑙,𝑚}+ℎ²12𝐵^{𝑗𝑦𝑦𝑙,𝑚}=𝐵^{𝑗𝑙,𝑚}−124𝐵^{𝑗𝑙,𝑚}𝐵^{𝑗𝑙,𝑚+1}−𝐵^{𝑗𝑙,𝑚−1}²−124𝐴^{𝑗𝑙,𝑚}+1𝐴^{𝑗𝑙+1,𝑚}−𝐴^{𝑗𝑙−1,𝑚}𝐵^{𝑗𝑙+1,𝑚}−𝐵^{𝑗𝑙−1,𝑚}12𝐵^{𝑗𝑙+1,𝑚}−2𝐵^{𝑗𝑙,𝑚}+𝐵^{𝑗𝑙−1,𝑚}+112𝐵^{𝑗𝑙,𝑚+1}−2𝐵^{𝑗𝑙,𝑚}+𝐵^{𝑗𝑙,𝑚−1}+𝑂ℎ⁴,…,andsoforth.(3.13)

Thus substituting the central difference approximations (3.11) into (2.8), we obtain the required numerical method of 𝑂(𝑘²+𝑘²ℎ²+ℎ⁴) 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 [20–24]).

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𝑦+16𝛿^2𝑥𝛿^2𝑦𝑈^{𝑗𝑙,𝑚}=𝑘²12𝑈^𝑗𝑡𝑡_𝑙+1,𝑚+𝑈^𝑗𝑡𝑡_{𝑙−1,𝑚}+𝑈^𝑗𝑡𝑡_𝑙,𝑚+1+𝑈^𝑗𝑡𝑡_{𝑙,𝑚−1}+8𝑈^𝑗𝑡𝑡_𝑙,𝑚−𝑘₂12𝐺^{𝑗𝑙+1,𝑚}+𝐺^{𝑗𝑙−1,𝑚}+𝐺^{𝑗𝑙,𝑚+1}+𝐺^{𝑗𝑙,𝑚−1}+8𝐺^{𝑗𝑙,𝑚}+𝑂𝑘⁴𝑙=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 𝑂(𝑘²+𝑘²ℎ²+ℎ⁴) 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/𝑥²,…, 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=𝐷_𝑙±ℎ𝐷_𝑥𝑙+ℎ²2𝐷^𝑙𝑥𝑥𝑓_{𝑗𝑙±1,𝑚}±𝑂ℎ³,=𝑓^{𝑗𝑙,𝑚}±ℎ𝑓^𝑗𝑥_𝑙,𝑚+ℎ²2𝑓^{𝑗𝑥𝑥𝑙,𝑚}𝑓_{𝑗𝑙,𝑚±1}±𝑂ℎ³,=𝑓^{𝑗𝑙,𝑚}±ℎ𝑓^𝑗𝑦_𝑙,𝑚+ℎ²2𝑓^{𝑗𝑦𝑦𝑙,𝑚}±𝑂ℎ³.(4.4)

Using the approximations (4.4) in the scheme (4.2) and neglecting local truncation error, we get1+112𝛿^2𝑥+𝑅₀2𝜇_𝑥𝛿_𝑥+𝛿^2𝑦12𝛿^2𝑡𝑢^{𝑗𝑙,𝑚}=𝑝²𝑅¹𝛿^2𝑥+𝑅₂2𝜇_𝑥𝛿_𝑥+𝛿^2𝑦+𝑅⁰6𝛿^2𝑦2𝜇_𝑥𝛿_𝑥+16𝛿^2𝑥𝛿^2𝑦𝑢^{𝑗𝑙,𝑚}+𝑓,(4.5) where𝑅₀=ℎ𝐷^𝑙2,𝑅¹=1+ℎ²122𝐷^𝑥𝑙+𝐷_𝑙²,𝑅²=𝑅₀+ℎ³24𝐷^𝑥𝑥𝑙𝑓=𝑘₂+𝐷_𝑙𝐷_𝑥𝑙,1212𝑓^{𝑗𝑙,𝑚}+ℎ²𝑓^{𝑗𝑥𝑥𝑙,𝑚}+𝑓^𝑗𝑦𝑦_𝑙,𝑚+𝐷_𝑙𝑓^𝑗𝑥_𝑙,𝑚.(4.6) The previous scheme in product form may be written as1+112𝛿^2𝑥+𝑅₀2𝜇_𝑥𝛿_𝑥1+1𝛿2^2𝑦𝛿^2𝑡𝑢^{𝑗𝑙,𝑚}=𝑝²𝑅¹𝛿^2𝑥+𝑅₂2𝜇_𝑥𝛿_𝑥+𝛿^2𝑦+𝑅⁰6𝛿^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]) as1+112𝑅¹𝛿^2𝑥+𝑅₂2𝜇_𝑥𝛿_𝑥1+1𝛿2^2𝑦𝛿^2𝑡𝑢^{𝑗𝑙,𝑚}=𝑝²𝑅¹𝛿^2𝑥+𝑅₂2𝜇_𝑥𝛿_𝑥+𝛿^2𝑦+16𝑅¹𝛿^2𝑥+𝑅₂2𝜇_𝑥𝛿_𝑥𝛿^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+𝜉⁻¹=−4sin²𝜙2=𝑝²𝑀₁+𝑀₂+(1/6)𝑀₁𝑀₂1+(1/12)𝑀₁1+(1/12)𝑀₂,(4.9) where𝑀₁=𝑅₁𝐴+𝐴⁻¹cos𝑀𝛽−2+𝑖𝐴−𝐴⁻¹sin𝛽+𝑅₂𝐴−𝐴⁻¹cos𝛽+𝑖𝐴+𝐴⁻¹sin𝛽₂,=𝐵+𝐵⁻¹cos𝛾−2+𝑖𝐵−𝐵⁻¹sin𝛾,(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 𝑅₁+𝑅₂𝐴𝐵−−𝐵𝑅₁₋₁−𝑅₂𝐴⁻¹=0,=0.(4.11) On solving, we get 𝐴=√(𝑅₁−𝑅₂)/(𝑅₁+𝑅₂), 𝐴⁻¹=√(𝑅₁+𝑅₂)/(𝑅₁−𝑅₂), and 𝐵=𝐵⁻¹=1, provided that 𝑅₁±𝑅₂>0,

Substituting the values of 𝐴, 𝐴⁻¹, 𝐵 and 𝐵⁻¹ in (4.10), we obtain𝑀₁=−2𝑅₁+𝑅²¹−𝑅²²2sin²𝛽2−1,𝑀²=−4sin²𝛾2.(4.12)

Let 𝑀₁=−2𝑀₃ and 𝑀₂=−4𝑀₄, where 𝑀₃=𝑅₁+𝑅²¹−𝑅²²2sin²𝛽2−1,𝑀⁴=sin²𝛾2.(4.13)

Substituting the values of 𝑀₁ and 𝑀₂ in (4.9), we obtain sin²𝜙2=𝑝²𝑀₃+2𝑀₄−(2/3)𝑀₃𝑀₄21−(1/6)𝑀₃1−(1/3)𝑀₄.(4.14)

Since sin²(𝜙/2)≤1, hence the method (4.8) is stable as long asmax𝑝²𝑀³+2𝑀₄−23𝑀³𝑀₄≤min21−16𝑀³1−31𝑀⁴.(4.15) this implies 0<𝑝²≤4−𝑅²+𝑅₂₂−𝑅₃₂32+𝑅₂+𝑅₂₂−𝑅₃₂,(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+𝛿12^𝑦2𝑢^∗𝑙,𝑚=𝑝²𝑅¹𝛿_𝑥2+𝑅₂2𝜇_𝑥𝛿_𝑥+𝛿_𝑦2+16𝑅¹1+1𝛿_𝑥2+𝑅₂2𝜇_𝑥𝛿_𝑥𝛿_𝑦2𝑢^{𝑗𝑙,𝑚}+𝑓,12𝑅¹𝛿_𝑥2+𝑅₂2𝜇_𝑥𝛿_𝑥𝛿_𝑡2𝑢^{𝑗𝑙,𝑚}=𝑢^∗𝑙,𝑚(4.17) or,1+112𝑅¹𝛿_𝑥2+𝑅₂2𝜇_𝑥𝛿_𝑥𝑢^{𝑗+1𝑙,𝑚}=𝑢^∗𝑙,𝑚+1+112𝑅¹𝛿^2𝑥+𝑅₂2𝜇_𝑥𝛿_𝑥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𝛼𝑢_𝑡+𝛽²𝑢=𝑢_𝑥𝑥+𝑢_𝑦𝑦+𝑓(𝑥,𝑦,𝑡),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𝜇_𝑡𝛿_𝑡𝑢^{𝑗𝑙,𝑚}+√𝑎12𝛿^2𝑥+𝛿^2𝑦2𝜇_𝑡𝛿_𝑡𝑢^{𝑗𝑙,𝑚}+𝑏𝑢^{𝑗𝑙,𝑚}+𝑏12−𝑝²𝛿^2𝑥+𝛿^2𝑦𝑢^{𝑗𝑙,𝑚}−𝑝²6𝛿^2𝑥𝛿^2𝑦𝑈^{𝑗𝑙,𝑚}+112𝛿^2𝑥+𝛿^2𝑦𝛿^2𝑡𝑈^{𝑗𝑙,𝑚}=𝐹𝐹,(5.2) where𝑎=𝛼²𝑘²,𝑏=𝛽²𝑘²,𝐹𝐹=𝑘²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 as1+𝜂𝑏²+112−𝛾𝑝²𝛿^2𝑥+𝛿^2𝑦𝛿^2𝑡𝑢^{𝑗𝑙,𝑚}+√𝑎1+1𝛿2^2𝑥+1𝛿2^2𝑦2𝜇^𝑡𝛿_𝑡𝑢^{𝑗𝑙,𝑚}=𝑝²−𝑏12𝛿^2𝑥+𝛿^2𝑦+𝑝²6𝛿^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:𝐴𝜉²+𝐵𝜉+𝐶=0,(5.5) where𝐴=1+𝜂𝑏²+√𝑎−√𝑎3sin²𝜃2+sin²𝜙2+4𝛾𝑝²−13sin²𝜃2+sin²𝜙𝐵=−2−2𝜂𝑏₂2,+4𝑝²−𝑏12sin²𝜃2+sin²𝜙+22+𝑏3−8𝛾𝑝²sin²𝜃2+sin²𝜙2−83𝑝²sin²𝜃2sin²𝜙𝐶=1+𝜂𝑏₂2,−√𝑎+√𝑎3sin²𝜃2+sin²𝜙2+4𝛾𝑝²−13sin²𝜃2+sin²𝜙2.(5.6)

Using the transformation 𝜉=(1+𝑧)/(1−𝑧) in (5.5), we obtain the transformed characteristic equation:(𝐴−𝐵+𝐶)𝑧²+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:𝐴+𝐵+𝐶=𝑏6sin²𝜃2+sin²𝜙2+2𝑏cos²𝜃2+cos²𝜙2+4𝑝²sin²𝜃2+sin²𝜙2−2sin²𝜃2sin²𝜙+16𝑝₂23sin²𝜃2sin²𝜙2>0for𝛼>0,𝛽≥0,andall𝜃,𝜙,except𝜃=𝜙=0or2𝜋.(5.9) We can treat this case separately.

The coefficient is as follows:𝐴−𝐶=√𝑎13sin²𝜃2+sin²𝜙2+cos²𝜃2+cos²𝜙2>0for𝛼>0,𝛽≥0,andall𝜃,𝜙.(5.10) For stability, it is required that the third coefficient is as follows:𝐴−𝐵+𝐶=4+4𝜂𝑏²−𝑏+𝑏3sin²𝜃2+sin²𝜙+4(4𝛾−1)𝑝₂2−13sin²𝜃2+sin²𝜙2+83𝑝²sin²𝜃2sin²𝜙2>0.(5.11) Multiplying throughout of (5.11) by 16𝜂, we get (64𝜂−1)+(8𝜂𝑏−1)²+16𝑏𝜂3sin²𝜃2+sin²𝜙+64𝜂(4𝛾−1)𝑝₂2−13sin²𝜃2+sin²𝜙2+1238𝑝²sin²𝜃2sin²𝜙2>0.(5.12) Thus the scheme is stable if𝜂≥164,𝛾≥11+2𝑝3²𝑝².(5.13) For 𝜃=𝜙=0 or 2𝜋 and 𝛽=0, we have the characteristic equation:1+√𝑎𝜉²−2𝜉+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𝑝²)/12𝑝² (which are independent of ℎ and 𝑘), the scheme (5.4) is superstable.

Further, the scheme (5.4) in product form may be written as1+=112𝑝−₂𝛾𝑝²𝛿^2𝑦1+𝜂𝑏²+112−𝛾𝑝²𝛿^2𝑥𝛿^2𝑡𝑢^{𝑗𝑙,𝑚}+√𝑎1+1𝛿2^2𝑥2𝜇^𝑡𝛿_𝑡𝑢^{𝑗𝑙,𝑚}−𝑏12𝛿^2𝑥+𝛿^2𝑦+𝑝²6𝛿^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:1+1+11𝜂2𝑏−₂𝛾𝑝²𝛿^2𝑦𝑢^∗𝑙,𝑚=𝑝²−1𝑏2𝛿^2𝑥+𝛿^2𝑦+𝑝6²𝛿^2𝑥𝛿^2𝑦−𝑏𝑢^{𝑗𝑙,𝑚}+𝐹𝐹,(5.16a)+112−𝛾𝑝²𝛿^2𝑥𝛿^2𝑡𝑢^{𝑗𝑙,𝑚}+√𝑎1+𝛿^2𝑥122𝜇^𝑡𝛿_𝑡𝑢^{𝑗𝑙,𝑚}=𝑢^∗𝑙,𝑚,(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⁻¹² 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 𝜎=𝑘/ℎ²=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 𝑂(𝑘²) 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 𝑂(𝑘²) at 𝑡=𝑘 may be written as𝑢^1𝑙,𝑚=𝑢^0𝑙,𝑚+𝑘𝑢^0𝑡_𝑙,𝑚+𝑘²2𝑢^𝑡𝑡^0𝑙,𝑚+𝑂𝑘³.(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−4𝛼+𝛽²𝑒^−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.

Table 1 

Example 6.1: the maximum absolute errors (using proposed 𝑂(𝑘²+𝑘²ℎ²+ℎ⁴)-method) (with CPU time).

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.

Table 2 

Example 6.2: the maximum absolute errors (with CPU time).

Example 6.3 (Van der Pol type nonlinear wave equation). 𝑢_𝑡𝑡=𝑢_𝑥𝑥+𝑢_𝑦𝑦+𝛾𝑢²−1𝑢_𝑡+2𝜋²+𝛾²𝑒^−2𝛾𝑡sin²(𝜋𝑥)sin²(𝜋0𝑦<)𝑥𝑒,𝑦^−𝛾𝑡s<in1(,𝜋𝑡𝑥>)si0n,(𝜋𝑦),(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.

Table 3 

Example 6.3: the maximum absolute errors (with CPU time).

Example 6.4 (dissipative nonlinear wave equation). 𝑢_𝑡𝑡=𝑢_𝑥𝑥+𝑢_𝑦𝑦−2𝑢𝑢_𝑡+2sin(𝜋𝑥)sin(𝜋𝑦)cos𝑡+2𝜋²−1s0in<(𝜋𝑥𝑥,)𝑦si<n1(𝜋,𝑦)s𝑡in>𝑡,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.

Table 4 

Example 6.4: the maximum absolute errors (with CPU time).

Example 6.5 (nonlinear wave equation with variable coefficients). 𝑢_𝑡𝑡=+1𝛾+𝑒𝑥²_−𝑡𝑢_𝑥𝑥+1+𝑦²𝑢_𝑦𝑦+𝛾𝑢𝑢_𝑥+𝑢_𝑦+𝑢_𝑡(cosh𝑥cosh𝑦−sinh(𝑥+𝑦))−01<−𝑥𝑥,²𝑦−<𝑦²1,𝑒𝑡^−𝑡>co0s,h𝑥cosh𝑦,(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.

Table 5 

Example 6.5: The maximum absolute errors (with CPU time).

Example 6.6 (quasilinear hyperbolic equation). 𝑢_𝑡𝑡=+1𝛾+𝑒𝑢²_−𝑡𝑢_𝑥𝑥+𝑢_𝑦𝑦+𝛾𝑢𝑢_𝑥+𝑢_𝑦+𝑢_𝑡−2𝑒_−𝑡(sinh𝑥cosh𝑦−cosh(𝑥+𝑦))−1𝑒^−𝑡sinh𝑥cosh𝑦sinh𝑥cosh𝑦³,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: log𝑒_ℎ1−log𝑒_ℎ2logℎ₁−logℎ₂,(6.9) where 𝑒_ℎ1and 𝑒_ℎ2 are maximum absolute errors for two uniform mesh widths ℎ₁ and ℎ₂, respectively. For computation of order of convergence of the proposed method, we have considered errors for last two values of ℎ, that is, ℎ₁=1/32, ℎ₂=1/64 for two linear problems, and ℎ₁=1/16, ℎ₂=1/32 for all non-linear and quasilinear problems, and results are reported in Table 7.

Table 6 

Example 6.6: The maximum absolute errors (with CPU time).

Table 7 

Order of convergence.

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 𝑂(𝑘²+𝑘²ℎ²+ℎ⁴) accuracy for the solution of quasilinear wave equation (1.1). Ultimately, we use less algebra for computation, and for a fixed parameter value 𝜎=𝑘/ℎ², 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.