# Cubic Spline Iterative Method for Poisson’s Equation in Cylindrical Polar Coordinates

Authors: R. K. Mohanty; Rajive Kumar; Vijay Dahiya
Journal: ISRN Mathematical Physics
 (2012)

---

## Abstract

Using nonpolynomial cubic spline approximation inx- and finite difference in y-direction, we discuss a numerical approximation of O(k2+h4) for the solutions of diffusion-convection equation, where k>0 and h>0 are grid sizes in y- and x-coordinates, respectively. We also extend our technique to polar coordinate system and obtain high-order numerical scheme for Poisson’s equation in cylindrical polar coordinates. Iterative method of the proposed method is discussed, and numerical examples are given in support of the theoretical results.


---

## Body

## 1. Introduction

We consider the two-dimensional elliptic equation of the form(1.1)∂2u∂x2+∂2u∂y2=D(x)∂u∂x+g(x,y),(x,y)∈Ω.
The Dirichlet boundary conditions are given by(1.2)u(x,y)=g(x,y),(x,y)∈Γ,
where Ω={(x,y)∣0<x,y<1}is the solution domain and Γ is its boundary. For g(x,y)=0 and D(x)=β, the above equation represents diffusion-convention equation. For D(x)=-1/x, the above equations represent the Poisson’s equation with singular coefficients in rectangular coordinates. Similarly, for D(x)=-1/x and replacing the variables x,y by r,z, we obtain a Poisson’s equation in cylindrical polar coordinates. We will assume that the 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 in discussing a new approximation based on cubic spline polynomial for the solution of elliptic equation (1.1). In many practical problems, coefficients of the second derivatives term are small compared to the coefficients of the first derivatives term. These problems are called singular perturbation problems. Singularly perturbed elliptic boundary value problems are mathematically models of diffusion-convections process or related physical phenomenon. The diffusion term is the term involving the second-order derivative, and convective term is that involving the first-order derivative. During last three decades, several numerical schemes for the solution of elliptic partial differential equations have been developed by many researchers. First-Lynch and Rice [1] have discussed high-accuracy finite difference approximations to the solutions of elliptic partial differential equations. Boisvert [2] has discussed a class of high-order accurate discretization for the elliptic boundary value problems. Yavneh [3] has reported the analysis of fourth-order compact scheme for convection diffusion equation. A fourth-order difference method for elliptic equations with non-linear first derivative terms has been discussed by Jain et al. [4, 5]. In 1997, Mohanty [6] has derived order h4 difference method for a class of 2D elliptic boundary value problems with singular coefficients. A new discretization method of order four for the numerical solution of 2D non-linear elliptic partial differential equations has been studied by Mohanty et al. [7–9]. The use of cubic spline polynomial and its approximation plays an important role for the formation of stable numerical methods. In the past, many authors (see [10–12]) have studied and analysed the use of cubic spline approximations in the solution of linear two-point boundary value problems. In 1983, Jain and Aziz [13] have developed a new method based on cubic spline approximations for the solution of two-point nonlinear boundary value problems. Later, Al-Said [14, 15] has discussed cubic spline methods for solving the system of second-order boundary value problems. Khan [16] has introduced parametric cubic spline approach for solving second order ordinary differential equations. Mohanty et al. [17, 18] have also reported a fourth-order accurate cubic spline alternating group explicit method for nonlinear singular two-point boundary value problems. Recently, Rashidinia et al. [19] have proposed a new cubic spline technique for two-point boundary value problems. Most recently, Mohanty and Dahiya [20] and Mohanty et al. [21] have used cubic spline polynomials and developed high order stable numerical methods for the solution of one space dimensional parabolic and hyperbolic partial differential equations. To the authors knowledge, no high-order method using cubic spline polynomial for the solution of 2D elliptic differential equations (1.1) has been discussed in the literature so far. In this paper, using nine-point compact cell (see Figure 1), we discuss a new compact cubic spline finite difference method of accuracy two in y- and four in x-coordinates for the solution of elliptic differential equation (1.1). In the 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 Poisson’s 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 solution retains its order and accuracy everywhere in the vicinity of the singularity. In Section 3, we discuss an iterative method. In Section 4, we compare the computed results with the results obtained by the method discussed in [8]. Concluding remarks are given in Section 5.Figure 1
9-point computational network.


## 2. The Approximation Based on Cubic Spline Polynomial

We consider our region of interest, a rectangular domainΩ=[0,1]×[0,1]. We choose grid spacing h>0 and k>0 in the directions x- and y- respectively, so that the mesh points (xl,ym) denoted by (l,m) are defined by xl=lh and ym=mk;l=0,1,…,N+1;m=0,1,…,M+1, where N and M are positive integers such that (N+1)h=1 and (M+1)k=1. The mesh ratio parameter is denoted by λ=(k/h). The notations ul,mandUl,m are used for the discrete approximation and the exact value of u(x,y) at the grid point (xl,ym), respectively. Similarly, we write D(xl)=Dl and D(xl±1)=Dl±1.At the grid point(xl,ym), we denote
(2.1)Wab=∂a+bW∂xla∂ymb,W=U,Dandg.LetSm(x) be the cubic spline interpolating polynomial of the function u(x,ym) between the grid point (xl-1,ym) and (xl,ym) and is given by(2.2)Sm(x)=(xl-x)36hMl-1,m+(x-xl-1)36hMl,m+(Ul-1,m-h26Ml-1,m)(xl-xh)+(Ul,m-h26Ml,m)⋅(x-xl-1h),xl-1≤x≤xl;l=1,2,…,N+1;m=0,1,2,…,M+1,
which satisfies at mth-line parallel to x-axis the following properties:(i)
Sm(x) coincides with a polynomial of degree three on each [xl-1,xl],l=1,2,…,N+1andm=1,2,…M,(ii)
Sm(x)∈C2[0,1],(iii)
Sm(xl)=Ul,m,l=0,1,2,…,N+1andm=1,2,…,M.The derivatives of cubic spline functionSm(x) are given by(2.3)Sm′(x)=-(xl-x)22hMl-1,m+(x-xl-1)22hMl,m+Ul,m-Ul-1,mh-h6[Ml,m-Ml-1,m],Sm′′(x)=(xl-x)hMl-1,m+(x-xl-1)hMl,m,
where(2.4)Ml,m=Sm′′(xl)=Uxxl,m=-Uyyl,m+DlUxl,m+gl,m,l=0,1,2,…,N+1,j=1,2,…,J,(2.5)ml,m=Sm′(xl)=Uxl,m=Ul,m-Ul-1,mh+h6[Ml-1,m+2Ml,m],xl-1≤x≤xl,
and replacing h by “-h”, we get(2.6)ml,m=Sm′(xl)=Uxl,m=Ul+1,m-Ul,mh-h6[Ml+1,m+2Ml,m],xl≤x≤xl+1,
Combining (2.5) and (2.6), we obtain(2.7)ml,m=Sm′(xl)=Uxl,m=Ul+1,m-Ul-1,m2h-h12[Ml+1,m-Ml-1,m],
Further, from (2.5), we have(2.8)ml+1,m=Sm′(xl+1)=Uxl+1,m=Ul+1,m-Ul,mh+h6[Ml,m+2Ml+1,m],
and from (2.6), (2.9)ml-1,m=Sm′(xl-1)=Uxl-1,m=Ul,m-Ul-1,mh-h6[Ml,m+2Ml-1,m],
We consider the following approximations:(2.10)U̅yyl,m=(Ul,m+1-2Ul,m+Ul,m-1)(k2)=Uyyl,m+O(k2),U̅yyl+1,m=(Ul+1,m+1-2Ul+1,m+Ul+1,m-1)(k2)=Uyyl+1,m+O(k2),U̅yyl-1,m=(Ul-1,m+1-2Ul-1,m+Ul-1,m-1)(k2)=Uyyl-1,m+O(k2),U̅xl,m=(Ul+1,m-Ul-1,m)(2h)=Uxl,m+h26U30+O(h4),U̅xl+1,m=(3Ul+1,m-4Ul,m+Ul-1,m)(2h)=Uxl+1,m-h23U30+O(h3),U̅xl-1,m=(-3Ul-1,m+4Ul,m-Ul+1,m)(2h)=Uxl-1,m-h23U30+O(h3).Since the derivative values ofSm(x) defined by (2.4), (2.7), (2.8), and (2.9) are not known at each grid point (xl,ym), we use the following approximations for the derivatives of Sm(x).Let(2.11)M̅l,m=-U̅yyl,m+DlU̅xl,m+gl,m,M̅l+1,m=-U̅yyl+1,m+Dl+1U̅xl+1,m+gl+1,m,M̅l-1,m=-U̅yyl-1,m+Dl-1U̅xl-1,m+gl-1,m,U̅̅xl+1,m=Ul+1,m-Ul,mh+h6[M̅l,m+2M̅l+1,m],U̅̅xl-1,m=Ul,m-Ul-1,mh-h6[M̅l,m+2M̅l-1,m],F̅l+1,m=Dl+1U̅xl+1,m+gl+1,m,F̅l-1,m=Dl-1U̅xl-1,m+gl-1,m,Ûxl,m=U̅xl,m-h12[F̅l+1,m-F̅l-1,m]+h12[U̅yyl+1,m-U̅yyl-1,m],F̂l,m=DlÛxl,m+gl,m.Then at each grid point(xl,ym), a cubic spline finite difference method of Numerov type with accuracy of O(k2+h4) for the solution of differential equation (1.1) may be written as(2.12)λ2δx2Ul,m+k212[U̅yyl+1,m+U̅yyl-1,m+10U̅yyl,m]=k212[Dl+1U̅̅xl+1,m+Dl-1U̅̅xl-1,m+10DlÛxl,m]+k212[gl+1,m+gl-1,m+10gl,m]+Tl,m,
where the local truncation error Tl,m=O(k4+k2h4). Note that the method (2.12) is of O(k2+h4) for the numerical solution of (1.1). However, the method (2.12) fails to compute at l=1, when D(x) and/or g(x,y) contains the singular terms like 1/x,1/x2, and so forth. For example, if D(x)=1/x, this implies Dl-1=1/xl-1, which cannot be evaluated at l=1 (since x0=0). We modify the method (2.12) in such a manner, so that the method retains its order and accuracy everywhere in the vicinity of the singularity.We use the following approximation:(2.13)Dl+1=Dl+hD10+h22D20+O(h3),Dl-1=Dl-hD10+h22D20-O(h3),gl+1,m=gl,m+hg10+h22g20+O(h3),gl-1,m=gl,m-hg10+h22g20-O(h3).Substituting the approximations (2.13) into (2.12) and neglecting higher-order terms and local truncation error, we get(2.14)(12+δx2)δy2ul,m+12λ2δx2ul,m=λ2h2(12D00+h2D20)(2μxδx)ul,m+λ2h2(D10-D002)δx2ul,m-h2D10δy2ul,m+hD002(δy22μxδx)ul,m+k2[12g00+h2(g20+D10g00-D00g10)],
where μxul=(1/2)(ul+1/2+ul-1/2) and δxul=(ul+1/2-ul-1/2). Note that the cubic spline method (2.14) is of O(k2+h4) for the numerical solution of elliptic equation (1.1), which is also free from the terms 1/xl±1 and hence can be computed for l=1(1)N,m=1(1)M.


## 3. Iterative Method

Now consider the convection-diffusion equation(3.1)∂2u∂x2+∂2u∂y2=β∂u∂x,0<x,y<1,
where β>0 is a constant, and magnitude of β determines the ratio of convection to diffusion term. Substituting D(x)=β and g(x,y)=0 into the difference scheme (2.14) and simplifying, we obtain a nine-point cubic spline difference scheme of O(k2+h4) accuracy for the solution of the convection-diffusion equation (3.1) given by(3.2)α0ul,m+α1ul+1,m+α2ul-1,m+α3ul,m+1+α4ul,m-1+α5ul+1,m+1+α6ul+1,m-1+α7ul-1,m+1+α8ul-1,m-1=0,
where R=βh/2 is called the cell Reynold number, and the coefficients αj,j=0,1,2,…,8 are defined by(3.3)α0=20+24λ2(1+R23),α1=-12λ2(1-R+R23)+2(1-R),α2=-12λ2(1+R+R23)+2(1+R),α3=α4=-10,α5=α6=-(1-R),α7=α8=-(1+R).The scheme (3.2) may be written in matrix form(3.4)Au=b,
where A is a square matrix of order NM×NM, u is the solution vector, and b is the right hand side vector consisting of boundary values.The coefficient matrixA has a block tridiagonal structureA=tri[-L,D,-U], with the submatrices -L,D, and -U given by(3.5)-L=tri[α8,α4,α6]=-U,D=tri[α2,α0,α1].We focus on line stationary iterative methods for solving the linear system (3.5). The coefficient matrix A can be written as A=D-L-U, where D is block tri-diagonal matrix of A, -L is strictly block lower triangular part, and -U is strictly block upper triangular part of matrix A. The iteration matrices of the block Jacobi and block Gauss-Seidel methods are described by
(3.6)GJ=D-1(L+U),GGS=(D-L)-1U.The matrixA has block tri-diagonal form and hence is block consistently ordered (see Varga [22]).It can be verified thatα0>0 and αj<0 for j=1,2,…,8 provided |R|<1. One can also verify that
(3.7)α0=∑i=18|αi|,
which implies that A is weakly diagonally dominant. Since A is reducible, we conclude that it is also an M-matrix (see Varga [22]).Applying the Jacobi iterative method to the scheme (3.2), we get the iterative scheme(3.8)(1-R)ul+1,m+1(k)+[12λ2(1-R+R23)-2(1-R)]ul+1,m(k)+(1-R)ul+1,m-1(k)+10ul,m+1(k)-[20+24λ2(1+R23)]ul,m(k+1)+10ul,m-1(k)+(1+R)ul-1,m+1(k)+[12λ2(1+R+R23)-2(1+R)]ul-1,m(k)+(1+R)ul-1,m-1(k)=0.The propagating factors for the Jacobi iterative methods are given by(3.9)ξj=15+6λ2(1+(R2/3))[5cos⁡⁡(πk)+(6λ2(1+R+R23)-(1+R)(1-cos⁡kπ))×(6λ2(1-R+R23)-(1-R)(1-cos⁡kπ))cos⁡⁡πh],
where (M+1)h=1,(N+1)k=1. Consequently, the spectral radii ρ of the Jacobi and Gauss-Seidel matrices are related by
(3.10)ρ(GGS)=ρ(GJ)2.Hence, the associated iteration(3.11)u(k+1)=Gu(k)+c
converges for any initial guess, where G is either Jacobi or Gauss-Seidel iteration matrix. The Jacobi and Gauss-Seidel splitting are regular for both the line and point versions, and hence, they converge for any initial guess.


## 4. Numerical Results

If we replace the partial derivatives in (1.1) by the central difference approximations at the grid point (l,m), we obtain a central difference scheme (CDS) which is of O(k2+h2). We now solve the following two benchmark problems whose exact solutions are known. The right hand side homogeneous function and boundary conditions may be obtained by using the exact solution as a test procedure. We use block Gauss-Seidel iterative method (see [22–24]) to solve the proposed scheme (2.14). In all cases, we have considered u(0)=0 as the initial guess, and the iterations were stopped when the absolute error tolerance |u(k+1)-u(k)|≤10-12 was achieved. In all cases, we have calculated maximum absolute errors (l∞-norm) for different grid sizes. All computations were performed using double-precision arithmetic.Example 4.1 (Convection-Diffusion Equation).
The problem is to solve (3.1) in the solution region 0<x,y<1 whose exact solution is given by
(4.1)u(x,y)=eβx/2sin⁡⁡πysinh⁡⁡σ[2e-β/2sinh⁡⁡σx+sinh⁡⁡σ(1-x)],
where σ2=π2+β2/4, β>0.
The maximum absolute errors foru are tabulated in Tables 1 and 2.Table 1
Example4.1: The maximum absolute errors.
(h,k)ProposedO(k2+h4) methodO(k2+h2) methodβ=10β=50β=10β=50(120,110)0.2645E−020.1982E−010.9212E−020.2240E+00(140,120)0.6560E−030.1823E−020.2279E−020.6511E−01(120,140)0.2167E−030.1875E−010.1125E−010.2253E+00(180,140)0.1635E−030.1511E−030.5682E−030.1416E−01(140,180)0.4352E−040.1604E−020.2755E−020.6534E−01Table 2
Example4.1: The maximum absolute errors (k/h2=64).
hProposedO(k2+h4) methodO(k4+h4) method discussed in [8]β=5β=10β=15β=5β=10β=151160.1808E−010.1636E−010.1494E−010.4408E−010.4212E−010.4018E−011320.1129E−020.1026E−020.9396E−030.2692E−020.2554E−020.2482E−021640.7054E−040.6419E−040.5881E−040.1645E−030.1566E−030.1512E−03Example 4.2 (Poisson’s Equation in Polar Cylindrical Coordinates).
Consider the following:(4.2)∂2u∂r2+∂2u∂z2+1r∂u∂r=f(r,z),0<r,z<1.
The exact solution is given byu(r,z)=r2sinhrcoshz. The maximum absolute errors for u are tabulated in Tables 3 and 4.Table 3
Example4.2: The maximum absolute errors.
(h, k)ProposedO(k2+h4) methodO(k2+h2) method(120,110)0.3986E−030.1133E−02(140,120)0.1089E−030.3003E−03(120,140)0.4117E−030.1132E−02(180,140)0.2893E−040.7836E−04(140,180)0.1116E−030.2984E−03Table 4
Example4.2: The maximum absolute errors (k/h2=64).
hProposedO(k2+h4) methodO(k4+h4) method discussed in [8]1160.2711E−020.4842E−021320.1654E−030.2988E−031640.1016E−040.1811E−04


## 5. Concluding Remarks

The available numerical methods based on spline approximations for the numerical solution of 2D Poisson’s equation are ofO(k2+h2) accurate, which require nine grid points. In this paper, using the same number of grid points, we have discussed a new stable compact nine point cubic spline finite difference method of O(k2+h4) accuracy for the solution of Poisson’s equation in polar cylindrical coordinates. For a fixed parameter γ=k/h2, the proposed method behaves like a fourth-order method, which is exhibited from the computed results.



---

*Source: 234516-2012-02-12.xml*