Cubic Spline Method for 1D Wave Equation in Polar Coordinates

Abstract

Using nonpolynomial cubic spline approximation in space and finite difference in time direction, we discuss three-level implicit difference scheme of 𝑂(𝑘²+ℎ⁴) 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𝜕²𝑢𝜕𝑡²=𝜕²𝑢𝜕𝑟²+𝛼𝑟𝜕𝜕𝑢𝑟+𝑓(𝑟,𝑡),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,𝑡)=𝑔₀(𝑡),𝑢(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𝑆_𝑗(𝑟)=𝑟^𝑙−𝑟³6ℎ𝑀^𝑗𝑙−1+𝑟−𝑟^𝑙−1³6ℎ𝑀^𝑗𝑙+𝑈^𝑗𝑙−1−ℎ²6𝑀^𝑗𝑙−1𝑟^𝑙−𝑟+𝑈_𝑗𝑙ℎ−ℎ²𝑟_𝑙−16𝑀^𝑗𝑙𝑟−ℎ𝑟^𝑙−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)𝑆_𝑗(𝑟)∈𝐶²[0,1], and(iii)𝑆_𝑗(𝑟_𝑙)=𝑈^𝑗𝑙, 𝑙=0,1,2,…,𝑁+1, 𝑗=1,2,…,𝐽.

The derivatives of cubic spline function 𝑆_𝑗(𝑟) are given by 𝑆^𝑗(𝑟)=−𝑟^𝑙−𝑟²2ℎ𝑀^𝑗𝑙−1+r−𝑟^𝑙−1²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−𝑈^𝑗𝑙−12ℎ−1ℎ2𝑀^𝑗𝑙+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𝑘=𝑈^𝑡𝑗𝑙+𝑂𝑘²,𝑈^{𝑗𝑡𝑙+1}=𝑈^𝑗+1𝑙+1−𝑈^{𝑗−1𝑙+1}2𝑘=𝑈^{𝑗𝑡𝑙+1}+𝑂𝑘²+𝑘²ℎ,𝑈^{𝑗𝑡𝑙−1}=𝑈^{𝑗+1𝑙−1}−𝑈^{𝑗−1𝑙−1}2𝑘=𝑈^{𝑗𝑡𝑙−1}+𝑂𝑘²−𝑘²ℎ,𝑈^{𝑗𝑡𝑡𝑙}=𝑈^𝑙𝑗+1−2𝑈^𝑗𝑙+𝑈^𝑙𝑗−1𝑘²=𝑈^{𝑗𝑡𝑡𝑙}+𝑂𝑘²,𝑈^{𝑗𝑡𝑡𝑙+1}=𝑈^𝑗+1𝑙+1−2𝑈^𝑗𝑙+1+𝑈^{𝑗−1𝑙+1}𝑘²=𝑈^{𝑗𝑡𝑡𝑙+1}+𝑂𝑘²+𝑘²ℎ,𝑈^{𝑗𝑡𝑡𝑙−1}=𝑈^{𝑗+1𝑙−1}−2𝑈^𝑗𝑙−1+𝑈^{𝑗−1𝑙−1}𝑘²=𝑈^{𝑗𝑡𝑡𝑙−1}+𝑂𝑘²−𝑘²ℎ,𝑈^𝑗𝑟𝑙=𝑈^𝑗𝑙+1−𝑈^𝑗𝑙−12ℎ=𝑈^𝑗𝑟𝑙+ℎ²6𝑈³⁰+𝑂𝑘²+ℎ⁴,𝑈^{𝑗𝑟𝑙+1}=3𝑈^𝑗𝑙+1−4𝑈^𝑗𝑙+𝑈^𝑗𝑙−1=𝑈_{𝑗𝑟𝑙+1}2ℎ−ℎ²3𝑈³⁰+𝑂𝑘²+𝑘²ℎ,𝑈^{𝑗𝑟𝑙−1}=−3𝑈^𝑗𝑙−1+4𝑈^𝑗𝑙−𝑈^𝑗𝑙+1=𝑈_{𝑗𝑟𝑙−1}2ℎ−ℎ²3𝑈³⁰+𝑂𝑘²−𝑘²ℎ.(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−𝑈^𝑗𝑙−12ℎ−1ℎ2𝑚_𝑗𝑙+1𝑀^𝑗𝑙+1−𝑀^𝑗𝑙−1,=𝑈^𝑗𝑙+1−𝑈^𝑗𝑙ℎ+ℎ6𝑚_𝑗𝑙−1𝑀^𝑗𝑙+2𝑀^𝑗𝑙+1,=𝑈^𝑗𝑙−𝑈^𝑗𝑙−1ℎ−ℎ6𝑀^𝑗𝑙+2𝑀^𝑗𝑙−1,(14) where 𝑚^𝑗𝑙 and 𝑀^𝑗𝑙 are approximations to 𝑢^𝑗𝑟𝑙 and 𝑢^{𝑗𝑟𝑟𝑙}, respectively.

We may rewrite the differential equation (1) as𝜕²𝑢𝜕𝑟²=𝜕²𝑢𝜕𝑡²−𝛼𝑟𝜕𝜕𝑢𝑟−𝑓(𝑟,𝑡)≡𝜙(𝑟,𝑡)(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=ℎ²12𝜙^𝑗𝑙+1+𝜙^𝑗𝑙−1+10𝜙^𝑗𝑙,(16) where𝜙^𝑗𝑙+1=𝑈^{𝑗𝑡𝑡𝑙+1}−𝛼𝑟_𝑙+1𝜙_𝑗𝑙−1𝑚^𝑗𝑙+1−𝑓^𝑗𝑙+1,=𝑈^{𝑗𝑡𝑡𝑙−1}−𝛼𝑟_𝑙−1𝜙_𝑗𝑙𝑚^𝑗𝑙−1−𝑓^𝑗𝑙−1,=𝑈^{𝑗𝑡𝑡𝑙}−𝛼𝑟_𝑙𝑚^𝑗𝑙−𝑓^𝑗𝑙.(17)

Multiplying throughout by 6𝜆² and by the help of the approximations (17), from (16), we obtain a cubic spline finite difference method with accuracy of 𝑂(𝑘²+ℎ⁴) for the solution of differential equation (1) as6𝜆²=𝑈𝑘^𝑗𝑙+1₂−2𝑈^𝑗𝑙+𝑈^𝑗𝑙−12−𝑘₂𝑈^{𝑗𝑡𝑡𝑙+1}+𝑈^{𝑗𝑡𝑡𝑙−1}+10𝑈^{𝑗𝑡𝑡𝑙}2𝑟𝛼_𝑙+1𝑚^𝑗𝑙+1+𝛼𝑟_𝑙−1𝑚^𝑗𝑙−1+10𝛼𝑟_𝑙𝑚^𝑗𝑙−𝑘₂2𝑓^𝑗𝑙+1+𝑓^𝑗𝑙−1+10𝑓^𝑗𝑙.(18)

Note that the method (18) is of 𝑂(𝑘²+ℎ⁴) 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𝑙+ℎ²𝑟^3𝑙1−𝑂ℎ³,𝑟_𝑙−1=1𝑟_𝑙+ℎ𝑟^2𝑙+ℎ²𝑟^3𝑙𝑓_𝑗𝑙+1+𝑂ℎ³,=𝑓^𝑗𝑙+ℎ𝑓₁₀+ℎ²2𝑓²⁰𝑓_𝑗𝑙−1+𝑂ℎ³,=𝑓^𝑗𝑙−ℎ𝑓₁₀+ℎ²2𝑓²⁰𝛼−𝑂ℎ³,𝑟_𝑙+1𝑈^{𝑗𝑟𝑙+1}+𝑓^𝑗𝑙+1=12ℎ𝑟𝛼_𝑙+12𝜇_𝑟𝛼𝛿_𝑟+2𝛿^2𝑟𝑈^𝑗𝑙+𝑓^𝑗𝑙+1,𝑟_𝑙−1𝑈^{𝑗𝑟𝑙−1}+𝑓^𝑗𝑙−1=12ℎ𝑟𝛼_𝑙−12𝜇_𝑟𝛼𝛿_𝑟−2𝛿^2𝑟𝑈^𝑗𝑙+𝑓^𝑗𝑙−1,𝑟_𝑙𝑈^𝑗𝑟𝑙+𝑓^𝑗𝑙=12ℎ𝑟𝛼_𝑙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12𝑟𝛼_𝑙+2ℎ²𝛼𝑟^3𝑙2𝜇^𝑟𝛿_𝑟𝑈^𝑗𝑙+𝜆²ℎ²𝛼(𝛼𝑟^2𝑙−1)𝛿^2𝑟𝑈^𝑗𝑙−ℎ²𝛼𝑟^2𝑙𝛿^2𝑡𝑈^𝑗𝑙−ℎ𝛼2𝑟_𝑙𝛿^2𝑡2𝜇_𝑟𝛿_𝑟𝑈^𝑗𝑙+𝑘²12𝑓⁰⁰+ℎ²𝑓²⁰+𝛼𝑟_𝑙𝑓₁₀−𝑓₀₀.(20)

Note that the cubic spline method (20) is of 𝑂(𝑘²+ℎ⁴) 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 𝑅⁰+112𝛿^2𝑟+𝑅₁2𝜇_𝑟𝛿_𝑟𝛿^2𝑡𝑈^𝑗𝑙=𝜆²𝑅₂𝛿^2𝑟+𝑅₃2𝜇_𝑟𝛿_𝑟𝑈^𝑗𝑙,(21) where 𝑅₀=1+ℎ²𝛼12𝑟²,𝑅₁=ℎ𝛼2𝑟,𝑅²=1+ℎ²𝑅₃12𝛼(𝛼𝑟²−1),=𝑅₁+ℎ³242𝑟𝛼³.(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𝑅⁰+112𝑅²𝛿^2𝑟+𝑅₃2𝜇_𝑟𝛿_𝑟𝛿^2𝑡𝑈^𝑗𝑙=𝜆²𝑅₂𝛿^2𝑟+𝑅₃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+𝜉⁻¹=−4sin²𝜙=𝜆₂2𝑅₂𝐴+𝐴⁻¹cos𝛽−2+𝑖𝐴−𝐴⁻¹sin𝛽+𝑅₃𝐴−𝐴⁻¹cos𝛽+𝑖𝐴+𝐴⁻¹sin𝛽𝑅₀+112𝑅²𝐴+𝐴⁻¹cos𝛽−2+𝑖𝐴−𝐴⁻¹sin𝛽+𝑅₃𝐴−𝐴⁻¹cos𝛽+𝑖𝐴+𝐴⁻¹sin𝛽(24,)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𝑅⟹₂𝐴−𝑅𝐴⁻¹+𝑅₃𝐴+𝐴⁻¹=0₂⟹𝐴=+𝑅₃𝐴−𝑅₂−𝑅₃𝐴⁻¹=0𝑅₂−𝑅₃𝑅₂+𝑅₃,𝐴⁻¹=𝑅₂+𝑅₃𝑅₂−𝑅₃(25) provided𝑅₂±𝑅₃>0,thatis,−𝑅₂<𝑅₃<𝑅₂.

Hence, 𝐴+𝐴⁻¹=2𝑅₂/𝑅²²−𝑅²³and𝐴−𝐴⁻¹=−2𝑅₃/𝑅²²−𝑅²³.

Substituting the values of 𝐴,𝐴⁻¹,𝐴+𝐴⁻¹,and𝐴−𝐴⁻¹in (24) and (25), we obtain−4sin²𝜙2=𝜆²2𝑅²²−𝑅²³cos𝛽−2𝑅₂𝑅₀+(1/6)𝑅²²−𝑅²³cos𝛽−𝑅₂=2𝜆₂𝑅²²−𝑅²³cos𝛽−𝑅₂𝑅₀+(1/6)𝑅²²−𝑅²³cos𝛽−𝑅₂⟹sin₂𝜙2=𝜆²𝑅²−𝑅²²−𝑅²³cos𝛽2𝑅₀+(1/3)𝑅²²−𝑅²³cos𝛽−𝑅₂=𝜆₂𝑅²+𝑅²²−𝑅²³2sin²𝛽/2−1≡2𝑁𝑅₀−(1/3)𝑅₂+𝑅²²−𝑅²³2sin²𝛽/2−1𝐷.(26)

Since 0≤sin²𝜙/2≤1, the method (23) is stable as long as 𝑁≤𝐷, which is true if 𝑁≤Max(𝑁)≤Min(𝐷)≤𝐷.

Hence, the required stability condition isMax𝜆²𝑅²+𝑅²²−𝑅²³2sin²𝛽≤Min2𝑅₀2−1−13𝑅²+𝑅²²−𝑅²³2sin²𝛽⟹𝜆₂2−1𝑅²+𝑅²²−𝑅²³≤2𝑅⁰−13𝑅²−𝑅²²−𝑅²³0<𝜆₂≤2𝑅⁰−(1/3)𝑅₂−𝑅²²−𝑅²³𝑅₂+𝑅²²−𝑅²³,(27) which is the required stability interval for the scheme (23).

4. Numerical Results

A difference method of 𝑂(𝑘²+ℎ²), 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 𝑂(𝑘²) 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 𝑂(𝑘²) at 𝑡=𝑘 may be written as𝑢^1𝑙=𝑢^0𝑙+𝑘𝑢^0𝑡𝑙+𝑘²2𝑢^𝑡𝑡^0𝑙+𝑂𝑘³.(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 (𝑢_𝑡𝑡)⁰_l, 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 𝛾=𝑘/ℎ²=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.

Table 1 

𝜆=𝑘/ℎ=0.8, 𝑡=5.0.

Table 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 𝐸_ℎ=Max||𝑢_exact−𝑢(ℎ)||≅𝐴ℎ^𝑝.(32)

Taking the logarithm of both sides of (32), we obtainlog𝐸_ℎ=log(𝐴)+𝑝log(ℎ).(33)

For two different refined mesh spacing ℎ₁ and ℎ₂, we have the following two relationslog𝐸_ℎ1=log(𝐴)+𝑝logℎ₁,(34a)log𝐸_ℎ2=log(𝐴)+𝑝logℎ₂.(34b)

Subtracting (34b) from (34a), we obtain the order (rate) of convergence 𝑝=log𝐸^ℎ1−log𝐸_ℎ2logℎ₁−logℎ₂,(35) where 𝐸_ℎ1 and 𝐸_ℎ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 ℎ₁=1/40 and ℎ₂=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 𝑂(𝑘²+ℎ²) 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 𝑂(𝑘²+ℎ⁴) accuracy for the solution of wave equation in polar coordinates. For a fixed parameter 𝛾=𝑘/ℎ², 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.