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

Abstract

Using nonpolynomial cubic spline approximation in 𝑥- and finite difference in 𝑦-direction, we discuss a numerical approximation of 𝑂(𝑘²+ℎ⁴) for the solutions of diffusion-convection equation, where 𝑘>0 and ℎ>0 are grid sizes in 𝑦- and 𝑥-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.

1. Introduction

We consider the two-dimensional elliptic equation of the form𝜕²𝑢𝜕𝑥²+𝜕²𝑢𝜕𝑦²=𝐷(𝑥)𝜕𝑢𝜕𝑥+𝑔(𝑥,𝑦),(𝑥,𝑦)∈Ω.(1.1) The Dirichlet boundary conditions are given by𝑢(𝑥,𝑦)=𝑔(𝑥,𝑦),(𝑥,𝑦)∈Γ,(1.2) where Ω={(𝑥,𝑦)∣0<𝑥,𝑦<1}is the solution domain and Γ is its boundary. For 𝑔(𝑥,𝑦)=0 and 𝐷(𝑥)=𝛽, the above equation represents diffusion-convention equation. For 𝐷(𝑥)=−1/𝑥, the above equations represent the Poisson’s equation with singular coefficients in rectangular coordinates. Similarly, for 𝐷(𝑥)=−1/𝑥 and replacing the variables 𝑥,𝑦 by 𝑟,𝑧, 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 ℎ⁴ 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 𝑦- and four in 𝑥-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 ℎ>0 and 𝑘>0 in the directions 𝑥- and 𝑦- respectively, so that the mesh points (𝑥_𝑙,𝑦_𝑚) denoted by (𝑙,𝑚) are defined by 𝑥_𝑙=𝑙ℎ and 𝑦_𝑚=𝑚𝑘;𝑙=0,1,…,𝑁+1;𝑚=0,1,…,𝑀+1, where 𝑁 and 𝑀 are positive integers such that (𝑁+1)ℎ=1 and (𝑀+1)𝑘=1. The mesh ratio parameter is denoted by 𝜆=(𝑘/ℎ). The notations 𝑢_𝑙,𝑚and𝑈_𝑙,𝑚 are used for the discrete approximation and the exact value of 𝑢(𝑥,𝑦) at the grid point (𝑥_𝑙,𝑦_𝑚), respectively. Similarly, we write 𝐷(𝑥_𝑙)=𝐷_𝑙 and 𝐷(𝑥_𝑙±1)=𝐷_𝑙±1.

At the grid point (𝑥_𝑙,𝑦_𝑚), we denote 𝑊_𝑎𝑏=𝜕^𝑎+𝑏𝑊𝜕𝑥_𝑙𝑎𝜕𝑦_𝑚𝑏,𝑊=𝑈,𝐷and𝑔.(2.1)

Let 𝑆_𝑚(𝑥) be the cubic spline interpolating polynomial of the function 𝑢(𝑥,𝑦_𝑚) between the grid point (𝑥_𝑙−1,𝑦_𝑚) and (𝑥_𝑙,𝑦_𝑚) and is given by𝑆_𝑚(𝑥)=𝑥^𝑙−𝑥³6ℎ𝑀^{𝑙−1,𝑚}+𝑥−𝑥^𝑙−1³6ℎ𝑀^𝑙,𝑚+𝑈_{𝑙−1,𝑚}−ℎ²6𝑀^{𝑙−1,𝑚}𝑥^𝑙−𝑥+𝑈_𝑙,𝑚ℎ−ℎ²6𝑀^𝑙,𝑚⋅𝑥−ℎ𝑥^𝑙−1,𝑥^𝑙−1≤𝑥≤𝑥_𝑙;𝑙=1,2,…,𝑁+1;𝑚=0,1,2,…,𝑀+1(2,.2) which satisfies at 𝑚th-line parallel to 𝑥-axis the following properties:(i)𝑆_𝑚(𝑥) coincides with a polynomial of degree three on each [𝑥_𝑙−1,𝑥_𝑙],𝑙=1,2,…,𝑁+1and𝑚=1,2,…𝑀,(ii)𝑆_𝑚(𝑥)∈𝐶²[0,1], (iii)𝑆_𝑚(𝑥_𝑙)=𝑈_𝑙,𝑚,𝑙=0,1,2,…,𝑁+1and𝑚=1,2,…,𝑀.

The derivatives of cubic spline function 𝑆_𝑚(𝑥) are given by𝑆^′_𝑚(𝑥)=−𝑥^𝑙−𝑥²2ℎ𝑀^{𝑙−1,𝑚}+𝑥−𝑥^𝑙−1²2ℎ𝑀^𝑙,𝑚+𝑈^𝑙,𝑚−𝑈_{𝑙−1,𝑚}ℎ−ℎ6𝑀^𝑙,𝑚𝑆_𝑚−𝑀_{𝑙−1,𝑚},(𝑥)=𝑥^𝑙−𝑥ℎ𝑀^{𝑙−1,𝑚}+𝑥−𝑥^𝑙−1ℎ𝑀^𝑙,𝑚,(2.3) where𝑀_𝑙,𝑚=𝑆^𝑚𝑥_𝑙=𝑈_{𝑥𝑥𝑙,𝑚}=−𝑈_{𝑦𝑦𝑙,𝑚}+𝐷_𝑙𝑈_{𝑥𝑙,𝑚}+𝑔_𝑙,𝑚,𝑙=0,1,2,…,𝑁+1,𝑗=1,2,…,𝐽,(2.4)𝑚_𝑙,𝑚=𝑆^𝑚𝑥_𝑙=𝑈_{𝑥𝑙,𝑚}=𝑈^𝑙,𝑚−𝑈_{𝑙−1,𝑚}ℎ+ℎ6𝑀^{𝑙−1,𝑚}+2𝑀_𝑙,𝑚,𝑥_𝑙−1≤𝑥≤𝑥_𝑙,(2.5) and replacing ℎ by “−ℎ”, we get𝑚_𝑙,𝑚=𝑆^𝑚𝑥_𝑙=𝑈_{𝑥𝑙,𝑚}=𝑈^𝑙+1,𝑚−𝑈_𝑙,𝑚ℎ−ℎ6𝑀^𝑙+1,𝑚+2𝑀_𝑙,𝑚,𝑥_𝑙≤𝑥≤𝑥_𝑙+1,(2.6) Combining (2.5) and (2.6), we obtain𝑚_𝑙,𝑚=𝑆^𝑚𝑥_𝑙=𝑈_{𝑥𝑙,𝑚}=𝑈^𝑙+1,𝑚−𝑈_{𝑙−1,𝑚}2ℎ−1ℎ2𝑀^𝑙+1,𝑚−𝑀_{𝑙−1,𝑚},(2.7) Further, from (2.5), we have𝑚_𝑙+1,𝑚=𝑆^𝑚𝑥_𝑙+1=𝑈_{𝑥𝑙+1,𝑚}=𝑈^𝑙+1,𝑚−𝑈_𝑙,𝑚ℎ+ℎ6𝑀^𝑙,𝑚+2𝑀_𝑙+1,𝑚,(2.8) and from (2.6), 𝑚_{𝑙−1,𝑚}=𝑆^𝑚𝑥_𝑙−1=𝑈_{𝑥𝑙−1,𝑚}=𝑈^𝑙,𝑚−𝑈_{𝑙−1,𝑚}ℎ−ℎ6𝑀^𝑙,𝑚+2𝑀_{𝑙−1,𝑚},(2.9) We consider the following approximations:𝑈_{𝑦𝑦𝑙,𝑚}=𝑈^𝑙,𝑚+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,𝑚}(2ℎ)=𝑈^{𝑥𝑙+1,𝑚}−ℎ²3𝑈³⁰+𝑂ℎ³,𝑈_{𝑥𝑙−1,𝑚}=−3𝑈^{𝑙−1,𝑚}+4𝑈_𝑙,𝑚−𝑈_𝑙+1,𝑚(2ℎ)=𝑈^{𝑥𝑙−1,𝑚}−ℎ²3𝑈³⁰+𝑂ℎ³.(2.10)

Since the derivative values of 𝑆_𝑚(𝑥) defined by (2.4), (2.7), (2.8), and (2.9) 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,𝑚−𝑈_𝑙,𝑚ℎ+ℎ6𝑀_𝑙,𝑚+2𝑀_𝑙+1,𝑚,𝑈_{𝑥𝑙−1,𝑚}=𝑈^𝑙,𝑚−𝑈_{𝑙−1,𝑚}ℎ−ℎ6𝑀_𝑙,𝑚+2𝑀_{𝑙−1,𝑚},𝐹_𝑙+1,𝑚=𝐷_𝑙+1𝑈_{𝑥𝑙+1,𝑚}+𝑔_𝑙+1,𝑚,𝑈_{𝑥𝑙,𝑚}𝐹_{𝑙−1,𝑚}=𝐷_𝑙−1𝑈_{𝑥𝑙−1,𝑚}+𝑔_{𝑙−1,𝑚},=𝑈_{𝑥𝑙,𝑚}−ℎ12𝐹_𝑙+1,𝑚−𝐹_{𝑙−1,𝑚}+ℎ12𝐹_𝑙,𝑚𝑈_{𝑦𝑦𝑙+1,𝑚}−𝑈_{𝑦𝑦𝑙−1,𝑚},=𝐷_𝑙𝑈_{𝑥𝑙,𝑚}+𝑔_𝑙,𝑚.(2.11)

Then at each grid point (𝑥_𝑙,𝑦_𝑚), a cubic spline finite difference method of Numerov type with accuracy of 𝑂(𝑘²+ℎ⁴) for the solution of differential equation (1.1) may be written as𝜆²𝛿^2𝑥𝑈_𝑙,𝑚+𝑘²12=𝑘₂𝑈_{𝑦𝑦𝑙+1,𝑚}+𝑈_{𝑦𝑦𝑙−1,𝑚}+10𝑈_{𝑦𝑦𝑙,𝑚}12𝐷^𝑙+1𝑈_{𝑥𝑙+1,𝑚}+𝐷_𝑙−1𝑈_{𝑥𝑙−1,𝑚}+10𝐷_𝑙𝑈_{𝑥𝑙,𝑚}+1𝑘2²𝑔^𝑙+1,𝑚+𝑔_{𝑙−1,𝑚}+10𝑔_𝑙,𝑚+𝑇_𝑙,𝑚,(2.12) where the local truncation error 𝑇_𝑙,𝑚=𝑂(𝑘⁴+𝑘²ℎ⁴). Note that the method (2.12) is of 𝑂(𝑘²+ℎ⁴) for the numerical solution of (1.1). However, the method (2.12) fails to compute at 𝑙=1, when 𝐷(𝑥) and/or 𝑔(𝑥,𝑦) contains the singular terms like 1/𝑥,1/𝑥², and so forth. For example, if 𝐷(𝑥)=1/𝑥, this implies 𝐷_𝑙−1=1/𝑥_𝑙−1, which cannot be evaluated at 𝑙=1 (since 𝑥₀=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:𝐷_𝑙+1=𝐷_𝑙+ℎ𝐷₁₀+ℎ²2𝐷²⁰𝐷_𝑙−1+𝑂ℎ³,=𝐷_𝑙−ℎ𝐷₁₀+ℎ²2𝐷²⁰𝑔_𝑙+1,𝑚−𝑂ℎ³,=𝑔_𝑙,𝑚+ℎ𝑔₁₀+ℎ²2𝑔²⁰𝑔_{𝑙−1,𝑚}+𝑂ℎ³,=𝑔_𝑙,𝑚−ℎ𝑔₁₀+ℎ²2𝑔²⁰−𝑂ℎ³.(2.13)

Substituting the approximations (2.13) into (2.12) and neglecting higher-order terms and local truncation error, we get12+𝛿^2𝑥𝛿^2𝑦𝑢_𝑙,𝑚+12𝜆²𝛿^2𝑥𝑢_𝑙,𝑚=𝜆²ℎ212𝐷⁰⁰+ℎ²𝐷₂₀2𝜇_𝑥𝛿_𝑥𝑢_𝑙,𝑚+𝜆²ℎ²𝐷₁₀−𝐷₂₀₀𝛿^2𝑥𝑢_𝑙,𝑚−ℎ²𝐷₁₀𝛿^2𝑦𝑢_𝑙,𝑚+ℎ𝐷⁰⁰2𝛿^2𝑦2𝜇_𝑥𝛿_𝑥𝑢_𝑙,𝑚+𝑘²12𝑔₀₀+ℎ²𝑔₂₀+𝐷₁₀𝑔₀₀−𝐷₀₀𝑔₁₀,(2.14) where 𝜇_𝑥𝑢_𝑙=(1/2)(𝑢_𝑙+1/2+𝑢_𝑙−1/2) and 𝛿_𝑥𝑢_𝑙=(𝑢_𝑙+1/2−𝑢_𝑙−1/2). Note that the cubic spline method (2.14) is of 𝑂(𝑘²+ℎ⁴) for the numerical solution of elliptic equation (1.1), which is also free from the terms 1/𝑥_𝑙±1 and hence can be computed for 𝑙=1(1)𝑁,𝑚=1(1)𝑀.

3. Iterative Method

Now consider the convection-diffusion equation 𝜕²𝑢𝜕𝑥²+𝜕²𝑢𝜕𝑦²=𝛽𝜕𝑢𝜕𝑥,0<𝑥,𝑦<1,(3.1) where 𝛽>0 is a constant, and magnitude of 𝛽 determines the ratio of convection to diffusion term. Substituting 𝐷(𝑥)=𝛽 and 𝑔(𝑥,𝑦)=0 into the difference scheme (2.14) and simplifying, we obtain a nine-point cubic spline difference scheme of 𝑂(𝑘²+ℎ⁴) accuracy for the solution of the convection-diffusion equation (3.1) given by𝛼₀𝑢_𝑙,𝑚+𝛼₁𝑢_𝑙+1,𝑚+𝛼₂𝑢_{𝑙−1,𝑚}+𝛼₃𝑢_𝑙,𝑚+1+𝛼₄𝑢_{𝑙,𝑚−1}+𝛼₅𝑢_{𝑙+1,𝑚+1}+𝛼₆𝑢_{𝑙+1,𝑚−1}+𝛼₇𝑢_{𝑙−1,𝑚+1}+𝛼₈𝑢_{𝑙−1,𝑚−1}=0,(3.2) where 𝑅=𝛽ℎ/2 is called the cell Reynold number, and the coefficients 𝛼_𝑗,𝑗=0,1,2,…,8 are defined by𝛼₀=20+24𝜆²1+𝑅3²,𝛼¹=−12𝜆²𝛼1−𝑅+𝑅3²+2(1−𝑅),₂=−12𝜆²1+𝑅+𝑅3²+2(1+𝑅),𝛼³𝛼=𝛼₄=−10,₅=𝛼₆=−(1−𝑅),𝛼₇=𝛼₈=−(1+𝑅).(3.3)

The scheme (3.2) may be written in matrix form𝐀𝐮=𝐛,(3.4) where 𝐀 is a square matrix of order 𝑁𝑀×𝑁𝑀, 𝐮 is the solution vector, and 𝐛 is the right hand side vector consisting of boundary values.

The coefficient matrix 𝐀 has a block tridiagonal structure𝐀=tri[−𝐋,𝐃,−𝐔], with the submatrices −𝐋,𝐃, and −𝐔 given by−𝐋=tri𝛼₈,𝛼₄,𝛼₆=−𝐔,𝐃=tri𝛼₂,𝛼₀,𝛼₁.(3.5)

We focus on line stationary iterative methods for solving the linear system (3.5). The coefficient matrix 𝐀 can be written as 𝐀=𝐃−𝐋−𝐔, where 𝐃 is block tri-diagonal matrix of 𝐀, −𝐋 is strictly block lower triangular part, and −𝐔 is strictly block upper triangular part of matrix 𝐀. The iteration matrices of the block Jacobi and block Gauss-Seidel methods are described by 𝐆_𝐽=𝐃⁻¹(𝐋+𝐔),𝐆_GS=(𝐃−𝐋)⁻¹𝐔.(3.6)

The matrix 𝐀 has block tri-diagonal form and hence is block consistently ordered (see Varga [22]).

It can be verified that 𝛼₀>0 and 𝛼_𝑗<0 for 𝑗=1,2,…,8 provided |𝑅|<1. One can also verify that 𝛼₀=⁸_𝑖=1||𝛼_𝑖||,(3.7) which implies that 𝐀 is weakly diagonally dominant. Since 𝐀 is reducible, we conclude that it is also an 𝐌-matrix (see Varga [22]).

Applying the Jacobi iterative method to the scheme (3.2), we get the iterative scheme(1−𝑅)𝑢^(𝑘)_{𝑙+1,𝑚+1}+12𝜆²1−𝑅+𝑅3²−2(1−𝑅)𝑢^{(𝑘)𝑙+1,𝑚}+(1−𝑅)𝑢^(𝑘)_{𝑙+1,𝑚−1}+10𝑢^(𝑘)_𝑙,𝑚+1−20+24𝜆²1+𝑅3²𝑢^{(𝑘+1)𝑙,𝑚}+10𝑢^(𝑘)_{𝑙,𝑚−1}+(1+𝑅)𝑢^(𝑘)_{𝑙−1,𝑚+1}+12𝜆²1+𝑅+𝑅3²−2(1+𝑅)𝑢^{(𝑘)𝑙−1,𝑚}+(1+𝑅)𝑢^(𝑘)_{𝑙−1,𝑚−1}=0.(3.8)

The propagating factors for the Jacobi iterative methods are given by𝜉_𝑗=15+6𝜆²1+𝑅²/3⎡⎢⎢⎣5cos(𝜋𝑘)+⎛⎜⎜⎝6𝜆²×⎛⎜⎜⎝1+𝑅+𝑅3²−(1+𝑅)(1−cos𝑘𝜋)⎞⎟⎟⎠6𝜆²1−𝑅+𝑅3²−(1−𝑅)(1−cos𝑘𝜋)⎞⎟⎟⎠cos𝜋ℎ⎤⎥⎥⎦,(3.9) where (𝑀+1)ℎ=1,(𝑁+1)𝑘=1. Consequently, the spectral radii 𝜌 of the Jacobi and Gauss-Seidel matrices are related by 𝜌𝐆_GS=𝜌(𝐆_𝐽)².(3.10)

Hence, the associated iteration𝐮^(𝑘+1)=𝐆𝐮^(𝑘)+𝐜(3.11) converges for any initial guess, where 𝐆 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 (𝑙,𝑚), we obtain a central difference scheme (CDS) which is of 𝑂(𝑘²+ℎ²). 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 𝐮⁽⁰⁾=0 as the initial guess, and the iterations were stopped when the absolute error tolerance |𝐮^(𝑘+1)−𝐮^(𝑘)|≤10⁻¹² was achieved. In all cases, we have calculated maximum absolute errors (𝑙_∞-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<𝑥,𝑦<1 whose exact solution is given by 𝑢(𝑥,𝑦)=𝑒^𝛽𝑥/2sin𝜋𝑦sinh𝜎2𝑒^−𝛽/2sinh𝜎𝑥+sinh𝜎(1−𝑥),(4.1) where 𝜎²=𝜋²+𝛽²/4, 𝛽>0.
The maximum absolute errors for 𝑢 are tabulated in Tables 1 and 2.

Table 1 

Example 4.1: The maximum absolute errors.

Table 2 

Example 4.1: The maximum absolute errors (𝑘/ℎ²=64).

Example 4.2 (Poisson’s Equation in Polar Cylindrical Coordinates). Consider the following: 𝜕²𝑢𝜕𝑟²+𝜕²𝑢𝜕𝑧²+1𝑟𝜕𝜕𝑢𝑟=𝑓(𝑟,𝑧),0<𝑟,𝑧<1.(4.2)
The exact solution is given by 𝑢(𝑟,𝑧)=𝑟²sinh𝑟cosh𝑧. The maximum absolute errors for 𝑢 are tabulated in Tables 3 and 4.

Table 3 

Example 4.2: The maximum absolute errors.

Table 4 

Example 4.2: The maximum absolute errors (𝑘/ℎ²=64).

5. Concluding Remarks

The available numerical methods based on spline approximations for the numerical solution of 2D Poisson’s equation are of 𝑂(𝑘²+ℎ²) 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 𝑂(𝑘²+ℎ⁴) accuracy for the solution of Poisson’s equation in polar cylindrical coordinates. For a fixed parameter 𝛾=𝑘/ℎ², the proposed method behaves like a fourth-order method, which is exhibited from the computed results.