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

we propose a new high-order approximation for the solution of two-space-dimensional quasilinear hyperbolic partial differential equation of the form utt A x, y, t, u uxx B x, y, t, u uyy g x, y, t, u, ux, uy, ut , 0 < x, y < 1, t > 0 subject to appropriate initial and Dirichlet boundary conditions , where k > 0 and h > 0 are mesh sizes in time and space directions, respectively. We use only five evaluations of the function g 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.


Introduction
We consider the following two-space dimensional quasilinear hyperbolic partial differential equation: u tt A x, y, t, u u xx B x, y, t, u u yy g x, y, t, u, u x , u y , u t , 0 < x, y < 1, t > 0 1.1 subject to the initial conditions u x, y, 0 φ x, y , u t x, y, 0 ψ x, y , 0 ≤ x, y ≤ 1, 1.2 Let h > 0 and k > 0 be the mesh spacing in the space and time directions, respectively.We replace the solution region Ω { x, y, t | 0 < x, y < 1, t > 0} by a set of grid points x l , y m , t j , where x l lh, l 0 1 N 1, y m mh, m 0 1 N 1, and t j jk, 0 < j < J, N and J being positive integers and N 1 h 1.Let p k/h > 0 be the mesh ratio parameter.The exact values of u x, y, t , A x, y, t , and B x, y, t at the grid point x l , y m , t j are denoted by U j l,m , A j l,m , and B j l,m , respectively.Similarly at the grid point x l , y m , t j , we denote A j x l,m ∂A j l,m /∂x, A j y l,m ∂A j l,m /∂y, A j xx l,m ∂ 2 A j l,m /∂x 2 , . . .and so forth.Let u j l,m be the approximate value of u x, y, t at the same grid point.
At the grid point x l , y m , t j , we consider the following approximations:

2.8
where . .and so forth and

3.2
Using Taylor expansion about the grid point x l , y m , t j , first we obtain

3.3
With the help of the approximations 2.2a -2.2e , and from 2.5a and 2.5b , we obtain

3.4
Now we define the following approximations: where a 1i and b 1i , i 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

3.6
where 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

3.9
In order to obtain the difference approximation of O k 2 k 2 h 2 h 4 the coefficient of k 2 h 2 in 3.9 must be zero, which gives the values of parameters:

3.10
Thus we obtain the difference method of O k 2 k 2 h 2 h 4 for the differential equation 4 for the solution of the quasilinear hyperbolic equation 1.1 .Whenever the coefficients are A A x, y, t, u and Advances in Mathematical Physics B B x, y, t, u , the difference scheme 2.8 needs to be modified.For this purpose, we use the following central differences:

3.11
where B x l , y m±1 , t j , U j l,m±1 .

3.12
By the help of the approximations 3.11 , it is easy to verify that . ., and so forth.

3.13
Thus substituting the central difference approximations 3.11 into 2.8 , we obtain the required numerical method of O k 2 k 2 h 2 h 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 20-24 .

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.

4.2
where we denote

4.3
Note that the linear scheme 4.2 is of O k 2 k 2 h 2 h 4 for the solution of differential equation 4.1 .However, the scheme 4.2 fails to compute at l 1, if the coefficient D x and/or the function f x, y, t contains the singular terms like 1/x, 1/x 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 x 0.
We need the following approximations:

4.4
Using the approximations 4.4 in the scheme 4.2 and neglecting local truncation error, we get where

4.6
The previous scheme in product form may be written as

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 f.

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 u j l,m A l B m ξ j e iβl e iγm where ξ e iφ such that |ξ| 1 in the homogeneous part of the scheme 4.8 , and we get where

4.10
where A and B 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
Substituting the values of A, A −1 , B and B −1 in 4.10 , we obtain

4.12
Let M 1 −2M 3 and M 2 −4M 4 , where Substituting the values of M 1 and M 2 in 4.9 , we obtain sin 2 φ 2 Since sin 2 φ/2 ≤ 1, 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 4.17 or, where u * l,m 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 y-and x-differences, respectively, which allows us to solve by sweeping first 4.17 in the y-and then 4.18 in the x-direction.It will be seen that these sweeps require only the solution of tri-diagonal systems.
Advances in Mathematical Physics 13

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 where

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 rewritten 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 u j l,m ξ j e iθl e iφm in the homogeneous part of the 5.4 , and we obtain the characteristic equation: where

5.6
Using the transformation ξ 1 z / 1 − z 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: and all θ, φ, except θ φ 0 or 2π.

5.9
We can treat this case separately.The coefficient is as follows: and all θ, φ.

5.10
For stability, it is required that the third coefficient is as follows:

5.11
Advances in Mathematical Physics 15 Multiplying throughout of 5.11 by 16η, we get

5.12
Thus the scheme is stable if
Further, the scheme 5.4 in product form may be written as FF.

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: where u * l,m is any intermediate value, and the intermediate boundary conditions required for the solution of u * l,m 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 < x, y < 1, t > 0 using a tridiagonal solver.The order of convergence may be obtained by using the following formula: where e h 1 and e h 2 are maximum absolute errors for two uniform mesh widths h 1 and h 2 , respectively.For computation of order of convergence of the proposed method, we have considered errors for last two values of h, that is, h 1 1/32, h 2 1/64 for two linear problems, and h 1 1/16, h 2 1/32 for all non-linear and quasilinear problems, and results are reported in Table 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 g as compared to seven evaluations of the function g discussed in 6, 7 , we have derived a new stable numerical method of O k 2 k 2 h 2 h 4 accuracy for the solution of quasilinear wave equation 1.1 .Ultimately, we use less algebra for computation, and for a fixed parameter value σ k/h 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.

Table 1 :
Example 6.1: the maximum absolute errors using proposed O k 2 k 2 h 2 h 4 -method with CPU time .The aforementioned equation represents the two-space dimensional wave equation in cylindrical polar coordinates.The exact solution is u cosh r cosh z sin t.The maximum absolute errors and CPU time in seconds are tabulated in Table2at t 1 and t 2. Van der Pol type nonlinear wave equation .u tt u xx u yy γ u 2 − 1 u t 2π 2 γ 2 e −2γt sin 2 πx sin 2 πy e −γt sin πx sin πy , −γt sin πx sin πy .The maximum absolute errors and CPU time in seconds are tabulated in Table 3 at t 2 for γ 1, 2, and 3.

Table 2 :
Example 6.2: the maximum absolute errors with CPU time .hProposed O k 2 k 2 h 2 h 4 -method O k 4 k 2 h 2 h 4 -method discussed in 6

Table 3 :
Example 6.3: the maximum absolute errors with CPU time .hProposedOk 2 k 2 h 2 h 4 -method O k 4 k 2 h 2 h 4 -method discussed in 6The maximum absolute errors and CPU time in seconds are tabulated in Table4at t 1. cosh x cosh y − sinh x y − 1 − x 2 − y 2 e −t cosh x cosh y, The maximum absolute errors and CPU time in seconds are tabulated in Table5at t 1 for γ 1, 2 and 5.
with exact solution u sin πx sin πy sin t. −t cosh x cosh y.

Table 4 :
Example 6.4: the maximum absolute errors with CPU time .hProposed O k 2 k 2 h 2 h 4 -method O k 4 k 2 h 2 h 4 -method discussed in 6

Table 5 :
Example 6.5:The maximum absolute errors with CPU time .hProposedOk 2 k 2 h 2 h 4 -method O k 4 k 2 h 2 h 4 -method discussed in 6The maximum absolute errors and CPU time in seconds are tabulated in Table6at t 1 for γ 1, 2, and 5.

Table 6 :
Example 6.6:The maximum absolute errors with CPU time .hProposed O k 2 k 2 h 2 h 4 -method O k 4 k 2 h 2 h 4 -method discussed in 7