2.1 Forcing Function 

Now, we are going to examine the first two-dimensional example, we consider the problem

By the Fourier sine series, we may let the exact solution as the form (4.7), and using (4.9) to determine the coefficient  . That is

Now the forcing function is a two-dimensional delta function, using the similar property of (2.9), we have a formula for the exact solution

We are going to use five-point finite difference approximations for this Poisson equation. Now, the square domain is covered by a grid with sides parallel to the coordinate axes and if we divide the interval [0, 1] into n equal parts, letting the grid spacing be h, then the number of internal grid points is  . Moreover, the coordinates of a typical internal grid point are  , and the function value of u at this grid point is denoted by  . For the Laplace's differential operator, let us replace each second derivative with the usual finite difference approximation at the grid point  except the point where the point source is acted,

In order to obtain the difference equation, we need to determine how to replace the forcing function. By using the Green's first identity, we have

where u and v are differentiable functions throughout the domain  . The expression  is simply the dot product of vector between two gradients, and  is normal derivatives, that is,

where  is the outward normal unit vector on the boundary. If we set v = 1 in Green's first identity (2.1.5), we have

Consider the uniform mesh of the problem (2.1.1) and the above figure. For each small mesh we have,

at the region where the point source is acted on. For each boundary of the small square region as shown in the above figure, we have

Using the usual finite difference approximation to replace each derivative, at the nodal point (x, y) = (1/2, 1/2), we have the following difference equation

Summarize the above analysis, when we impose Dirichlet boundary conditions along the square domain and adopt the row by row ordering, the totality of equations at the  internal nodes of the unit square leads to the matrix equation

where u is a vector of order  such that

and A is a matrix of order  given by

with I the identity matrix of order (n - 1) ,  a zero matrix of order (n - 1) and B a square matrix of order (n - 1) given by

and the elements of vector b depend on the boundary values and the forcing function. According to (2.1.10), we have

The Table 6.1 and Table 6.2 show that the solutions at the nodal point (1/4, 1/4) are convergent and the rate of convergence is satisfactory. Although the data which obtained from the numerical results are not so absolute as the results in one-dimensional problems, we also obtain the asymptotic error expansion as

TABLE 6.1

h=1/4	.76374815x10-2	.28454037x10-3	.26556058x10-5	.13540857x10-7
h=1/8	.21227757x10-2	.20273403x10-4	.54823121x10-7	.33494879x10-9
h=1/16	.54589897x10-3	.13184844x10-5	.11863265x10-8	.25417585x10-11
h=1/32	.13746361x10-3	.83517455x10-7	.21038393x10-10	.27730596x10-12
h=1/64	.34428539x10-4	.52395644x10-8	.60169925x10-12
h=1/128	.86110645x10-5	.32803686x10-9
h=1/256	.21530122x10-5

TABLE 6.2

r(1)	.27794184	.071249656	.020644299	.024736160
r(2)	.25716282	.065035178	.021639163	.075884987
r(3)	.25181144	.063343530	.017734067
r(4)	.25045567	.062736160	.028600057
r(5)	.25011414	.062607659
r(6)	.25002857

We also see that extrapolating from the three solutions corresponding to h = 1/4, 1/8 and 1/16, we get a numerical solution with error  which is slightly smaller than the numerical solution error  corresponding to h = 1/256. But the most important thing is that we can get a more accurate solution in a very short time. At other nodal points, we have similar improvement on the numerical solution accuracy by using extrapolation [See Appendix F.6a-F.6c].

Wed Sep 15 10:03:39 HKT 1999