we propose a new high-order approximation for the solution of two-space-dimensional quasilinear hyperbolic partial differential equation of the form ,
, , subject to appropriate initial and Dirichlet boundary conditions , where
and 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.
We consider the following two-space dimensional quasilinear hyperbolic partial differential equation:
subject to the initial conditions
and the boundary conditionswhere (1.1) is assumed to satisfy the hyperbolicity condition and in the solution region . 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.  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.  on Cartesian meshes, and by Dumbser and Käser  on general two- and three-dimensional unstructured meshes. In 1996, Mohanty et al.  have developed a high-order numerical method for the solution of two-space dimensional second-order nonlinear hyperbolic equation. Later, Mohanty et al.  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  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:
Let and be the mesh spacing in the space and time directions, respectively. We replace the solution region by a set of grid points , where , , , , and , , and being positive integers and . Let 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:
Then at each grid point , , , an approximation with accuracy of for the solution of differential equation (2.1) may be written as
where , and so forth and .
3. Derivation Procedure of the Approximation
At the grid point , let us denote
Further at the grid point , the exact solution value satisfies
Using Taylor expansion about the grid point , first we obtain
With the help of the approximations (2.2a)–(2.2e), and from (2.5a) and (2.5b), we obtain
Now we define the following approximations:
where and , , are free parameters to be determined.
By using the approximations (2.2a)–(2.4e) and (3.4), from (3.5), we get
Thus from (2.7), we obtain
Finally by the help of the approximations (3.4) and (3.8), from (2.8) and (3.3), we get
In order to obtain the difference approximation of the coefficient of in (3.9) must be zero, which gives the values of parameters:
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:
By the help of the approximations (3.11), it is easy to verify that
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:
Applying the approximation (2.8) to the differential equation (4.1), we obtain
where we denote
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 , if the coefficient and/or the function contains the singular terms like , 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 .
We need the following approximations:
Using the approximations (4.4) in the scheme (4.2) and neglecting local truncation error, we get
The previous scheme in product form may be written as
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 ) as
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 ) in the homogeneous part of the scheme (4.8), and we get
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
On solving, we get , , and , provided that ,
Substituting the values of , , and in (4.10), we obtain
Let and , where
Substituting the values of and in (4.9), we obtain
Since , hence the method (4.8) is stable as long as
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
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:
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
The previous scheme is conditionally stable (see ). In order to obtain a superstable method, we simply follow the ideas given by Chawla  and the scheme (5.2) may be re-written as
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:
Using the transformation in (5.5), we obtain the transformed characteristic equation:
For stability of the scheme (5.4), we must have the following conditions:
It is easy to verify that the coefficient is as follows:
We can treat this case separately.
The coefficient is as follows:
For stability, it is required that the third coefficient is as follows:
Multiplying throughout of (5.11) by , we get
Thus the scheme is stable if
For or and , we have the characteristic equation:
whose roots are , .
Hence for , , , (which are independent of and ), the scheme (5.4) is superstable.
Further, the scheme (5.4) in product form may be written as
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: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 , , 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 for Example 6.1, and for other examples we have chosen the value of . 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 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 , this implies that all their successive tangential derivatives are known at ; that is, the values of , and so forth are known at .
Thus using the initial values and their successive tangential derivative values, from (6.2) we can obtain the value of , 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).
where and are real parameters. The exact solution is . The maximum absolute errors and CPU time (in seconds) are tabulated in Table 1 at for and .
Table 1: Example 6.1: the maximum absolute errors (using proposed -method) (with CPU time).
Example 6.2 (wave equation in cylindrical polar coordinates).
The aforementioned equation represents the two-space dimensional wave equation in cylindrical polar coordinates. The exact solution is . The maximum absolute errors and CPU time (in seconds) are tabulated in Table 2 at and .
Table 2: Example 6.2: the maximum absolute errors (with CPU time).
Example 6.3 (Van der Pol type nonlinear wave equation).
with exact solution . The maximum absolute errors and CPU time (in seconds) are tabulated in Table 3 at for , 2, and 3.
Table 3: Example 6.3: the maximum absolute errors (with CPU time).
Example 6.4 (dissipative nonlinear wave equation).
with exact solution . The maximum absolute errors and CPU time (in seconds) are tabulated in Table 4 at .
Table 4: Example 6.4: the maximum absolute errors (with CPU time).
Example 6.5 (nonlinear wave equation with variable coefficients).
with exact solution . The maximum absolute errors and CPU time (in seconds) are tabulated in Table 5 at for , 2 and 5.
Table 5: Example 6.5: The maximum absolute errors (with CPU time).
Example 6.6 (quasilinear hyperbolic equation).
with exact solution . The maximum absolute errors and CPU time (in seconds) are tabulated in Table 6 at for , 2, and 5. The order of convergence may be obtained by using the following formula:
where and 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, , for two linear problems, and , 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.
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.