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

Using nonpolynomial cubic spline approximation in xand finite difference in y-direction, we discuss a numerical approximation ofO k2 h4 for the solutions of diffusion-convection equation, where k > 0 and h > 0 are grid sizes in yand 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. Iterativemethod of the proposedmethod is discussed, and numerical examples are given in support of the theoretical results.


Introduction
We consider the two-dimensional elliptic equation of the form The Dirichlet boundary conditions are given by u x, y g x, y , x, y ∈ Γ, 1.2 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 secondorder 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 h 4 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  .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 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.

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 x l , y m denoted by l, m are defined by x l lh and y m 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 u l,m and U l,m are used for the discrete approximation and the exact value of u x, y at the grid point x l , y m , respectively.Similarly, we write D x l D l and D x l±1 D l±1 .At the grid point x l , y m , we denote W ab ∂ a b W ∂x l a ∂y m b , W U, D and g.

2.1
Let S m x be the cubic spline interpolating polynomial of the function u x, y m between the grid point x l−1 , y m and x l , y m and is given by which satisfies at mth-line parallel to x-axis the following properties: i S m x coincides with a polynomial of degree three on each iii S m x l U l,m , l 0, 1, 2, . . ., N 1 and m 1, 2, . . ., M.

ISRN Mathematical Physics
The derivatives of cubic spline function S m x are given by where and replacing h by "−h", we get Combining 2.5 and 2.6 , we obtain Further, from 2.5 , we have and from 2.6 , We consider the following approximations:

2.10
Since the derivative values of S m x defined by 2.4 , 2.7 , 2.8 , and 2.9 are not known at each grid point x l , y m , we use the following approximations for the derivatives of S m x . Let

2.11
Then at each grid point x l , y m , a cubic spline finite difference method of Numerov type with accuracy of O k 2 h 4 for the solution of differential equation 1.1 may be written as where the local truncation error T l,m O k 4 k 2 h 4 .Note that the method 2.12 is of O k 2 h 4 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/x 2 , and so forth.For example, if D x 1/x, this implies D l−1 1/x l−1 , which cannot be evaluated at l 1 since x 0 0 .

ISRN Mathematical Physics
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
Substituting the approximations 2.13 into 2.12 and neglecting higher-order terms and local truncation error, we get where μ x u l 1/2 u l 1/2 u l−1/2 and δ x u l u l 1/2 − u l−1/2 .Note that the cubic spline method 2.14 is of O k 2 h 4 for the numerical solution of elliptic equation 1.1 , which is also free from the terms 1/x l±1 and hence can be computed for l 1 1 N, m 1 1 M.

Iterative Method
Now consider the convection-diffusion equation 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 k 2 h 4 accuracy for the solution of the convection-diffusion equation 3.1 given by where R βh/2 is called the cell Reynold number, and the coefficients α j , j 0, 1, 2, . . ., 8 are defined by

3.3
The scheme 3.2 may be written in matrix form Au b, 3.4 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 matrix A has a block tridiagonal structure A tri −L, D, −U , with the submatrices −L, D, and −U given by 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
The matrix A 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 which implies that A is weakly diagonally dominant.Since A is reducible, we conclude that it is also an M-matrix see Varga 22 .

ISRN Mathematical Physics
Applying the Jacobi iterative method to the scheme 3.2 , we get the iterative scheme

3.8
The propagating factors for the Jacobi iterative methods are given by 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
Hence, the associated iteration 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.

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 k 2 h 2 .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.The maximum absolute errors for u are tabulated in Tables 1 and 2. ∂ 2 u ∂z 2   1 r ∂u ∂r f r, z , 0 < r, z < 1.

4.2
The exact solution is given by u r, z r 2 sinh r cosh z.The maximum absolute errors for u are tabulated in Tables 3 and 4.

Concluding Remarks
The available numerical methods based on spline approximations for the numerical solution of 2D Poisson's equation are of O k 2 h 2 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 k 2 h 4 accuracy for the solution of Poisson's equation in polar cylindrical coordinates.For a fixed parameter γ k/h 2 , the proposed method behaves like a fourth-order method, which is exhibited from the computed results.

Table 4 :
Example 4.2: The maximum absolute errors k/h 2 64 .Poisson's Equation in Polar Cylindrical Coordinates .Consider the following: