# New Nonpolynomial Spline in Compression Method ofOk2+h4 for the Solution of 1D Wave Equation in  Polar Coordinates

Authors: Venu Gopal; R. K. Mohanty; Navnit Jha
Journal: Advances in Numerical Analysis
 (2013)

---

## Abstract

We propose a three-level implicit nine point compact finite difference formulation of order two in time and four in space direction, based on nonpolynomial spline in compression approximation inr-direction and finite difference approximation in t-direction for the numerical solution of one-dimensional wave equation in polar coordinates. We describe the mathematical formulation procedure in detail and also discussed the stability of the method. Numerical results are provided to justify the usefulness of the proposed method.


---

## Body

## 1. Introduction

We consider the one-dimensional wave equation in polar forms:(1)utt=urr+D(r)ur+E(r)u+f(r,t),0<r<1,t>0,
with the  following initial  conditions:
(2)u(r,0)=ϕ(r),ut(r,0)=φ(r),0≤r≤1,
and the following boundary conditions:
(3)u(0,t)=p0(t),u(1,t)=p1(t),t≥0,
where D(r)=γ/r and E(r)=-γ/r2.We assume that the conditions (2) and (3) are given with sufficient smoothness to maintain the order of accuracy in the numerical method under consideration.The study of wave equation in polar form is of keen interest in the fields like acoustics, electromagnetic, fluid dynamics, mathematical physics, and so forth.  Efforts are being made to develop efficient and high accuracy finite difference methods for such types of PDEs. During the last three decades, there has been much effort to develop stable numerical methods based on spline approximations for the solution of time-dependent partial differential equations. But so far in the literature, very limited spline methods are there for the wave equation in polar coordinates. In 1968-69, Bickley [1] and Fyfe [2] studied boundary value problems using cubic splines.  In 1973, Papamichael and Whiteman [3], and the next year, Fleck [4] and Raggett and Wilson [5] have used a cubic spline technique of lower order accuracy to solve one-dimensional heat conduction equation and wave equation, respectively. Then, Jain et al. [6–9] have derived cubic spline solution for the differential equations including fourth order cubic spline method for solving the nonlinear two point boundary value problems with significant first derivative terms. Recently, Kadalbajoo et al. [10, 11] and Khan et al. [12, 13] have studied parametric cubic spline technique for solving two point boundary value problems. In recent years, Rashidinia et al. [14], Ding and Zhang [15], and Mohanty et al. [16–21] have discussed spline and high order finite difference methods for the solution of hyperbolic equations. In this present paper, we follow the idea of Jain and Aziz [7] by using nonpolynomial spline in compression approximation to develop order four method in space direction for the wave equation in polar co-ordinates. We have shown that our method is in general of order four, but for the sake of computations, we have used the consistency of the first order continuity condition.In this paper, using nine grid points (see Figure1), we discuss a new three-level implicit non-polynomial spline finite difference method of accuracy two in time and four in space for the solution of one-dimensional wave equation in polar forms. In this method, we require only three evaluation of function G (which is defined in Section 2). In the next section, we discuss the non-polynomial spline in compression finite difference method. Difficulties were experienced in the past for the high order spline solution of wave equation in polar coordinates. The solution usually deteriorates in the vicinity of the singularity. In this section, we modify our technique in such a way that the solution retains its order and accuracy everywhere in the solution region. In this section, we also discussed the stability analysis of the proposed method. In Section 4, we discuss the higher order approximation at first time level in order to compute the proposed numerical method of same accuracy and compare the numerical results of proposed high accuracy non-polynomial spline in compression finite difference method with the corresponding second order accuracy non-polynomial spline in compression method. Concluding remarks are given in Section 5.Figure 1
Schematic representation of three-level implicit scheme.


## 2. The Numerical Method Based on Nonpolynomial Spline in Compression

The solution domain[0,1]×[t>0] is divided into (N+1)×J mesh with the spatial step size h=1/(N+1) in r-direction and the time step size k>0 in t-direction, respectively, where N and J are positive integers. The mesh ratio parameter is given by λ=(k/h)>0. Grid points are defined by (rl,tj)=(lh,jk), l=0,1,2,…,N+1, and j=0,1,2,…,J. The notations ulj and Ulj are used for the discrete approximation and the exact value of u(r,t) at the grid point (rl,tj), respectively.For the derivation of the non-polynomial spline in compression finite difference method for the solution of differential equation (1), we follow the ideas given by Jain and Aziz [7]. We use the non-polynomial spline in compression approximations in r-direction and second order finite difference approximation in t-direction.At the grid point(rl,tj), we may write the differential equation (1) as
(4)Uttlj-Urrlj=G(rl,tj,Ulj,Urlj)≡Glj(say),
where G(r,t,u,ur)=D(r)ur+E(r)u+f(r,t).LetSj(r) be the non-polynomial spline in compression interpolating function of the value ulj at the grid point (rl,tj) and is given by
(5)Sj(r)=alj+blj(r-rl)+cljsinτ(r-rl)+dljcosτ(r-rl),aaaaarl-1≤r≤rl;l=1,2,…,N+1;j=1,2,…,J,
where alj, blj, clj, and dlj are constants and ω is arbitrary parameter. Sj(r) is a class of C2[0,1] which interpolates u(r,t) at the grid point (rl,tj).The derivatives of non-polynomial spline in compression functionSj(r) are given by
(6)Sj′(r)=blj+τcljcosτ(r-rl)-τdljsinτ(r-rl),aaaaaaaaaaal=1,2,…,N+1;j=1,2,…,J,(7)Sj′′(r)=-τ2[cljsinτ(r-rl)+dljcosτ(r-rl)],aaaaaaaaaal=1,2,…,N+1;j=1,2,…,J,
where
(8)mlj=Sj′(rl)=Urlj,Mlj=Sj′′(rl)=Urrlj=Uttlj-D(rl)Urlj-E(rl)Ulj-f(rl,tj),l=0,1,2,…,N+1;j=1,2,…,J.
Substituting r=rl in (5) and r=rl-1  in (4) and (5), we obtain
(9)Sj′(rl)=h6[Ml-1j+2Mlj]+(ulj-ul-1jh)=mlj,(10)Sj′(rl-1)=-h[βMl-1j+αMlj]+(Ulj-Ul-1jh)=Url-1j=ml-1j,rl-1≤r≤rl.
By considering Sj(r) and Sj′(r) in rl≤r≤rl+1, we have
(11)Sj′(rl)=h6[Ml+1j+2Mlj]+(ul+1j-uljh)=mlj,(12)Sj′(rl+1)=h[βMl+1j+αMlj]+(Ul+1j-Uljh)=Url+1j=ml+1j,rl≤r≤rl+1.
Adding (9) and (11), we get
(13)mlj=Sj′(rl)=-h12[Ml+1j-Ml-1j]+(Ul+1j-Ul-1jh).
To derive expression for the coefficients of (5) in terms of  Ulj, Ul+1j, Mlj, and Ml+1j, we use
(14)Sj(rl)=Ulj,Sj(rl+1)=Ul+1j,Mlj=Sj′′(rl),Ml+1j=Sj′′(rl+1).
From algebraic manipulation, we get
(15)alj=Ulj+Mljω2,blj=Ul+1j-Uljh+Ml+1j-Mljωθ,clj=Mljcosθ-Ml+1jω2sinθ,dlj=-Mljω2,
where θ=ωh and l=0,1,2,…,N+1.Using the continuity of the first derivative at(rl,tj), that is, Sj′(rl-)=Sj′(rl+), we obtain the following relation for l=1,2,…,N-1:
(16)Ul+1j-2Ulj+Ul-1jh2=αMl+1j+2βMlj+αMl-1j,
where
(17)α=1θ2(θcosecθ-1),β=1θ2(1-θcotθ),θ=ωh.
When ω→0, that is, θ→0, then (α,β)→(1/6,1/3), and the relation (16) reduces to ordinary cubic spline relation:
(18)Ul+1j-2Ulj+Ul-1j=h26(Ml+1j+2Mlj+Ml-1j).
Note that (10), (12), (13), and (16) are important properties of the non-polynomial cubic spline in compression function Sj(r).We consider the following approximations:(19a)U¯ttlj=Ulj+1-2Ulj+Ulj-1k2=Uttlj+O(k2),(19b)U¯ttl+1j=Ul+1j+1-2Ul+1j+Ul+1j-1k2=Uttl+1j+O(k2+k2h),(19c)U¯ttl-1j=Ul-1j+1-2Ul-1j+Ul-1j-1k2=Uttl-1j+O(k2-k2h),(20a)
U
¯
r
l
j
=
U
l
+
1
j
-
U
l
-
1
j
2
h
=
U
r
l
j
+
h
2
6
U
r
r
r
l
j
+
O
(
h
4
)
,
(20b)
U
¯
r
l
+
1
j
=
3
U
l
+
1
j
-
4
U
l
j
+
U
l
-
1
j
2
h
=
U
r
l
+
1
j
-
h
2
3
U
r
r
r
l
j
+
O
(
h
3
)
,
(20c)
U
¯
r
l
-
1
j
=
-
3
U
l
-
1
j
+
4
U
l
j
-
U
l
+
1
j
2
h
=
U
r
l
-
1
j
-
h
2
3
U
r
r
r
l
j
-
O
(
h
3
)
,(21a)
G
¯
l
j
=
G
(
r
l
,
t
j
,
U
l
j
,
U
¯
r
l
j
)
,
(21b)
G
¯
l
+
1
j
=
G
(
r
l
+
1
,
t
j
,
U
l
+
1
j
,
U
¯
r
l
+
1
j
)
,
(21c)
G
¯
l
-
1
j
=
G
(
r
l
-
1
,
t
j
,
U
l
-
1
j
,
U
¯
r
l
-
1
j
)
.
Since the derivative values of Sj(r) defined by (10), (12) and (13) are not known at each grid point (rl,tj), we use the following approximations for the derivatives of Sj(r). Let(22a)M¯lj=(U¯ttlj-G¯lj),(22b)M¯l+1j=(U¯ttl+1j-G¯l+1j),(22c)M¯l-1j=(U¯ttl-1j-G¯l-1j),(23a)
m
^
l
j
=
U
l
+
1
j
-
U
l
-
1
j
2
h
-
α
h
2
[
M
¯
l
+
1
j
-
M
¯
l
-
1
j
]
,
(23b)
m
^
l
+
1
j
=
U
l
+
1
j
-
U
l
j
h
+
h
[
β
M
¯
l
+
1
j
+
α
M
¯
l
j
]
,
(23c)
m
^
l
-
1
j
=
U
l
j
-
U
l
-
1
j
h
-
h
[
β
M
¯
l
-
1
j
+
α
M
¯
l
j
]
.
Now we define the following approximations:(24a)G^lj=G(xl,tj,Ulj,m^lj),(24b)G^l+1j=G(xl+1,tj,Ul+1j,m^l+1j),(24c)G^l-1j=G(xl-1,tj,Ul-1j,m^l-1j),in which we use the non-polynomial spline in compression function Ulj=Sj(xl), approximation of its first order space derivative  defined by (23a)–(23c) in r-direction.With the help of the approximations (20a), from (21a), we obtain(25a)G¯lj=G(rl,tj,Ulj,Urlj+h26Urrrlj+O(h4))=Glj+h26Urrrlj(∂G∂Ur)lj+O(h4).
Similarly,
(25b)G¯l+1j=Gl+1j-h23Urrrlj(∂G∂Ur)lj+O(k2+h4),(25c)G¯l-1j=Gl-1j-h23Urrrlj(∂G∂Ur)lj+O(k2+h4).Now, with the help of the approximations (22a), (23a), and (25a), from (24a), we obtain(26a)G^lj=G(rl,tj,Ulj,mlj+O(k2+h4))=G(rl,tj,Ulj,mlj)+O(k2+h4)=Glj+O(k2+h4).
Similarly,
(26b)G^l+1j=Gl+1j+O(k2+h4),(26c)G^l-1j=Gl-1j+O(k2+h4).Then, at each grid point (rl,tj), a non-polynomial spline in compression finite difference method  with accuracy of O(k2+h4) for the solution of differential equation (1) may be written as
(27)6λ2[Ul+1j-2Ulj+Ul-1j]=k22[U¯ttl+1j+U¯ttl-1j+10U¯ttlj]-k22[G^l+1j+G^l-1j+10G^lj]+T^lj.
Using the approximations (19a)–(19c) and (26a)–(26c), from (27), we obtain the local truncation error:
(28)T^lj=6λ2[Ul+1j-2Ulj+Ul-1j]-k22[Uttl+1j+Uttl-1j+10Uttlj]+k22[Gl+1j+Gl-1j+10Glj]+O(k4+k2h4).
Now, substituting the values Glj=Uttlj-Urrlj and Gl±1j=Uttl±1j-Urrl±1j in (28) and then using Taylor’s expansion of Ul±1j, Uttl±1j and Urrl±1j at the grid point (rl,tj) in (28), that is, using the following values:
(29)Ul±1j=Ulj±hU10+h22U20±h36U30+h424U40±O(h5),Urrl±1j=Urrlj±hU30+h22U40±O(h3),…,
we obtain the local truncation error T^lj=O(k4+k2h4).Note that the initial and Dirichlet boundary conditions are given by (2) and (3), respectively. Incorporating the initial and boundary conditions, we can write the method (27) in a tri-diagonal matrix form. Since the differential equation (1) is linear, we can solve the linear system using the Gauss-elimination (tri-diagonal solver) method [22].


## 3. Stability Analysis

We can write the finite difference method based on non-polynomial spline in compression approximation (27) as follows by neglecting the LTE:
(30)6λ2[ul+1j-2ulj+ul-1j]=k22[u¯ttl+1j+u¯ttl-1j+10u¯ttlj]-k22[Dl+1u^rl+1j+Dl-1u^rl-1j+10Dlu^rlj]-k22[El+1u^l+1j+El-1u^l-1j+10Elu^lj]-k22[fl+1j+fl-1j+10flj],l=1(1)N,j=1,2,…,J,
where the approximations associated with (30) are defined in Section 2.Note that the scheme (30) is of O(k2+h4) accuracy for the solution of wave equation (1). Since r0=0, the scheme (30) fails to compute at l=1 due to zero division. In order to get a stable non-polynomial cubic spline in compression scheme of O(k2+h4) accuracy, we need the following approximations:(31a)Dl±1=Dl±hDrl+h2Drrl±O(h3),(31b)El±1=El±hErl+h2Errl±O(h3),(31c)fl±1j=flj±hfrlj+h22frrlj±O(h3),where
(32)Dl=D(rl),Drl=Dr(rl),Drrl=Drr(rl),El=E(rl),Erl=Er(rl),Errl=Err(rl),flj=f(rl,tj),frlj=fr(rl,tj),frrlj=frr(rl,tj),….
Now, with the help of the approximations defined in Section 2 and (31a)–(31c), neglecting high order terms, we can rewrite the scheme (30) in three-level operator compact implicit form:
(33)[R0+112(δr2+R1(2μrδr))]δt2ulj=λ2[R2δr2+R3(2μrδr)+2R4]ulj+∑f,l=1(1)N,j=1(1)J,
where
(34)R0=1+γ12l2,R1=12γl,R2=1+γ(γ-2)12l2,R3=R1+γ(6-γ)24l3,R4=-γ2l2+γ(6-γ)24l4,∑f=k212[(12+γl2)flj+γhlfrlj+h2frrjj],
and μrulj=(1/2)(ul+1/2j+ul-1/2j) and δrulj=(ul+1/2j-ul-1/2j)  are averaging and central difference operators with respect to r-direction, and so forth. This implies that (2μrδr)ulj=ul+1j-ul-1j, δr2ulj=ul+1j-2ulj+ul-1j, δt2ulj=ulj+1-2ulj+ulj-1, and so forth. The non-polynomial spline in compression finite difference scheme (33) has a local truncation error of O(k2+h4) and is free from the terms 1/(l±1), and hence, it can be solved for l=1(1)N, j=1(1)J in the region 0<r<1, t>0.For stability of the method (33), we follow the technique used by Mohanty [19]. We may rewrite (33) as
(35)[R0+112(R2δr2+R3(2μrδr))]δt2ulj=λ2[R2δr2+R3(2μrδr)+2R4]ulj+∑f.
The additional terms are of high orders and do not affect the accuracy of the scheme. The exact value Ulj=u(rl,tj) satisfies
(36)[R0+112(R2δr2+R3(2μrδr))]δt2Ulj=λ2[R2δr2+R3(2μrδr)+2R4]Ulj+∑f+O(k4+k2h4).
We assume that there exists an error εlj=Ulj-ulj at the grid point (xl,tj). Subtracting (35) from (36), we obtain the folowing error equation:
(37)[R0+112(R2δr2+R3(2μrδr))]δt2εlj=λ2[R2δr2+R3(2μrδr)+2R4]εlj+O(k4+k2h4).
For stability of the modified scheme (35), we assume that εlj=Aleiϕjeiθl (where ξ=eiϕ such that |ξ|=1) at the grid point (rl,tj), where ξ is in general complex, θ is an arbitrary real number, and A is a nonzero real parameter to be determined. Substituting εlj=Aleiϕjeiθl in the homogeneous part of the error equation (37), we obtain the amplification factor:(38)-4sin2(ϕ2)=(λ2[R2{(A+A-1)cosθ-2+i(A-A-1)sinθ}+R3{(A-A-1)cosθ+i(A+A-1)sinθ}+2R4])×(R0+112[R2{(A+A-1)cosθ-2+i(A-A-1)sinθ}+R3{(A-A-1)cosθ+i(A+A-1)sinθ}]112)-1.Since the left-hand side of (38) is a real quantity, the imaginary part of the right-hand side of (38) must be zero, from which we obtain
(39)R2(A-A-1)+R3(A+A-1)=0,
or
(40)A=R2-R3R2+R3,
where R2±R3>0. Substituting the values of A and A-1 in (38), we get(41)sin2(ϕ2)=λ2[R2+(R22-R32)(2sin2(θ/2)-1)-R4]2R0-(1/3)[R2+(R22-R32)(2sin2(θ/2)-1)].Since 0≤sin2(ϕ/2)≤1, maxsin2(θ/2)=1 and minsin2(θ/2)=0, it follows from (41) that the non-polynomial spline in compression finite difference scheme (35) is stable if
(42)0<λ2≤2R0-(1/3)[R2-R22-R32]R2-R4+R22-R32,
leading to |ξ|=1. It is easy to verify that as l→∞, 0<λ2≤1.


## 4. Numerical Illustrations

In this section, we have solved the problem (1)–(3) using the method described by (27) and compared our results with those obtained by the numerical method of O(k2+h2) accuracy based on non-polynomial spline in compression approximations and the method O(k4+h4) derived in [16] for the solution of 1D wave equations in polar form in different cases. The exact solution is provided.  The difference equation has been solved using a tri-diagonal solver. In order to demonstrate the fourth order convergence of the proposed method, throughout the computation we, have chosen the fixed value of the parameter σ=k/h2=3.2.  All computations were carried out using double precision arithmetic.Note that the proposed non-polynomial spline in compression finite difference method (27) for the 1D wave equations in polar form is a three-level scheme. The value of u at t=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 t=k. In this section, we discuss an explicit scheme of O(k2) for u at first time level, that is, at t=k in order to solve the differential equation (1) using the method (27).Since the values ofu and ut are known explicitly at t=0, this implies that all their successive tangential derivatives are known at t=0, that is, the values of u,ur,urr,…, ut,utr,…, and so forth are known at t=0.An approximation foru of O(k2) at t=k may be written as
(43)ul1=ul0+kutl0+k22(utt)l0+O(k3).
From (1), we have
(44)(utt)l0=[urr+G(r,t,u,ur)]l0.
Thus, using the initial values and their successive tangential derivative values, from (44), we can obtain the value of (utt)l0, and then ultimately, from (43), we can compute the value of u at first time level, that is, at t=k.Relation (16) is suitable for solving (1) provided it satisfies the consistency condition. That is, if ω is a root of the equation tan(ω/2)=ω/2. This equation has an infinite number of roots, the smallest positive nonzero root being given by ω=8.986818916⋯ [6].We solve (1) using the method (33) in the region bounded by 0<r<1, t>0 subject to the following initial conditions:
(45)u(r,0)=0,ut(r,0)=r2,0≤r≤1,
and the following boundary  conditions:
(46)u(0,t)=0,u(1,t)=sinht,t≥0.
The exact solution is given by u(r,t)=r2sinht. The maximum absolute errors (MAE) [23] are tabulated in Tables 1 and 2 at t=5.0 and for γ=1, γ=2. The exact and the numerical solutions are plotted in Figure 2(a) at γ=2, t=5.Table 1
The maximum absolute error atγ=1, t=5.
h
O
(
k
2
+
h
4
)-method
O
(
k
2
+
h
2
)-method
1/8
0.4559(−04)
0.2842(−03)
1/16
0.2850(−05)
0.7133(−04)
1/32
0.1781(−06)
0.1751(−04)
1/64
0.1113(−07)
0.4230(−05)Table 2
The maximum absolute error atγ=2, t=5.
h
O
(
k
2
+
h
4
)-method
O
(
k
2
+
h
2
)-method
1/8
0.3865(−04)
0.3736(−03)
1/16
0.2416(−05)
0.9379(−04)
1/32
0.1510(−06)
0.2267(−04)
1/64
0.9442(−08)
0.5137(−05)(a) Exact and numerical solution atγ=2 for t=5; (b) exact and numerical solution at γ=1 for t=1.
(a)
(b)We also solve (1) when E(r)=0 using the method (33) in the region bounded by 0<r<1, t>0 subject to the following initial conditions:
(47)u(r,0)=0,ut(r,0)=coshr,0≤r≤1,
and the following boundary  conditions:
(48)u(0,t)=sint,u(1,t)=12(e+e-1)sint,t≥0.
The exact solution is given by u(r,t)=coshrsint. The maximum absolute errors (MAE) [23] are tabulated in Tables 3 and 4 at t=2.0  for γ=1, γ=2. The exact and the numerical solutions are plotted in Figure 2(b) at γ=1, t=1.Table 3
The maximum absolute error atγ=1, t=2.
h
O
(
k
2
+
h
4
)-method
O
(
k
4
+
h
4
)-method [16]
1/8
0.3786(−04)
0.4848(−02)
1/16
0.2359(−05)
0.3113(−03)
1/32
0.1443(−06)
0.1677(−04)
1/64
0.8679(−08)
0.8892(−06)Table 4
The maximum absolute error atγ=2, t=2.
h
O
(
k
2
+
h
4
)-method
O
(
k
4
+
h
4
)-method [16]
1/8
0.1889(−04)
0.5488(−02)
1/16
0.1137(−05)
0.3862(−03)
1/32
0.6943(−06)
0.2400(−04)
1/64
0.4169(−07)
0.1515(−05)


## 5. Final Remarks

Available numerical methods based on non-polynomial spline in compression approximations for the numerical solution of the 1D wave equations in polar form are ofO(k2+h2) accuracy only and require 9-grid points. In this paper, using the same number of grid points and three evaluations of the function G (which is defined in Section 2), we have derived a new stable non-polynomial spline in compression finite difference method  of O(k2+h4) accuracy for the solution of the 1D wave equation (1). For a fixed parameter σ=k/h2, the proposed method behaves like a fourth order method, which is exhibited from the computed results. The proposed numerical method for the wave equation in polar coordinates is conditionally stable.



---

*Source: 470480-2013-09-30.xml*