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 (see Griffith's E&M textbook). 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 potential (e.g. 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

By replacing the 2nd-order derivatives of Poisson's equation with 2nd-order central finite differences, we have the discretized version of Poisson's equation in 1D (See Leveque, Finite Difference Methods for Ordinary and Partial Differential Equations),

\frac{1}{h}(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.

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 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 the Laplace Equation. Since this is such an ubiquitious equation in physics, and a subset of the Poisson equation, I felt it was important to review some of the properties behind it.

• Boundary conditions
• Minimizing distance b/w boundaries (in 1d this is a line, in 2d this is soap bubble)
• Physical interprestaion - no charge within the domain, but elsewhere.

Poisson Equation Physical Meaning

While the interpreation of the solution propertoies are not as clear cut and intuitive as Laplace, the solution can be thoguht of as composing green's functions (cite references).

The code

Explain what my thing does, post it for download

Tests and Results

Test 1 - Laplace's Equation

Test 2 - Poisson's Equation with simple forcing

Test 3 - Poisson's Equation with complicated forcing

Properties - aka only an approximate solution now . Verify 2nd order (how?)

Cite References

References