Numerical Solution of a Two Dimensional Poisson Equation with Dirichlet Boundary Conditions

In this paper we have introduced Numerical techniques to solve a two dimensional Poisson equation together with Dirichlet boundary conditions. Specifically two methods are used for the purpose of numerical solution, viz. Finite difference method and Finite element method. The implementation of the solutions is done using Microsoft Office Excel worksheet or spreadsheet. The numerical solutions obtained by these two methods are also compared with each other graphically in two and three dimension.


Introduction
Differential equations arise in many areas of applications such as in Sciences and Engineering. Most of the problems in almost all the scientific and engineering areas are modeled using partial differential equations. Partial differential equations arise when a dependent variable depends on two or more independent variables. Usually analytic solutions for these partial differential equations are difficult to compute and hence alternative methods are used. Numerical methods are widely used to solve partial differential equations within the given domain together with boundary conditions.
In general the domain is divided in to finite number of elements or the domain is made discretization into small triangles or rectangles. The values of the dependent variable at the nodes or grid points of these triangles or rectangles are computed. In literature finite difference and finite element methods are applied to solve two dimensional Laplace equations with Dirichlet boundary conditions using spreadsheet [8]. This motivates us to use spreadsheet implementation for solving two dimensional Poisson equations with Dirichlet boundary condition.
In Section 2, we introduced and discussed the two dimensional Poisson equation with Dirichlet boundary conditions. In Section 3, the finite difference method is applied for solving Poisson equation with Dirichlet boundary conditions. Also, we have made the use of spreadsheet implementation. In Section 4, the finite element method is applied for solving Poisson equation with Dirichlet boundary conditions. Also, we have made the use of spreadsheet implementation. In Section 5, we have compared the results obtained in Sections 3 and 4. The paper ends with concluding remarks in Section 6.

Poisson Equation and the Problem Definition
The Poisson equation is given by [1 -2, 4, 6, and 10] Here the Laplace operator denotes ∆= + . Equation (1) is elliptic second order linear partial differential equation. In (1), the dependent variable , is a function of its arguments and depends on the independent variables and . The function , is known as a source function.For the purpose of numerical solution of the Poisson equation (1) we consider a plane region defined by = , : 0 ≤ ≤ 1 and 0 ≤ ≤ 1 .We will also impose Dirichlet boundary conditions for the dependent variable , on the region defined as 0, = 0, , 0 = 0, The region together with the Dirichlet boundary conditions (2) is illustrated in Figure 1. The -variable increases along the positive horizontal axis while thevariable increases along the positive vertical axis.

Figure 1. Rectangular region R with boundary conditions.
We now split the region R in Figure 1 in to a finite number of rectangular elements. The selection of step lengths 1 4 ⁄ in -direction and ! 1 4 ⁄ in -direction will split the plane region R into 16 rectangular elements.
The nodes and the sides in and on the region R can be classified in to two groups, viz. interior and exterior. The interior nodes and sides lie inside the region R while the exterior nodes and sides lie on the boundary of the region R. Further all the interior nodes and sides are common to adjacent rectangular elements. Of course this fact is exempted for the exterior nodes and sides. The nodes of the region R are numbered and shown in Figure 2. Here our problem is to find numerical values of the function , at the interior node points of the region provided that the Poisson equation (1) and the Dirichlet conditions given in (2) are satisfied.

Finite Difference Method
The finite difference techniques are based on the replacement of differential equations with approximately equivalent finite difference equations. Whenever solving differential equations analytically is not easy then the differential equations are replaced by corresponding finite difference equations and those will be solved. Further, the differential equations provide exact values while the difference equations provide approximate values for the variables. Approximation of the solution relates to the number of grids given region is divided into. More are the grids and better is the approximation.
The finite difference equations are algebraic in nature and their approximate solutions are related to number of grid points of the region .A finite difference solution basically involves the following three steps: (i) Divide the solution region using grids and nodes into small rectangles (ii) Approximate the given differential equation by equivalent finite difference equations that relate the solutions to the grid points. (iii) Solve the difference equations subject to the prescribed boundary and/or initial conditions. To find the solution of the function , on the region , we divide the solution region R of -plane into equal rectangles or meshes of sides and ! along and directions respectively, as shown in Figure 2.We now construct a finite difference equation to represent equation (1) and boundary conditions (2).This will help us to calculate the values of , at the nodes of the rectangles in the region . The detailed derivation of finite difference equation from (1) is given in Appendix 1. Here we directly consider the difference equation as Also, the boundary conditions as Equation (3) represents a difference equation and (4) represents the boundary conditions equivalent respectively to the differential equation (1) and the boundary conditions (2).Here in (3) we have used the notations as ",# & " , # + and ",# & " , # +.Also, the set $ , * , 2 , % , 0 is the set of partition points of the interval [0, 1] along -axis with step length . Similarly, the set $ , * , 2 , % , 0 is the set of partition points of the interval [0, 1] along -axis with step length!. Of course in our present study we choose !.The numerical solutions for , is computed at all the nine free nodes of the region located at the locations , , -: , 2, 3, 4 345-2, 3, 4 . These nine free nodes are numbered as shown in Figure 3.
Gauss Seidel method is applied to solve (3) numerically for the specific value of the function , 1 .The function , in fact admits any constant or expression in terms of and . For simplicity and to illustrate the procedure of numerical solution here we have selected that , 1 . Other expressions for , will be considered as and when it demands in our further work.
At the nodal points as shown in Figure 3 we compute the successive approximate solutions for , using spreadsheet. The iterations and the approximations of the solutions are provided in a tabular form in Table 1.The Gauss -Seidel iterations for every node are continued till the approximate solution values converge to a constant.

Finite Element Method
Finite element methods are essentially methods for finding approximate solutions of partial differential equations defined in a finite region or domain. Finite element method involves the following steps in solving a partial differential equation together with boundary conditions: i. Divide the domain in to a finite number of elements.
ii. Drive the weak formulation corresponding to the given problem. iii. Calculate the stiffness matrix and load vector for each element in the domain. iv. Calculate global stiffness matrix and assemble.
Calculate global load vector and assemble v. Solve the resulting system of algebraic linear equations subject to the satisfaction of the boundary conditions. The given domain is divided into 32 congruent right angled and isosceles triangles. The nodes of the triangles are represented using unenclosed numbers. These unenclosed numbers are called global node numbers. Similarly, the triangles are represented by en-rectangle numbers. These enclosed numbers are called element numbers. The enrectangle number is a number that is enclosed by a rectangle. These details are shown in Figure 4.  (1) is given in appendix 2.
The finite element problem for the approximation of Poisson equation (1) together with boundary condition (2) requires finding C ∈ D C satisfying the following condition: −: = ∇ C • ∇< C 5 = : = < C 5 ∀< C ∈ D C .
Here in (5)  The function ] " is called shape function. The function < C ∈ D C can be expressed through a linear combination of the basis functions of D C in the following way < C , = ∑ < # ] # , , < # = < W \ # 2 #c$ (6) By expressing the discrete solution C in terms of the basis ] " , we have C , = ∑ " ] " , , < " = < W \ " 2 "c$ We now assign numbers 1, 2, and 3 to the nodes of each triangle element. Starting with right angled vertex, number the nodes in anti clockwise direction using 1, 2, and 3 respectively. This numbering procedure will help us in having all elements the same stiffness matrix since these elements have the same geometry. The detailed derivation of equation of stiffness matrix and load vector is given in appendix 3.Here we consider only equitation of stiffness matrix and load vector. The stiffness matrix for each triangular element is calculated using the relation Hereg = 1, 2, 3.
Thus the stiffness matrix for each element is The h, i entry of the global matrix of element T U takes the ,,entry of the stiffness matrix of element T U where h, i = 1,2, … ,25 are global node numbers, and,, -= 1,2,3 are corresponding local node numbers for the global node numbers andO = 1,2, … ,32,. The sum of the 32 global matrices gives the assembled global matrix and the assembled global matrix is given in Table 3. Table 3. Assembled global stiffness matrix. Here j # K L denotes theentry of the load vector j K L .

Global node numbers
Also,] f g = 1, 2, 3 represents a shape function at the three nodes of any triangular element T U where O = 1, 2, … , 32.
Using two dimensional four point rule Gauss quadrature formula, each element has the same load vector The h entry of the global load vector of element T U takes the , entry of the load vector of element T U where h = 1,2, … ,25 are global node numbers, , = 1,2,3 are corresponding local node numbers for the global node numbers andO = 1,2, … ,32,. The sum of the 32 global load vectors gives the assembled global load vector and the assembled global load vector is Applying the boundary conditions (2) the assembled global matrix and the load vector becomes respectively Table 6. Reduced assembled global matrix and load vector. Using Gauss Seidel method successive approximate solution for , at the free nodes of the region is obtained. The results are given in Table 7.

Comparison Between the Two Solution
For the differential equation (1) together with boundary condition equations (2) we have found numerical solutions using two methods, viz. Finite difference and Finite element methods. Here in this section we compare these numerical solutions with each other. The solutions of , at the free nodes of the region obtained using the two methods are given in the Table 8. It can be observed from Table 8 that the maximum difference of the two solutions is 0.0039 at node 7. This difference can be minimized by taking more number of nodal points. The two solutions can be described graphically in two dimensions and three dimensions as shown in Figure 5 and Figure 6 respectively.
As it is seen in Figure 5 the finite element and finite difference solutions are close to each other at the free nodes. As it is seen from Figure 6, the finite difference solution in (a) and finite element solution in (b) are somewhat the same.

Conclusions
In this paper we have introduced spreadsheet implementations of numerical methods for solving Poisson equation in two dimensions with Dirichlet boundary conditions. Other types of boundary conditions and right hand side function of the Poisson equation could be considered to solve the same problem using spreadsheet. Application of finite difference method is somewhat complex in comparison with that of the finite element method. Many physical problems have irregularly shaped boundaries and boundary conditions are expressed using derivatives. Boundary conditions which are expressed using differential equations are difficult to handle using finite -difference techniques because each boundary condition involving a derivative must be approximated by a difference quotient at the grid points. Also, if the shape of the boundary is irregular then placing the grid points on the boundary is difficult. The Finite element method alleviates all difficulties appears in finite difference method. It is independent of boundary conditions of the problem.