About this Journal Submit a Manuscript Table of Contents
ISRN Computational Mathematics
VolumeΒ 2012Β (2012), Article IDΒ 302923, 6 pages
http://dx.doi.org/10.5402/2012/302923
Research Article

Cubic Spline Method for 1D Wave Equation in Polar Coordinates

1Department of Mathematics, Faculty of Mathematical Sciences, University of Delhi, Delhi 110 007, India
2Department of Mathematics, Deenbandhu Chhotu Ram University of Science & Technology, Murthal 131039, India

Received 10 September 2011; Accepted 28 September 2011

Academic Editors: P. B.Β Vasconcelos and J. G.Β Zhou

Copyright Β© 2012 R. K. Mohanty et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Using nonpolynomial cubic spline approximation in space and finite difference in time direction, we discuss three-level implicit difference scheme of 𝑂(π‘˜2+β„Ž4) for the numerical solution of 1D wave equations in polar coordinates, where π‘˜>0 and β„Ž>0 are grid sizes in time and space coordinates, respectively. The proposed method is applicable to problems with singularity. Stability theory of the proposed method is discussed, and numerical examples are given in support of the theoretical results.

1. Introduction

We consider the 1D linear singular hyperbolic equation of the formπœ•2π‘’πœ•π‘‘2=πœ•2π‘’πœ•π‘Ÿ2+π›Όπ‘Ÿπœ•π‘’πœ•π‘Ÿ+𝑓(π‘Ÿ,𝑑),0<π‘Ÿ<1,𝑑>0,𝛼=1,2(1) subject to the initial conditions𝑒(π‘Ÿ,0)=πœ™(π‘Ÿ),𝑒𝑑(π‘Ÿ,0)=πœ“(π‘Ÿ),0β‰€π‘Ÿβ‰€1(2) and with boundary conditions at π‘Ÿ=0andπ‘Ÿ=1of the form𝑒(0,𝑑)=𝑔0(𝑑),𝑒(1,𝑑)=𝑔1(𝑑),𝑑β‰₯0,(3) where 𝑒=𝑒(π‘Ÿ,𝑑),𝑑 and π‘Ÿ are time and distance variable, respectively. For 𝛼=1 and 2, the equation above represent, 1D wave equation in cylindrical and spherical polar coordinates, respectively. We shall assume that the initial and boundary conditions are given with sufficient smoothness to maintain the order of accuracy of the difference scheme and spline functions under consideration.

In this paper, we are interested to discuss a new approximation based on cubic spline polynomial for the solution to singular hyperbolic equation (1). During last three decades, several numerical schemes for the solution of two-point boundary value problems and partial differential equations have been developed by many researchers. First Bickley [1] and Fyfe [2] have discussed the second-order accurate spline method for the solution of linear two-point boundary value problems. Jain and Aziz [3] have derived fourth-order cubic spline method for solving the nonlinear two-point boundary value problems with significant first derivative terms. Further, Khan and Aziz [4] have studied parametric cubic spline method to the solution of a system of second-order boundary value problem. In 1974, Raggett and Wilson [5] have used a cubic spline technique of lower order accuracy to solve the wave equation. Later, Fleck Jr [6] has proposed a cubic spline method for solving the wave equation of linear optics. In recent years, Mohanty [7], Gao and Chi [8], Rashidinia et al. [9, 10], Mohebbi and Dehghan [11], Ding and Zhang [12, 13], and Liu and Liu [14] have derived various numerical methods for solution of nonsingular 2D hyperbolic equations. A difference method of accuracy two in time and four in space for the solution of differential equation (1) has been studied by Mohanty [15] using finite difference method. A new discretization method of order four for the numerical solution of one-space dimensional second-order nonlinear hyperbolic equations has been studied by Mohanty et al. [16, 17]. Recently, Mohanty and Dahiya [18] have derived highly accurate cubic spline method for the solution of parabolic equations. In this paper, using nine grid points, we discuss a new three-level implicit cubic spline finite difference method of accuracy two in time and four in space for the solution of differential equation (1). To the authors’ knowledge, no cubic spline method of accuracy two in time and four in space for the solution of (1) has been discussed in the literature so far. In next section, we discuss the derivation of the proposed cubic spline method. It has been experienced in the past that the cubic spline solutions for the wave equation in polar coordinates usually deteriorate in the vicinity of the singularity. We overcome this difficulty by modifying the method in such a way that the solutions retain its order and accuracy everywhere in the vicinity of the singularity. In Section 3, we discuss the linear stability analysis of the proposed cubic spline method. In Section 4, we compare the computed results with the one obtained by the method discussed in [17]. Final remarks are given in Section 5.

2. The Method Based on Cubic Spline Approximation

The solution domain [0,1]Γ—[𝑑>0] is divided into (𝑁+1)×𝐽 mesh with the spatial step size β„Ž=1/(𝑁+1)in π‘Ÿ-direction and the time step size π‘˜>0 in 𝑑-direction, respectively, where 𝑁 and 𝐽 are positive integers. The mesh ratio parameter is given by πœ†=π‘˜/β„Ž>0. Grid points are defined by (π‘Ÿπ‘™,𝑑𝑗)=(π‘™β„Ž,π‘—π‘˜),𝑙=0,1,2,…,𝐽. The notations 𝑒𝑗𝑙andπ‘ˆπ‘—π‘™ are used for the discrete approximation and the exact value of 𝑒(π‘Ÿ,𝑑) at the grid point (π‘Ÿπ‘™,𝑑𝑗), respectively.

At the grid point (π‘Ÿπ‘™,𝑑𝑗), we denoteπ‘ˆπ‘Žπ‘=πœ•π‘Ž+π‘π‘ˆπœ•π‘Ÿπ‘Žπ‘™πœ•π‘‘π‘π‘—,π‘“π‘Žπ‘=πœ•π‘Ž+π‘π‘“πœ•π‘Ÿπ‘Žπ‘™πœ•π‘‘π‘π‘—.(4)

Let 𝑆𝑗(π‘Ÿ) be the cubic spline interpolating polynomial of the function 𝑒(π‘Ÿ,𝑑𝑗) between the grid points (π‘Ÿπ‘™βˆ’1,𝑑𝑗) and (π‘Ÿπ‘™,𝑑𝑗) and be given byπ‘†π‘—ξ€·π‘Ÿ(π‘Ÿ)=π‘™ξ€Έβˆ’π‘Ÿ3𝑀6β„Žπ‘—π‘™βˆ’1+ξ€·π‘Ÿβˆ’π‘Ÿπ‘™βˆ’1ξ€Έ3𝑀6β„Žπ‘—π‘™+ξ‚΅π‘ˆπ‘—π‘™βˆ’1βˆ’β„Ž26π‘€π‘—π‘™βˆ’1ξ‚Άξ‚€π‘Ÿπ‘™βˆ’π‘Ÿβ„Žξ‚+ξ‚΅π‘ˆπ‘—π‘™βˆ’β„Ž26π‘€π‘—π‘™ξ‚Άξ‚€π‘Ÿβˆ’π‘Ÿπ‘™βˆ’1β„Žξ‚,π‘Ÿπ‘™βˆ’1β‰€π‘Ÿβ‰€π‘Ÿπ‘™,𝑙=1,2,…,𝑁+1,𝑗=1,2,…,𝐽,(5) which satisfies at 𝑗th level the following properties:(i)𝑆𝑗(π‘Ÿ) coincides with a polynomial of degree three on each [π‘Ÿπ‘™βˆ’1,π‘Ÿπ‘™], 𝑙=1,2,…,𝑁+1, 𝑗=1,2,…𝐽,(ii)𝑆𝑗(π‘Ÿ)∈𝐢2[0,1], and(iii)𝑆𝑗(π‘Ÿπ‘™)=π‘ˆπ‘—π‘™, 𝑙=0,1,2,…,𝑁+1, 𝑗=1,2,…,𝐽.

The derivatives of cubic spline function 𝑆𝑗(π‘Ÿ) are given by π‘†ξ…žπ‘—βˆ’ξ€·π‘Ÿ(π‘Ÿ)=π‘™ξ€Έβˆ’π‘Ÿ2𝑀2β„Žπ‘—π‘™βˆ’1+ξ€·rβˆ’π‘Ÿπ‘™βˆ’1ξ€Έ2𝑀2β„Žπ‘—π‘™+π‘ˆπ‘—π‘™βˆ’π‘ˆπ‘—π‘™βˆ’1β„Žβˆ’β„Ž6ξ€Ίπ‘€π‘—π‘™βˆ’π‘€π‘—π‘™βˆ’1ξ€»,π‘†π‘—ξ…žξ…žξ€·π‘Ÿ(π‘Ÿ)=π‘™ξ€Έβˆ’π‘Ÿβ„Žπ‘€π‘—π‘™βˆ’1+ξ€·π‘Ÿβˆ’π‘Ÿπ‘™βˆ’1ξ€Έβ„Žπ‘€π‘—π‘™,(6) where𝑀𝑗𝑙=π‘†π‘—ξ…žξ…žξ€·π‘Ÿπ‘™ξ€Έ=π‘ˆπ‘—π‘Ÿπ‘Ÿπ‘™=π‘ˆπ‘—π‘‘π‘‘π‘™βˆ’π›Όπ‘Ÿπ‘™π‘ˆπ‘—π‘Ÿπ‘™βˆ’π‘“π‘—π‘™,π‘šπ‘™=0,1,2,…,𝑁+1,𝑗=1,2,…,𝐽,(7)𝑗𝑙=π‘†ξ…žπ‘—ξ€·π‘Ÿπ‘™ξ€Έ=π‘ˆπ‘—π‘Ÿπ‘™=π‘ˆπ‘—π‘™βˆ’π‘ˆπ‘—π‘™βˆ’1β„Ž+β„Ž6ξ€Ίπ‘€π‘—π‘™βˆ’1+2𝑀𝑗𝑙,π‘Ÿπ‘™βˆ’1β‰€π‘Ÿβ‰€π‘Ÿπ‘™,(8) and replacing β„Ž by β€œβˆ’β„Ž,” we getπ‘šπ‘—π‘™=π‘†ξ…žπ‘—ξ€·π‘Ÿπ‘™ξ€Έ=π‘ˆπ‘—π‘Ÿπ‘™=π‘ˆπ‘—π‘™+1βˆ’π‘ˆπ‘—π‘™β„Žβˆ’β„Ž6𝑀𝑗𝑙+1+2𝑀𝑗𝑙,π‘Ÿπ‘™βˆ’1β‰€π‘Ÿβ‰€π‘Ÿπ‘™.(9) Combining (8) and (9), we obtain π‘šπ‘—π‘™=π‘†ξ…žπ‘—ξ€·π‘Ÿπ‘™ξ€Έ=π‘ˆπ‘—π‘Ÿπ‘™=π‘ˆπ‘—π‘™+1βˆ’π‘ˆπ‘—π‘™βˆ’1βˆ’β„Ž2β„Žξ€Ίπ‘€12𝑗𝑙+1βˆ’π‘€π‘—π‘™βˆ’1ξ€».(10) Further, from (8), we haveπ‘šπ‘—π‘™+1=π‘†ξ…žπ‘—ξ€·π‘Ÿπ‘™+1ξ€Έ=π‘ˆπ‘—π‘Ÿπ‘™+1=π‘ˆπ‘—π‘™+1βˆ’π‘ˆπ‘—π‘™β„Ž+β„Ž6𝑀𝑗𝑙+2𝑀𝑗𝑙+1ξ€»,(11) and, from (9), we haveπ‘šπ‘—π‘™βˆ’1=π‘†ξ…žπ‘—ξ€·π‘Ÿπ‘™βˆ’1ξ€Έ=π‘ˆπ‘—π‘Ÿπ‘™βˆ’1=π‘ˆπ‘—π‘™βˆ’π‘ˆπ‘—π‘™βˆ’1β„Žβˆ’β„Ž6𝑀𝑗𝑙+2π‘€π‘—π‘™βˆ’1ξ€».(12) We consider the following approximations:π‘ˆπ‘—π‘‘π‘™=ξ‚€π‘ˆπ‘™π‘—+1βˆ’π‘ˆπ‘™π‘—βˆ’12π‘˜=π‘ˆπ‘‘π‘—π‘™ξ€·π‘˜+𝑂2ξ€Έ,π‘ˆπ‘—π‘‘π‘™+1=ξ‚€π‘ˆπ‘—+1𝑙+1βˆ’π‘ˆπ‘—βˆ’1𝑙+12π‘˜=π‘ˆπ‘—π‘‘π‘™+1ξ€·π‘˜+𝑂2+π‘˜2β„Žξ€Έ,π‘ˆπ‘—π‘‘π‘™βˆ’1=ξ‚€π‘ˆπ‘—+1π‘™βˆ’1βˆ’π‘ˆπ‘—βˆ’1π‘™βˆ’12π‘˜=π‘ˆπ‘—π‘‘π‘™βˆ’1ξ€·π‘˜+𝑂2βˆ’π‘˜2β„Žξ€Έ,π‘ˆπ‘—π‘‘π‘‘π‘™=ξ‚€π‘ˆπ‘™π‘—+1βˆ’2π‘ˆπ‘—π‘™+π‘ˆπ‘™π‘—βˆ’1ξ‚π‘˜2=π‘ˆπ‘—π‘‘π‘‘π‘™ξ€·π‘˜+𝑂2ξ€Έ,π‘ˆπ‘—π‘‘π‘‘π‘™+1=ξ‚€π‘ˆπ‘—+1𝑙+1βˆ’2π‘ˆπ‘—π‘™+1+π‘ˆπ‘—βˆ’1𝑙+1ξ‚π‘˜2=π‘ˆπ‘—π‘‘π‘‘π‘™+1ξ€·π‘˜+𝑂2+π‘˜2β„Žξ€Έ,π‘ˆπ‘—π‘‘π‘‘π‘™βˆ’1=ξ‚€π‘ˆπ‘—+1π‘™βˆ’1βˆ’2π‘ˆπ‘—π‘™βˆ’1+π‘ˆπ‘—βˆ’1π‘™βˆ’1ξ‚π‘˜2=π‘ˆπ‘—π‘‘π‘‘π‘™βˆ’1ξ€·π‘˜+𝑂2βˆ’π‘˜2β„Žξ€Έ,π‘ˆπ‘—π‘Ÿπ‘™=ξ€·π‘ˆπ‘—π‘™+1βˆ’π‘ˆπ‘—π‘™βˆ’1ξ€Έ2β„Ž=π‘ˆπ‘—π‘Ÿπ‘™+β„Ž26π‘ˆ30ξ€·π‘˜+𝑂2+β„Ž4ξ€Έ,π‘ˆπ‘—π‘Ÿπ‘™+1=ξ€·3π‘ˆπ‘—π‘™+1βˆ’4π‘ˆπ‘—π‘™+π‘ˆπ‘—π‘™βˆ’1ξ€Έ2β„Ž=π‘ˆπ‘—π‘Ÿπ‘™+1βˆ’β„Ž23π‘ˆ30ξ€·π‘˜+𝑂2+π‘˜2β„Žξ€Έ,π‘ˆπ‘—π‘Ÿπ‘™βˆ’1=ξ€·βˆ’3π‘ˆπ‘—π‘™βˆ’1+4π‘ˆπ‘—π‘™βˆ’π‘ˆπ‘—π‘™+1ξ€Έ2β„Ž=π‘ˆπ‘—π‘Ÿπ‘™βˆ’1βˆ’β„Ž23π‘ˆ30ξ€·π‘˜+𝑂2βˆ’π‘˜2β„Žξ€Έ.(13) Since the derivative values of 𝑆𝑗(π‘Ÿ) defined by (7), (10), (11), and (12) are not known at each grid point (π‘Ÿπ‘™,𝑑𝑗), we use the following approximations for the derivatives of 𝑆𝑗(π‘Ÿ).

Let 𝑀𝑗𝑙=π‘ˆπ‘—π‘‘π‘‘π‘™βˆ’π›Όπ‘Ÿπ‘™π‘ˆπ‘—π‘Ÿπ‘™βˆ’π‘“π‘—π‘™,𝑀𝑗𝑙+1=π‘ˆπ‘—π‘‘π‘‘π‘™+1βˆ’π›Όπ‘Ÿπ‘™+1π‘ˆπ‘—π‘Ÿπ‘™+1βˆ’π‘“π‘—π‘™+1,π‘€π‘—π‘™βˆ’1=π‘ˆπ‘—π‘‘π‘‘π‘™βˆ’1βˆ’π›Όπ‘Ÿπ‘™βˆ’1π‘ˆπ‘—π‘Ÿπ‘™βˆ’1βˆ’π‘“π‘—π‘™βˆ’1,ξπ‘šπ‘—π‘™=π‘ˆπ‘—π‘™+1βˆ’π‘ˆπ‘—π‘™βˆ’1βˆ’β„Ž2β„Žξ‚ƒ12𝑀𝑗𝑙+1βˆ’π‘€π‘—π‘™βˆ’1ξ‚„,ξπ‘šπ‘—π‘™+1=π‘ˆπ‘—π‘™+1βˆ’π‘ˆπ‘—π‘™β„Ž+β„Ž6𝑀𝑗𝑙+2𝑀𝑗𝑙+1ξ‚„,ξπ‘šπ‘—π‘™βˆ’1=π‘ˆπ‘—π‘™βˆ’π‘ˆπ‘—π‘™βˆ’1β„Žβˆ’β„Ž6𝑀𝑗𝑙+2π‘€π‘—π‘™βˆ’1ξ‚„,(14) where ξπ‘šπ‘—π‘™ and 𝑀𝑗𝑙 are approximations to π‘’π‘—π‘Ÿπ‘™ and π‘’π‘—π‘Ÿπ‘Ÿπ‘™, respectively.

We may rewrite the differential equation (1) asπœ•2π‘’πœ•π‘Ÿ2=πœ•2π‘’πœ•π‘‘2βˆ’π›Όπ‘Ÿπœ•π‘’πœ•π‘Ÿβˆ’π‘“(π‘Ÿ,𝑑)β‰‘πœ™(π‘Ÿ,𝑑)(say).(15)

Let us denote πœ™π‘—π‘™=πœ™(π‘Ÿπ‘™,𝑑𝑗).

A fourth-order method (see [3]) based on cubic spline approximations for the differential equation (15) may be written asπ‘ˆπ‘—π‘™+1βˆ’2π‘ˆπ‘—π‘™+π‘ˆπ‘—π‘™βˆ’1=β„Ž2ξ‚ƒξπœ™12𝑗𝑙+1+ξπœ™π‘—π‘™βˆ’1ξπœ™+10𝑗𝑙,(16) whereξπœ™π‘—π‘™+1=π‘ˆπ‘—π‘‘π‘‘π‘™+1βˆ’π›Όπ‘Ÿπ‘™+1ξπ‘šπ‘—π‘™+1βˆ’π‘“π‘—π‘™+1,ξπœ™π‘—π‘™βˆ’1=π‘ˆπ‘—π‘‘π‘‘π‘™βˆ’1βˆ’π›Όπ‘Ÿπ‘™βˆ’1ξπ‘šπ‘—π‘™βˆ’1βˆ’π‘“π‘—π‘™βˆ’1,ξπœ™π‘—π‘™=π‘ˆπ‘—π‘‘π‘‘π‘™βˆ’π›Όπ‘Ÿπ‘™ξπ‘šπ‘—π‘™βˆ’π‘“π‘—π‘™.(17)

Multiplying throughout by 6πœ†2 and by the help of the approximations (17), from (16), we obtain a cubic spline finite difference method with accuracy of 𝑂(π‘˜2+β„Ž4) for the solution of differential equation (1) as6πœ†2ξ€Ίπ‘ˆπ‘—π‘™+1βˆ’2π‘ˆπ‘—π‘™+π‘ˆπ‘—π‘™βˆ’1ξ€»=π‘˜22ξ‚ƒπ‘ˆπ‘—π‘‘π‘‘π‘™+1+π‘ˆπ‘—π‘‘π‘‘π‘™βˆ’1+10π‘ˆπ‘—π‘‘π‘‘π‘™ξ‚„βˆ’π‘˜22ξ‚Έπ›Όπ‘Ÿπ‘™+1ξπ‘šπ‘—π‘™+1+π›Όπ‘Ÿπ‘™βˆ’1ξπ‘šπ‘—π‘™βˆ’1𝛼+10π‘Ÿπ‘™ξπ‘šπ‘—π‘™ξ‚Ήβˆ’π‘˜22𝑓𝑗𝑙+1+π‘“π‘—π‘™βˆ’1+10𝑓𝑗𝑙.(18)

Note that the method (18) is of 𝑂(π‘˜2+β„Ž4) for the numerical solution of (1). However, the method fails to compute at 𝑙=1. We modify the method in such a manner that it retains its order and accuracy in the vicinity of the singularity.

We use the following approximations:1π‘Ÿπ‘™+1=1π‘Ÿπ‘™βˆ’β„Žπ‘Ÿ2𝑙+β„Ž2π‘Ÿ3π‘™ξ€·β„Žβˆ’π‘‚3ξ€Έ,1π‘Ÿπ‘™βˆ’1=1π‘Ÿπ‘™+β„Žπ‘Ÿ2𝑙+β„Ž2π‘Ÿ3π‘™ξ€·β„Ž+𝑂3ξ€Έ,𝑓𝑗𝑙+1=𝑓𝑗𝑙+β„Žπ‘“10+β„Ž22𝑓20ξ€·β„Ž+𝑂3ξ€Έ,π‘“π‘—π‘™βˆ’1=π‘“π‘—π‘™βˆ’β„Žπ‘“10+β„Ž22𝑓20ξ€·β„Žβˆ’π‘‚3ξ€Έ,π›Όπ‘Ÿπ‘™+1π‘ˆπ‘—π‘Ÿπ‘™+1+𝑓𝑗𝑙+1=1𝛼2β„Žπ‘Ÿπ‘™+1ξ€·2πœ‡π‘Ÿπ›Ώπ‘Ÿ+2𝛿2π‘Ÿξ€Έπ‘ˆπ‘—π‘™+𝑓𝑗𝑙+1,π›Όπ‘Ÿπ‘™βˆ’1π‘ˆπ‘—π‘Ÿπ‘™βˆ’1+π‘“π‘—π‘™βˆ’1=1𝛼2β„Žπ‘Ÿπ‘™βˆ’1ξ€·2πœ‡π‘Ÿπ›Ώπ‘Ÿβˆ’2𝛿2π‘Ÿξ€Έπ‘ˆπ‘—π‘™+π‘“π‘—π‘™βˆ’1,π›Όπ‘Ÿπ‘™π‘ˆπ‘—π‘Ÿπ‘™+𝑓𝑗𝑙=1𝛼2β„Žπ‘Ÿπ‘™ξ€·2πœ‡π‘Ÿπ›Ώπ‘Ÿξ€Έπ‘ˆπ‘—π‘™+𝑓𝑗𝑙.(19) where πœ‡π‘Ÿπ‘’π‘™=(1/2)(𝑒1+(1/2)+𝑒1βˆ’(1/2)) and π›Ώπ‘Ÿπ‘’π‘™=(𝑒1+(1/2)βˆ’π‘’1βˆ’(1/2)). Substituting the approximations (19) into (18) and neglecting higher order terms, we getξ€·12+𝛿2π‘Ÿξ€Έπ›Ώ2π‘‘π‘ˆπ‘—π‘™βˆ’12πœ†2𝛿2π‘Ÿπ‘ˆπ‘—π‘™=πœ†2β„Ž2𝛼12π‘Ÿπ‘™+2β„Ž2π›Όπ‘Ÿ3𝑙ξƒͺξ€·2πœ‡π‘Ÿπ›Ώπ‘Ÿξ€Έπ‘ˆπ‘—π‘™+πœ†2β„Ž2𝛼(π›Όβˆ’1)π‘Ÿ2𝑙ξƒͺ𝛿2π‘Ÿπ‘ˆπ‘—π‘™βˆ’β„Ž2π›Όπ‘Ÿ2𝑙𝛿2π‘‘π‘ˆπ‘—π‘™βˆ’β„Žπ›Ό2π‘Ÿπ‘™ξ€·π›Ώ2𝑑2πœ‡π‘Ÿπ›Ώπ‘Ÿξ€Έπ‘ˆπ‘—π‘™+π‘˜2ξ‚Έ12𝑓00+β„Ž2𝑓20+π›Όπ‘Ÿπ‘™ξ€·π‘“10βˆ’π‘“00ξ€Έ.ξ‚Άξ‚Ή(20)

Note that the cubic spline method (20) is of 𝑂(π‘˜2+β„Ž4) for the numerical solution of differential equation (1), which is free from the terms 1/π‘Ÿπ‘™Β±1 hence can be computed for 𝑙=1(1)𝑁,𝑗=0,1,2,….

3. Stability Analysis

Now we discuss the stability analysis for the scheme (20). For stability, we consider the homogeneous part of the scheme (20), which can be written as 𝑅0+1𝛿122π‘Ÿ+𝑅1ξ€·2πœ‡π‘Ÿπ›Ώπ‘Ÿξ‚„π›Ώξ€Έξ€Έ2π‘‘π‘ˆπ‘—π‘™=πœ†2𝑅2𝛿2π‘Ÿ+𝑅3ξ€·2πœ‡π‘Ÿπ›Ώπ‘Ÿπ‘ˆξ€Έξ€»π‘—π‘™,(21) where 𝑅0β„Ž=1+2𝛼12π‘Ÿ2,𝑅1=β„Žπ›Ό2π‘Ÿ,𝑅2β„Ž=1+2ξ‚Έ12𝛼(π›Όβˆ’1)π‘Ÿ2ξ‚Ή,𝑅3=𝑅1+β„Ž3ξ‚€242π›Όπ‘Ÿ3.(22)

It is difficult to find the stability region for the scheme (21). In order to obtain a valid stability region, we may modify the scheme (21) (see [19]) as𝑅0+1𝑅122𝛿2π‘Ÿ+𝑅3ξ€·2πœ‡π‘Ÿπ›Ώπ‘Ÿξ‚„π›Ώξ€Έξ€Έ2π‘‘π‘ˆπ‘—π‘™=πœ†2𝑅2𝛿2π‘Ÿ+𝑅3ξ€·2πœ‡π‘Ÿπ›Ώπ‘Ÿπ‘ˆξ€Έξ€»π‘—π‘™.(23)

The additional terms added in (23) are of high order and do not affect the accuracy of the scheme.

Put π‘ˆπ‘—π‘™=π΄π‘™πœ‰π‘—π‘’π‘–π›½π‘™=π΄π‘™π‘’π‘–πœ™π‘—π‘’π‘–π›½π‘™ (where πœ‰=π‘’π‘–πœ™ such that |πœ‰|=1) in the modified scheme (23), we getπœ‰βˆ’2+πœ‰βˆ’1=βˆ’4sin2πœ™2=πœ†2𝑅2𝐴+π΄βˆ’1ξ€Έξ€·cosπ›½βˆ’2+π‘–π΄βˆ’π΄βˆ’1ξ€Έξ€Ύsin𝛽+𝑅3ξ€½ξ€·π΄βˆ’π΄βˆ’1ξ€Έξ€·cos𝛽+𝑖𝐴+π΄βˆ’1ξ€Έsin𝛽𝑅0+1𝑅122𝐴+π΄βˆ’1ξ€Έξ€·cosπ›½βˆ’2+π‘–π΄βˆ’π΄βˆ’1ξ€Έξ€Ύsin𝛽+𝑅3ξ€½ξ€·π΄βˆ’π΄βˆ’1ξ€Έξ€·cos𝛽+𝑖𝐴+π΄βˆ’1ξ€Έ,sin𝛽(24)where 𝐴 is a nonzero real parameter to be determined. Left-hand side of (24) is a real quantity. Thus, the imaginary part of right-hand side of (24) must be zero.

Thus, we obtain𝑅2ξ€·π΄βˆ’π΄βˆ’1ξ€Έ+𝑅3𝐴+π΄βˆ’1ξ€ΈβŸΉξ€·π‘…=02+𝑅3ξ€Έξ€·π‘…π΄βˆ’2βˆ’π‘…3ξ€Έπ΄βˆ’1ξƒŽ=0⟹𝐴=𝑅2βˆ’π‘…3𝑅2+𝑅3,π΄βˆ’1=ξƒŽπ‘…2+𝑅3𝑅2βˆ’π‘…3(25) provided𝑅2±𝑅3>0,thatis,βˆ’π‘…2<𝑅3<𝑅2.

Hence, 𝐴+π΄βˆ’1=2𝑅2/𝑅22βˆ’π‘…23andπ΄βˆ’π΄βˆ’1=βˆ’2𝑅3/𝑅22βˆ’π‘…23.

Substituting the values of 𝐴,π΄βˆ’1,𝐴+π΄βˆ’1,andπ΄βˆ’π΄βˆ’1in (24) and (25), we obtainβˆ’4sin2πœ™2=πœ†2ξ‚Έ2𝑅22βˆ’π‘…23cosπ›½βˆ’2𝑅2𝑅0+(1/6)𝑅22βˆ’π‘…23cosπ›½βˆ’π‘…2ξ‚Ή=2πœ†2𝑅22βˆ’π‘…23cosπ›½βˆ’π‘…2𝑅0+(1/6)𝑅22βˆ’π‘…23cosπ›½βˆ’π‘…2ξ‚ΉβŸΉsin2πœ™2=πœ†2𝑅2βˆ’ξ”π‘…22βˆ’π‘…23ξ‚Ήcos𝛽2𝑅0+(1/3)𝑅22βˆ’π‘…23cosπ›½βˆ’π‘…2ξ‚Ή=πœ†2𝑅2+𝑅22βˆ’π‘…23ξ€·ξ€·2sin2𝛽/2βˆ’12𝑅0ξ‚Έπ‘…βˆ’(1/3)2+𝑅22βˆ’π‘…23ξ€·ξ€·2sin2≑𝑁𝛽/2βˆ’1𝐷.(26)

Since 0≀sin2πœ™/2≀1, the method (23) is stable as long as 𝑁≀𝐷, which is true if 𝑁≀Max(𝑁)≀Min(𝐷)≀𝐷.

Hence, the required stability condition isξ‚΅πœ†Max2𝑅2+𝑅22βˆ’π‘…23ξ‚΅2sin2𝛽2ξ‚΅βˆ’1≀Min2𝑅0βˆ’13𝑅2+𝑅22βˆ’π‘…23ξ‚΅2sin2𝛽2βˆ’1ξ‚Άξ‚Ήξ‚ΆβŸΉπœ†2𝑅2+𝑅22βˆ’π‘…23≀2𝑅0βˆ’13𝑅2βˆ’ξ”π‘…22βˆ’π‘…23ξ‚Ή0<πœ†2≀2𝑅0ξ‚Έπ‘…βˆ’(1/3)2βˆ’ξ”π‘…22βˆ’π‘…23𝑅2+𝑅22βˆ’π‘…23,(27) which is the required stability interval for the scheme (23).

4. Numerical Results

A difference method of 𝑂(π‘˜2+β„Ž2), for the differential equation (1) may be written as𝑒𝑗𝑑𝑑𝑙=π‘’π‘—π‘Ÿπ‘Ÿπ‘™+π›Όπ‘Ÿπ‘’π‘—π‘Ÿπ‘™ξ€·π‘Ÿ+𝑓𝑙,𝑑𝑗,𝑙=1(1)𝑁,𝑗=0,1,2,….(28)

Note that the proposed cubic spline method (23) and the difference method (28) for second-order hyperbolic equation (1) are three-level schemes. The value of u at 𝑑=0 is known from the initial condition. To start any computation, it is necessary to know the numerical value of u of required accuracy at 𝑑=π‘˜. In this section, we discuss an explicit scheme of 𝑂(π‘˜2) for u at first time level, that is, at 𝑑=π‘˜, in order to solve the differential equation (1) using the proposed method (23) and second-order method (28).

Since the values of 𝑒 and 𝑒𝑑 are known explicitly at 𝑑=0, this implies 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ξ€Έ.(29)

From (1), we have𝑒𝑑𝑑0𝑙=ξ‚ƒπ‘’π‘Ÿπ‘Ÿ+π›Όπ‘Ÿπ‘’π‘Ÿξ‚„+𝑓(π‘Ÿ,𝑑)0𝑙.(30)

Thus, using the initial values and their successive tangential derivative values, from (30), we can obtain the value of (𝑒𝑑𝑑)0l, and, then ultimately, from (29), we can compute the value of 𝑒 at first time level, that is, at 𝑑=π‘˜.

We solve the differential equation (1) in the region 0<π‘Ÿ<1,𝑑>0, whose exact solution is given by 𝑒(π‘Ÿ,𝑑)=cosh(π‘Ÿ)β‹…sin(𝑑). The maximum absolute errors at 𝑑=5.0 are tabulated in Table 1 for various values of πœ†=π‘˜/β„Ž=0.8 and in Table 2 for 𝛾=π‘˜/β„Ž2=4.0. We have compared the numerical results of the proposed method with the results obtained by using the method discussed in [17] in terms of accuracy. All computations were performed using double precision arithmetic.

tab1
Table 1: πœ†=π‘˜/β„Ž=0.8, 𝑑=5.0.
tab2
Table 2: 𝛾=π‘˜/β„Ž2=4.0, 𝑑=5.0.

A relation between the exact solution 𝑒𝑒π‘₯π‘Žπ‘π‘‘ and the approximate numerical solution 𝑒(β„Ž) as given in the following equation:𝑒exact=𝑒(β„Ž)+π΄β„Žπ‘+β‹―higher-orderterms,(31) where β„Ž is the measure of the mesh discretization, 𝐴 is a constant, and 𝑝 is the order (rate) of convergence. If the meshes to be considered are sufficiently refined, the higher-order terms can be neglected. Then, the maximum absolute errors πΈβ„Ž can be approximated as πΈβ„Ž||𝑒=Maxexact||βˆ’π‘’(β„Ž)β‰…π΄β„Žπ‘.(32)

Taking the logarithm of both sides of (32), we obtain𝐸logβ„Žξ€Έ=log(𝐴)+𝑝log(β„Ž).(33)

For two different refined mesh spacing β„Ž1 and β„Ž2, we have the following two relations𝐸logβ„Ž1ξ€Έξ€·β„Ž=log(𝐴)+𝑝log1ξ€Έ,(34a)𝐸logβ„Ž2ξ€Έξ€·β„Ž=log(𝐴)+𝑝log2ξ€Έ.(34b)

Subtracting (34b) from (34a), we obtain the order (rate) of convergence 𝐸𝑝=logβ„Ž1ξ€Έξ€·πΈβˆ’logβ„Ž2ξ€Έξ€·β„Žlog1ξ€Έξ€·β„Žβˆ’log2ξ€Έ,(35) where πΈβ„Ž1 and πΈβ„Ž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 β„Ž1=1/40 and β„Ž2=1/80 in Table 2, and we found the order of convergence of the proposed method for 𝛼=1is 4.08 and for 𝛼=2 is 4.01.

5. Final Remarks

Available numerical methods based on spline approximations for the numerical solution of 1D wave equation in polar coordinates are of 𝑂(π‘˜2+β„Ž2) accurate, which require nine grid points. In this paper, using the same number of grid points, we have discussed a new stable three-level implicit cubic spline finite difference method of 𝑂(π‘˜2+β„Ž4) accuracy for the solution of wave equation in polar coordinates. For a fixed parameter 𝛾=π‘˜/β„Ž2, the proposed method behaves like a fourth-order method, which is exhibited from the computed results.

Acknowledgment

The authors thank the anonymous reviewers for their constructive suggestions, which substantially improved the standard of the paper.

References

  1. W. G. Bickley, β€œPiecewise cubic interpolation and two-point boundary value problems,” The Computer Journal, vol. 11, pp. 206–208, 1968.
  2. D. J. Fyfe, β€œThe use of cubic splines in the solution of two-point boundary value problems,” The Computer Journal, vol. 12, pp. 188–192, 1969. View at Publisher Β· View at Google Scholar
  3. M. K. Jain and T. Aziz, β€œCubic spline solution of two-point boundary value problems with significant first derivatives,” Computer Methods in Applied Mechanics and Engineering, vol. 39, no. 1, pp. 83–91, 1983. View at Publisher Β· View at Google Scholar
  4. A. Khan and T. Aziz, β€œParametric cubic spline approach to the solution of a system of second-order boundary-value problems,” Journal of Optimization Theory and Applications, vol. 118, no. 1, pp. 45–54, 2003. View at Publisher Β· View at Google Scholar Β· View at MathSciNet
  5. G. F. Raggett and P. D. Wilson, β€œA fully implicit finite difference approximation to the one-dimensional wave equation using a cubic spline technique,” Journal of the Institute of Mathematics and its Applications, vol. 14, pp. 75–77, 1974. View at Publisher Β· View at Google Scholar
  6. J. A. Fleck Jr., β€œA cubic spline method for solving the wave equation of nonlinear optics,” Journal of Computational Physics, vol. 16, pp. 324–341, 1974. View at Publisher Β· View at Google Scholar
  7. R. K. Mohanty, β€œAn unconditionally stable difference scheme for the one-space-dimensional linear hyperbolic equation,” Applied Mathematics Letters, vol. 17, no. 1, pp. 101–105, 2004. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH Β· View at MathSciNet
  8. F. Gao and C. Chi, β€œUnconditionally stable difference schemes for a one-space-dimensional linear hyperbolic equation,” Applied Mathematics and Computation, vol. 187, no. 2, pp. 1272–1276, 2007. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH Β· View at MathSciNet
  9. J. Rashidinia, R. Jalilian, and V. Kazemi, β€œSpline methods for the solutions of hyperbolic equations,” Applied Mathematics and Computation, vol. 190, no. 1, pp. 882–886, 2007. View at Publisher Β· View at Google Scholar Β· View at MathSciNet
  10. J. Rashidinia, R. Mohammadi, and R. Jalilian, β€œSpline methods for the solution of hyperbolic equation with variable coefficients,” Numerical Methods for Partial Differential Equations, vol. 23, no. 6, pp. 1411–1419, 2007. View at Publisher Β· View at Google Scholar Β· View at MathSciNet
  11. A. Mohebbi and M. Dehghan, β€œHigh order compact solution of the one-space-dimensional linear hyperbolic equation,” Numerical Methods for Partial Differential Equations, vol. 24, no. 5, pp. 1222–1235, 2008. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH
  12. H. Ding and Y. Zhang, β€œParameters spline methods for the solution of hyperbolic equations,” Applied Mathematics and Computation, vol. 204, no. 2, pp. 938–941, 2008. View at Publisher Β· View at Google Scholar
  13. H. Ding and Y. Zhang, β€œA new unconditionally stable compact difference scheme of O(τ2+h4) for the 1D linear hyperbolic equation,” Applied Mathematics and Computation, vol. 207, no. 1, pp. 236–241, 2009. View at Publisher Β· View at Google Scholar
  14. H.-W. Liu and L.-B. Liu, β€œAn unconditionally stable spline difference scheme of O(k2+h4) for solving the second-order 1D linear hyperbolic equation,” Mathematical and Computer Modelling, vol. 49, no. 9-10, pp. 1985–1993, 2009. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH
  15. R. K. Mohanty, β€œNew unconditionally stable difference schemes for the solution of multi-dimensional telegraphic equations,” International Journal of Computer Mathematics, vol. 86, no. 12, pp. 2061–2071, 2009. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH
  16. R. K. Mohanty, M. K. Jain, and K. George, β€œOn the use of high order difference methods for the system of one space second order nonlinear hyperbolic equations with variable coefficients,” Journal of Computational and Applied Mathematics, vol. 72, no. 2, pp. 421–431, 1996. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH Β· View at MathSciNet
  17. R. K. Mohanty and U. Arora, β€œA new discretization method of order four for the numerical solution of one-space dimensional second-order quasi-linear hyperbolic equation,” International Journal of Mathematical Education in Science and Technology, vol. 33, no. 6, pp. 829–838, 2002. View at Publisher Β· View at Google Scholar Β· View at MathSciNet
  18. R. K. Mohanty and V. Dahiya, β€œAn O(k2+kh2+h4) accurate two-level implicit cubic spline method for one space dimensional Quasi-linear parabolic equations,” American Journal of Computational Mathematics, vol. 1, pp. 11–17, 2011.
  19. R. K. Mohanty, β€œStability interval for explicit difference schemes for multi-dimensional second-order hyperbolic equations with significant first-order space derivative terms,” Applied Mathematics and Computation, vol. 190, no. 2, pp. 1683–1690, 2007. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH Β· View at MathSciNet