Context Navigation

← Previous Version
View Latest Version
Next Version →

Version 35 (modified by Erica Kaminski, 11 years ago) ( diff )
—

Elliptic equations background
Equation Discretization
Matrix form, relaxation form
1. Matrix form (not followed)
2. Relaxation form (followed)
Laplace Equation (source term = 0) Physical Meaning
Poisson Equation Physical Meaning
Boundary Conditions
The code
Tests and Results
Questions
References

Elliptic equations background

Elliptic equations can be thought of as being the steady state limit of the diffusion equation,

u_t = (\kappa u_x)_x + (\kappa u_y)_y + \psi

when both the boundary conditions and the source (or "forcing") term are time in-dependent. Under these conditions, it can be expected that the time dependent terms vanish as the system relaxes to steady state (i.e. in the limit t \rightarrow \infty, \partial t \rightarrow 0), and we are led to the elliptic equation (here in 2D),

(\kappa u_x)_x + (\kappa u_y)_y = f

where u is the dependent variable we are solving for, and f is the forcing term.

The solution to this equation needs to simultaneously 1) satisfy this equation at all points within a bounding region and 2) satisfy the boundary conditions on that region. This then lends an "instantaneous" feel to the system, much different than the wave-like solutions of hyperbolic equations which travel through space with finite speed.

In fact, apart from their classification as either a boundary value problem for the elliptic equation, or initial value problem for the hyperbolic, they can be considered of either a "static" or "time-evolution" nature, respectively. This is a more helpful classification in terms of designing numerical methods to solve these different types of equations (see Fortran Numerical Recipes, Press et al, Vol. 2, Chapt. 19 - Partial Differential Equations). The following figure from that book illustrates this concept:

Some special cases for elliptic equations, of relevance here, occur when \kappa = 1; when f is non-zero, we have the Poisson equation,

\triangledown ^2 u = f

when f is zero, we have Laplace's equation,

\triangledown ^2 u = 0

The solutions to both Poisson's and Laplace's equations satisfy the Uniqueness theorem. That is, if the solution satisfies the aforementioned 2 conditions, it is the unique solution to the equation.

As is, Poisson's equation is set up to solve for f (be it charge density, or matter density, etc.) given some u (e.g. potential- either electrical or gravitational, etc.). However, in most situations we have the opposite information about a system — given the density, we seek the potential. This requires some numerical techniques for solving this 2nd order differential equation, especially for those systems that do not admit closed-form solutions.

Equation Discretization

The 2nd-order discretized version of Poisson's equation in 1D for fixed cell size is,

\frac{1}{h^2}(u_{i+1} + u_{i-1} - 2 u_i) = f

Here, h is the internode size, which for a 1D line of mesh points is h=L/m+1, where L is domain length, and m is the number of grid points.

Matrix form, relaxation form

Generally speaking, discretization leads to systems of equations. For the Poisson equation, there are 2 broad ways of solution. The first is called direct methods, and these relate to solving a large matrix (see next section). The second type of technique are known as relaxation or iteration methods. This was the type I pursued.

Matrix form (not followed)

I did not write a code that uses a direct matrix method for solving the system of equations, but include this section here for completeness.

Given the discretized Poisson equation for a region actually constitutes a system of m unknowns (one equation for each of the m grid points, and boundary conditions are specified for flanking points), we can write the system in matrix form,

\hat{A} \vec{x} = \vec{b}

To get this form, first specify the ordering of unknowns in the vector x. This will fix the ordering of the source terms in the vector b. Then, following the discretized equation you can fill in the matrix A.

Depending on your ordering, your matrix will have a different look. The properties of the system can be studied by studying this matrix. See Numerical Recipes for details.

Relaxation form (followed)

This is the type of method I used to solve Poisson's equation numerically.

Note that the above discretized Poisson equation can be rearranged for u at a given grid point,

u_i = \frac{1}{2}(u_{i-1}+ u_{i+1} - h^2 f_i)

This then can be rewritten as an iteration equation, where the left hand side (LHS) is the iterated, or updated solution, at grid point i, and the RHS is the value of the grid point's neighbors prior to the iteration step.

Thus, starting with initial guesses at the solution for each point on the grid, one can loop through each point on the grid and 'update' it so that it solves this relationship for its neighbors (and source f).

We say that the iteration converges when the LHS-RHS < tol, where tol is some small number for the tolerance.

By evaluating the matrix spoke of earlier, one can prove that this process will both always converge (in some finite number of iterations) and will provides an approximation to how long convergence will take (see Numerical Recipes for nice overview on this process, again chapt 19).

Jacobi iteration (the specific scheme I used), applies this equation to each grid point on the grid, stores each value of ui, and then uses all of the updated ui's in the next iteration step, if convergence was not yet reached. This is different than Gauss-Siedel, which updates half the grid in the first ½ of the iteration step, and the next half in the second (see Leveque). The amount of time to converge is inversely proportional to the grid point spacing. The finer your grid, the longer it takes to converge. A helpful formula is given, here

Laplace Equation (source term = 0) Physical Meaning

Note that the discretization equation with f=0 describes each point as being the average of its neighbors. This is in fact one of the key properties of the harmonic functions, which are the class of solution to Laplace's Equation. Since this is such an ubiquitous equation in physics, and a subset of the Poisson equation, I felt it was important to review some of its defining properties.

The most interesting property of the harmonic function is that it minimizes the distance (in 1D), or the surface area (in 2D), (or the analogous topology in higher dimensions). So in 1D, the solution to Laplace's equation is a line, and in 2D it is a soap bubble (or a waxy film, pulled taught over box with curvy top edges). Thus, solutions to Laplace's equation DOES NOT contain any local mins or maxs.

The equation describes regions that contain no charge (be it electrical, or gravitational), but describes how the potential due to charge elsewhere effects the region of interest.

Poisson Equation Physical Meaning

While the interpretation of Laplace's equation and the harmonic functions are nice, the same does not seem to apply for the more general Poisson's equation. Instead, the solutions to Poisson's equation are called Green's functions, and they are stitched together out of the solutions to the homogenous part of the PDE (namely, Laplace's equation). Here are some references:

constructing Green's functions

related to Poisson's equation

Boundary Conditions

In 1D, either the value of u is to be given at both boundaries (Dirchlet boundary condition), or EITHER a normal derivative (von Neumann) plus the value of u (but not 2 normal derivatives, for this would not satisfy the boundary conditions).

In 2D (or higher), one can have either Dirchlet or von Neumann conditions all around the region.

The code

Here is an html copy of my code and problem.data.

Tests and Results

I followed the suggested exercises here for the 1D problem.

Test 1 - Laplace's Equation

This was just a basic test of the ability of the solver to converge to the correct solution. Given Laplace's equation in 1D is very easy to solve by hand, I could verify the code's ability to reproduce the analytical result. This it did well:

Test 2 - Poisson's Equation with simple forcing

To test the code with a non-zero forcing function, I used the forcing function, f=6x

This also has a closed-form solution which I then used to check the code. It performed correctly.

Here is the forcing function:

Here is the solution:

Test 3 - Poisson's Equation with complicated forcing

Lastly, I used a more complicated forcing function (refer to the reference). It also performed well.

Here is the forcing function:

Here is the solution:

Attached to this page is the .out file for my program which shows the number of iterations until convergence for the different mesh sizes, and the max error of the final iteration.