Once again, we consider solving the equation , where is nonsingular. When we use a factorization method for dense matrices, such as LU, Cholesky, SVD, and QR, we make operations and store matrix entries. As a result, it’s hard to accurately model a real-life problem by adding more variables.
Luckily, when the matrix is sparse (e.g. because we used finite difference method or finite element method), we can use an iterative method to approximate the solution efficiently. Typically, an iterative method incurs operations, if not fewer, and requires storage.
Over the next few posts, we will look at 5 iterative methods:
- Successive Over-Relaxation (SOR)
- Symmetric SOR (SSOR)
- Conjugate Gradient (CG).
In addition, we will solve 2D Poisson’s equation using these methods. We will compare the approximate solutions to the exact to illustrate the accuracy of each method.
1. Poisson’s equation
We consider a partial differential equation (PDE), called Poisson’s equation, in 2D. This way, we can capture any features that we might miss in 1D, but learn iterative methods more easily than in 3D. For simplicity, we set zero Dirichlet boundary condition (BC). However, with additional work, we can also solve Poisson’s equation with nonzero BCs, Neumann BCs, or mixed BCs.
2D Poisson’s equation with zero Dirichlet BC
Given a function defined over the domain , find such that,
Again, for simplicity, we choose a square domain and a nice RHS function:
where are positive integers.
Under these conditions, we have the exact solution:
As an aside, we can generalize the notion of an eigenvalue problem from matrices to functions. We see that is an eigenfunction of the operator corresponding to the eigenvalue . This coincidence will later come into play when we study Conjugate Gradient method.
Poisson’s equation models several engineering problems. The solution can represent the temperature in a steady-state heat flow, the displacement of an elastic membrane under a vertical load, or the electric potential given a charge distribution. The equation is also a good one to study, since the Laplace operator,
shows up often in PDEs. We know its properties well.
That said, it’s difficult to solve Poisson’s equation analytically. One, we do not know the exact solution for a general and . Even if we do, it may be costly to evaluate it at a point , e.g. because is cast as an infinite series. Two, the problem is too strict. It asks for a solution that satisfies the differential equation and BC at infinitely many points in the domain and on its boundary. On top of that, the solution has to be twice-differentiable inside and continuous outside. Talk about demanding! We call our original problem, the strong form.
We bypass these difficulties by solving the problem numerically and approximately. One approach is to satisfy the differential equation and BC at finitely many points. This leads to a finite difference method or a collocation method. Another is to make the residual function equal to 0 on average. This gives us a finite element method. All of these methods will result in a sparse matrix. We call our easier, nicer problem, the weak form.
b. Finite difference method
To solve a PDE with a finite difference method, we need to discretize the domain and the differential operators in the PDE. In our problem, the domain is square and the only differential operator is the Laplace operator, .
First, we discretize the domain with a uniform mesh. Place points inside the domain along each axis, equally spaced.
For brevity, we denote the values that the functions and take at the point by and . In other words,
Often, the mesh size, denoted and defined to be the maximum distance between two neighboring points along any axis, determines how fast our approximate solution can approach the exact when we increase the problem size (i.e. add more points). Since the points are equally spaced, the mesh size is,
The second step is to discretize the differential operator . We must do so, since the continuous version does not make sense in the absence of neighborhoods of points. The most simple approximation is given by the 5-point stencil:
Notice that we approximated the second derivative of at the point by taking a linear combination of the values of at the point and its four immediate neighbors. We can visualize the points and their contributions with a graph:
We have all the ingredients to turn the Poisson’s equation into a matrix equation now, so let’s review what we just learned.
First, we required that the differential equation is satisfied at points (finitely many) in the domain:
Then, we replaced the Laplace operator with a discrete version by using the 5-point stencil:
Afterwards, we apply the zero Dirichlet BC. We can simply set when either of its indices is 0 or . However, if we have nonzero Dirichlet BCs, we must move the known variables to the RHS.
Finally, we gather the unknowns in a natural order (e.g. traverse along the -direction, then ). We can then write a matrix equation,
where the entries of include . The matrix is sparse, because each row has at most 5 nonzero entries.
Note that, when we write a program for an iterative method, it is easier if we store the entries of in a double array (with extra size to handle boundaries). The double array also feels natural, since it mimics the layout of points on the mesh. The same goes for the entries of .
2. Classical iterative methods
Iterative methods such as Jacobi, Gauss-Seidel, SOR, and SSOR rely on decomposing (splitting) the matrix into a sum of matrices. They are so-called classical, and differ by how we split the matrix. On the other hand, Conjugate Gradient relies on creating subspaces, so we will study this method separately.
Given the equation , we can write as a sum of two matrices,
with one condition that is nonsingular.
Since , we see that satisfies the equation,
We call an equation of the form , such as that above, a fixed-point equation. Fixed-point equations are nice because we can solve them iteratively:
The question is, as , will converge to the exact solution ? Better yet, can we stop after finite iterations and be confident that the error is small?
Yes, if the matrix
meets a certain condition.
Thanks to the Contraction Mapping Theorem, we know when the fixed-point method converges to the exact solution.
Theorem 1. Convergence Criterion I
If for some induced matrix norm , then as for any initial guess . The solution exists and is unique.
In fact, we can show that,
i.e. with each iteration, we can reduce the initial absolute error by a factor of .
However, this criterion is not useful in practice, since we do not know which induced matrix norm results in . It is a sufficient condition for convergence, but not a necessary one. If we happen to get with the 1-norm, 2-norm, and -norm, then we cannot conclude that the fixed-point method will not work.
Luckily, thanks to eigenvalues, we have another criterion, one that is a necessary and sufficient condition for convergence. The spectral radius, denoted , is the size of the largest eigenvalue of . The largest eigenvalue is something that we can compute easily. These lemmas motivate using the spectral radius to decide convergence:
For any induced matrix norm , we have .
For any tolerance , there exists an induced matrix norm such that .
Proof for Lemma (i).
Let be the largest eigenvalue of , and a corresponding unit eigenvector. Then,
We see that the spectral radius is a very good estimate of an induced matrix norm. Furthermore, we can say the following:
Theorem 2. Convergence Criterion II
The fixed-point method converges to for any initial guess , if and only if, .
Proof for Theorem 2.
() If , we can find an such that (density of real numbers). Lemma (ii) guarantees the existence of an induced matrix norm such that . Finally, we apply Theorem 1.
() We prove by contrapositive. We assume that and will find an initial guess for which the fixed-point method does not converge.
Select the initial guess to be,
where is the exact solution and is a unit eigenvector corresponding to the largest eigenvalue of , denoted .
Therefore, as ,
Whew. That was probably a lot to read and digest. I suggest that we stop here today and meet each classical iterative method next time. Much like how indicates the rate of convergence, the spectral radius can also. We will analyze the spectral radius for each method and see its effect on convergence for Poisson’s equation.
We can also discretize the operator with a different stencil to study how it affects the matrix equation and the approximate solution . For example, there is a 9-point stencil:
Finite difference method is great for solving PDEs on a simple domain such as a line, a rectangle, or a box, but not so much on a complex domain with slant edges or curves. In general, we use finite element method to solve PDEs. Nonetheless, finite difference method is great for studying iterative methods!