Some of the formulas on this page are "pre-formatted". You may need to widen your browser window.

                     f(x-h) - 2 f(x) + f(x+h)
          f''(x) =  --------------------------  + O(h^2)
                                h^2

1. Define a grid --- a set of points in the given domain [0, 1].
2. Let U_j be our approximation to the u(x_j), the value of the solution at the grid point x_j.
3. Replace the differential equation by a difference equation at each grid point x_j. Since the underlying differential equation is linear, the equation at each x_j will be linear.
4. If we have n grid points and n unknowns, U_1, ..., U_n, we have an n x n linear system of equations, Ax=b.
5. Solve the linear system using a method that is suitable for the structure of A.

1. Pick a positive integer n, and define the mesh spacing to be h = 1/(n+1).
   For j = 0, 1, 2, ..., n, n+1, let x_j = j * h.
2. Now U_j is our approximation to u(x_j), j = 1, ..., n.
   Note that U_0 = u(o) = 0, and U_(n+1) = u(1) = 0.
   
3. At x_j, we use the finite difference approximation to the second derivative to replace the differential equation by the following difference equation:
   
   

   
                   -U_(j-1) + 2 U_j - U_(j+1)
                  ---------------------------- = f(x_j).
                               h^2
   
   

   Note that for j = 1, the equation involves U_0, which is given by the boundary condition, and is not an unknown.
   Similarly, for j = n, the equation involves U_(n+1), which is 0 in this problem.
4. For n = 6, we get the following linear system in our discretization of the original ODE:

   
                   | 2 -1  0  0  0  0 |   |U_1|          |f(x_1)|
                   |-1  2 -1  0  0  0 |   |U_2|          |f(x_2)|
                   | 0 -1  2 -1  0  0 |   |U_3|          |f(x_3)|
   	        | 0  0 -1  2 -1  0 | * |U_4| = h^2 *  |f(x_4)|.
   	        | 0  0  0 -1  2 -1 |   |U_5|          |f(x_5)|
   	        | 0  0  0  0 -1  2 |   |U_6|          |f(x_6)|
   
   

   The tridiagonal structure of the matrix comes from the finite difference approximation we have used, as do the matrix entries, -1, 2, -1.
5. A tridiagonal system of linear equations can be solved efficiently by a specialized variant of Gaussian elimination in O(n) opertions (sometimes referred to as the Thomas algorithm).

1. We now consider the Poisson equation in two dimensions on the unit square (0,1) x (0,1).
   Let u = u(x,y) be the unknown function, defined on the unit square.
   u_xx is the second partial derivative with respect to x, and u_yy is the same thing, w.r.t. y.
   
   -u_xx - u_yy = f(x,y)
   
   u = 0 on the boundary.
   
2. We lay down a uniform grid on the unit square.
   Pick a positive integer n, and define the mesh spacing to be h = 1/(n+1).
   For j, k = 0, 1, 2, ..., n, n+1, let x_j = j * h and y_k = k * h.
   Then (x_j, y_k), j = 0, 1, ..., n, n+1 and k = 0, 1, ..., n, n+1
   is a uniform mesh on the unit square.
   
3. At (x_j,y_k) we use the finite difference approximation to the second derivative to replace the partial differential equation by the following difference equation:
   
   

   
   [ -U_(j-1,k)+2 U_(j,k)-U_(j+1,k)]+[-U_(j,k-1)+2 U_(j,k)-U_(j,k+1)]
   ----------------------------------------------------------------- = f(x_j,y_k).
                                       h^2
   

   
   Note that for j or k = 0 or n+1, U_(j,k) is a boundary value, and is not an unknown.
   
   If we label the (j,k) point with a P, and its four neighbors by S, W, E, and N, we get a simplified notation for the above equation:

   
                    4 U_P - (U_S + U_W + U_E + U_N)
                   -------------------------------- = f_P
                                    h^2
   
   

4. In two dimensions, we need to introduce a global ordering on the grid points and the unknowns:
   (In the previous equation, uppercase P was a label for the (j,k)-th grid point. Here lowercase p is an integer.)
   
   (j,k) <--> p
   with
   p = (k-1) * n + j
   This is a common choice, and is called natural ordering by rows. Note that as j and k range from 1 through n, p ranges from 1 through N = n*n.
   The relation between the (j,k) indexing and the global p indexing is depicted for the n = 4 case:
   
   With this global ordering, the equation at the p-th grid point can be written as:

   
              4 U_p - U_(p-nx) - U_(p-1) - U_(p+1) - U_(p+nx)
             -------------------------------------------------- = f_p
                                     h^2
   
   

5. The 2D discrete Laplacian matrix has been analysed extensively. For example, its eigenvalues and eigenvectors are known analytically. Numerical analysts around the world refer to this matrix commonly as The Model Problem.
   
   The 2D discrete Laplacian matrix derived from the 4 x 4 grid, is the following 16 x 16 matrix:

   
             |  4 -1  0  0 -1  0  0  0  0  0  0  0  0  0  0  0 |
   	  | -1  4 -1  0  0 -1  0  0  0  0  0  0  0  0  0  0 |
   	  |  0 -1  4 -1  0  0 -1  0  0  0  0  0  0  0  0  0 |
   	  |  0  0 -1  4  0  0  0 -1  0  0  0  0  0  0  0  0 |
   	  | -1  0  0  0  4 -1  0  0 -1  0  0  0  0  0  0  0 |
   	  |  0 -1  0  0 -1  4 -1  0  0 -1  0  0  0  0  0  0 |
   	  |  0  0 -1  0  0 -1  4 -1  0  0 -1  0  0  0  0  0 |
   	  |  0  0  0 -1  0  0 -1  4  0  0  0 -1  0  0  0  0 |
   	  |  0  0  0  0 -1  0  0  0  4 -1  0  0 -1  0  0  0 |
   	  |  0  0  0  0  0 -1  0  0 -1  4 -1  0  0 -1  0  0 |
   	  |  0  0  0  0  0  0 -1  0  0 -1  4 -1  0  0 -1  0 |
   	  |  0  0  0  0  0  0  0 -1  0  0 -1  4  0  0  0 -1 |
   	  |  0  0  0  0  0  0  0  0 -1  0  0  0  4 -1  0  0 |
   	  |  0  0  0  0  0  0  0  0  0 -1  0  0 -1  4 -1  0 |
   	  |  0  0  0  0  0  0  0  0  0  0 -1  0  0 -1  4 -1 |
   	  |  0  0  0  0  0  0  0  0  0  0  0 -1  0  0 -1  4 |
   
   

   
   
   This matrix is generated in Matlab's sparse format by the function sp_laplace that is available here.
   
   In Matlab, try
   
   nx = 4;
   [A,N] = sp_laplace(nx,nx);
   spy(A);
   and you should see the sparsity structure of the discrete laplacian.
   
   

   Fun with sp_laplace

   A theoretical and an empirical comparison of elliptic solvers.

   	Complexity
   Jacobi	N^2 log(N)
   Gauss-Seidel	N^2 log(N)
   SOR	N^(1.5) log(N)
   PCG	N^(1.5)
   Sparse GE	N^(1.5)
   Dense GE	N^3

   
   
   

   More fun with sp_laplace

   Potential surface for a randomly generated charge distribution.
   
   Equi-potentials for this charge distribution.