Tension spline method for the solution of elliptic equations

In this paper, two classes of methods are developed for the solution of two-dimensional elliptic partial differential equations. We have used tension spline function approximation in both x and y spatial directions and a new scheme of order has been obtained. The convergence analysis of the methods has been carried out. Numerical examples are given to illustrate the applicability and accurate nature of our approach.


Introduction
Elliptic partial differential equations often arise in many fields of engineering and theoretical physics such as electric potential in electrostatics, electromagnetic scattering theory [1], micro and nano electronic devise physics, membrane vibration and radar scattering [2,3].
We consider the linear second-order two-dimensional elliptic equation of the form subjected to Dirichlet boundary conditions u(x, y) = g(x, y), (x, y) ∈ ∂ , where A, B and C are positive constants and = {(x, y); 0 < x, y < 1} with the boundary ∂ . We assume that the solution u(x, y) and forcing function f are sufficiently smooth and have required continuous partial derivatives. The well-known Poisson equation and modified Helmholtz equation are special cases of Equation (1).
Various numerical methods have been discussed and developed for solving elliptic partial differential equations. For the solution of two-dimensional Poisson equation, fourth-order compact difference schemes have discussed in [4,5] and also a new family of high-order finite difference schemes is proposed in [6] by using implicit finite difference formulas. Sixthorder compact scheme is developed for the Poisson equation in [7] based on the fourth-order nine-point difference approximation, combined with Richardson extrapolation. Fourth-order finite difference methods for the solution of Helmholtz equation are presented in [8,9]. Cubic spline collocation method was introduced for Poisson equation by Romenski [10]. Bialecki used orthogonal spline collocation method for the solution of Poisson equation [11]. The method of collocation at nodal points based on the tensor product of cubic splines for Helmholtz equation with Dirichlet boundary conditions was first analysed independently by Cavendish and Ito in their PhD theses and second-order convergence of the method was proved. Christara has developed quadratic spline collocation methods for elliptic problems in the unit square [12]. These methods which have third order of accuracy, collocate a perturbed differential equation. Houstis et al. analysed a new class of collocation methods using cubic spline for the general form of elliptic PDEs and derived fourth order of accuracy [13]. Quadratic spline collocation methods for two-dimensional Helmholtz problem are developed in [14], it is shown that the methods are fourth-order accurate at the nodes and collocation points. Rashidinia et al. proposed the cubic spline method for the singular elliptic boundary value problems in [15,16]. Islam Khan et al. used tension spline method to solve second-order singularly perturbed boundary-value problems [17] also tension spline method for the solution of nonlinear partial differential equations have discussed in [18,19]. Mohanty et al. have used a new finite difference scheme for twodimensional elliptic problems [20], also they proposed cubic spline finite difference method of order O(k 2 + h 4 ) for two-dimensional quasi-linear elliptic equations [21]. Higher-order accurate finite volume schemes for two-dimensional Helmholtz equations are developed in [22]. For numerically solving inhomogeneous elliptic partial differential equations, a Chebyshev polynomial scheme combined with the method of fundamental solutions (MFS) and the equilibrated collocation Trefftz method have presented in [23] and a particular solution using the standard polynomial basis function has been derived in [24]. A new wavelet multigrid method based on Daubechies filter coefficients is proposed in [25] for the numerical solution of elliptic type equations.
In this paper, we have presented the non-polynomial tension spline method for two-dimensional linear elliptic equation (1). We discussed the derivation of nonpolynomial tension spline approximation and by using tension spline approximation in both spatial directions, two classes of nine-point scheme have been obtained. The truncation error of the scheme is carried out and convergence of our approach is proved.

Non-polynomial spline functions
We consider the solution domain = [0, 1] × [0, 1] and choose grid spacing h > 0 and k > 0, in the xdirection and y-direction, respectively. The mesh points are defined by (x i , y j ) = (ih, jk), i = 0, 1, . . . , m + 1, j = 0, 1, . . . , n + 1, where m and n are positive integers. In each segment [x i , x i+1 ] let s 1j (x) be the tension spline function which interpolates u(x, y j ) and define where a 1i , b 1i , c 1i and d 1i are unknown coefficients and λ 1 is free parameter. Also, we let s 2i (y) be the tension spline function interpolating u(x i , y j ) in each segment [y j , y j+1 ] and define where a 2i , b 2i , c 2i and d 2i are unknown coefficients and λ 2 is an arbitrary parameter. To determine explicit expressions for all of coefficients we first denote By using (4) and (5) and from algebraic manipulations we derive where From the continuity of the first derivatives of spline functions s 1j (x) and s 2i (y) at (x i , y j ), we get the following consistency relations: where When λ 1 → 0 that ω 1 → 0,then (α 1 , β 1 ) → (1/6, 1/3) and also when λ 2 → 0, then (α 2 , β 2 ) → (1/6, 1/3) so consistency relations of tension splines defined in (6) and (7) reduce to the following cubic spline relations, respectively:

Spline numerical methods
At the grid point (x i , y j ), Equation (1) may be discretized as We next develop an approximation for Equation (1) in which we use the non-polynomial tension spline functions approximations of second derivatives in xdirection and also in y-direction: From (10) and by neglecting the truncation error, we get Similar to (6), for the (j − 1)-th and (j + 1)-th levels in y-direction we have Similar to (7), for the (i − 1)-th and (i + 1)-th levels in x-direction we have By using Equations (6), (7) and (13)- (19) and after manipulating them we obtain where p 1 = Ak 2 α 2 + Bh 2 α 1 − Cr 1 ,

Truncation error
The partial differential equation (10) may be discretized as where D 2 x = (∂ 2 /∂x 2 ) and D 2 y = (∂ 2 /∂y 2 ). Also we have Now,by using the above notations and expanding both sides of (20) in Taylor series in terms of u(x i , y j ) and its derivatives, we obtain the truncation error of the method as follows: From (24) we conclude that • If we choose we obtain various schemes of order O(h 2 + k 2 ). In particular we can choose α 1 = α 2 = (1/6) and β 1 = β 2 = (1/3), which we denoted this scheme by T 1 . • If we choose we get a new scheme of O(h 4 + k 4 ), which we denoted by T 2 .

Convergence analysis
The scheme (20) may be written in the following matrix form: where U is the solution vector in the form U = [u 1,1 , u 2,1 , . . . , u m,1 , . . . , u 1,n , u 2,n , . . . , u m,n ] T and F is the right-hand side vector. P is a square block tridiagonal matrix of order n × n and can be expressed as From relations (21) and (22) we have Then using (25) and (26) we can obtain From (28) we conclude that P is a diagonally dominant matrix and the system (27) has unique solution. For solving the system (27) by P = L+D+R, the block Jacobi method is as follows: where M J , the Jacobi iteration matrix is defined by M J = −D −1 (L + R). Also the block Gauss-Seidel method is For the nine-point scheme (20) the Jacobi iteration formula is The error equation of the iteration is where the values of the error are equal to zero on the boundary of the unit square. Equation (30) can be solved by the method of separation of variables, we let where ξ is the propagation factor.By substituting (31) and (32) in (30), we get and where γ is an arbitrary parameter to be determine and and boundary conditions are V 0 = V m+1 = 0 and W 0 = W n+1 = 0. The solution of difference equation (33) subject to homogeneous boundary conditions is with γ = 2 cos(kπq) and w 1 an arbitrary constant. Also, the solution of (34) with homogeneous boundary conditions is and v 1 is an arbitrary constant. By eliminating l 1 and l 2 from (35) and (38), we get − 2p 3 cos(kπq)).
From (28-(iv)) and (39), we conclude that |ξ | < 1. Therefore, the Jacobi method converges for any initial guess and for any values of h and k. Also the spectral radius ρ of Jacobi and Gauss-Seidel iteration matrices are related by Thus the Gauss-Seidel method converges too.

Numerical results
In this section, we consider the following examples with the known exact solutions. We applied our methods to solve these examples, with various values of h and k.
The computed solutions at grid points are compared with the exact solutions, the maximum absolute errors and convergence orders are tabulated in Tables 1-3.
The orders of convergence are established from the following relation: where u is the exact solution and u 1 and u 2 are computed solutions in two consequence steps (h 1 , k 1 ) and (h 2 , k 2 ). Example 6.1: Consider the equation subjected to Dirichlet boundary condition u = 0 in boundaries of unit square. The exact solution of this problem is u(x, y) = x sin(π y)(1 + cos(π x)).
We applied our methods T 1 and T 2 to solve this problem.
The computed solutions are compared with exact solutions at grid points and the maximum absolute errors of solutions are tabulated in Table 1. The graphs of the exact and computed solutions are given in Figure 1.

Example 6.2: Consider the equation
subjected to Dirichlet boundary condition u = 0 in boundaries of unit square. The exact solution is We applied method T 2 to solve this problem. The maximum absolute errors and orders of convergence are tabulated in Table 2 and are compared with the results in reference [13].
We have solved this problem with methods T 1 and T 2 for two values of c. The maximum absolute errors and orders of convergence of methods are tabulated in Table 3.

Conclusion
We applied non-polynomial tension spline approximation in both space directions for the solution of secondorder elliptic partial differential equations. Using the parameters adopted in cubic tension spline, we have obtained a scheme of order O(h 2 + k 2 ). Also by appropriate choices of the parameters we have increased the order of accuracy to new scheme of order O(h 4 + k 4 ).
The results of the numerical examples show that our methods are accurate and easy applicable.