The Jacobi method is the simplest of the iterative methods, and relies on the fact that the matrix is diagonally dominant. Starting from the problem definition:
\[ A\mathbf{x} = \mathbf{b} \]
we decompose $A$ in to $A = L + D + U$, where $L$ is lower triangular, $D$ is diagonal, $U$ is upper triangular.
\[ A\mathbf{x} = L\mathbf{x} + D\mathbf{x} + U\mathbf{x} = \mathbf{b} \]
We then assume that we have an initial guess at the solution $\mathbf{x}^0$, and try to find a new estimate $\mathbf{x}^1$. Assuming that the diagonal $D$ dominates over $L$ and $U$, a sensible choice would be to insert $x^0$ and the unknown $x^1$ into the equation like so:
\[ L\mathbf{x}^0 + D\mathbf{x}^1 + U\mathbf{x}^0 = \mathbf{b} \]
we can rearrange to get an equation for $x^1$. This is easily solved as we can take the inverse of the diagonal matrix by simply inverting each diagonal element individually:
\[ D\mathbf{x}_1 = \mathbf{b} - (L+U)\mathbf{x}_0 \]
Thus we end up with the general Jacobi iteration:
\[ \mathbf{x}_{k+1} = D^{-1}(\mathbf{b} - (L+U)\mathbf{x}_k) \]
The Jacobi method is an example of a relaxation method, where the matrix $A$ is split into a dominant part $M$ (which is easy to solve), and the remainder $N$. That is, $A = M - N$
\[M\mathbf{x}_{k+1} = N\mathbf{x}_k + \mathbf{b}\]
\[\mathbf{x}_{k+1} = M^{-1}N\mathbf{x}_k + M^{-1}\mathbf{b}\]
This can be rearranged in terms of the residual $\mathbf{r}_k = \mathbf{b} - A \mathbf{x}_k$ to the update equation
\[\mathbf{x}_{k+1} = \mathbf{x}_{k} + M^{-1}\mathbf{r}_k\]
For the Jacobi method $M = D$ and $N = -(L + U)$. Other relaxation methods include Gauss-Seidel, where $M = (D + L)$ and $N = -U$, and successive over-relaxation (SOR), where $M = \frac{1}{\omega} D + L$ and $N = -(\frac{\omega - 1}{\omega} D + U)$, where $\omega$ is the relaxation parameter that is within the range $0 \le \omega \le 2$.
For any relaxation method to converge we need $\rho(M^{-1}N) < 1$, where $\rho()$ is the spectral radius of $M^{-1} N$, which is defined as the largest eigenvalue $\lambda$ of a a given matrix $G$:
\[ \rho(G) = \max{|\lambda|: \lambda \in \lambda(G)} \]
For the SOR method, the relaxation parameter $\omega$ is generally chosen to minimise $\rho(M^{-1}N)$, so that the speed of convergence is maximised. In some cases this optimal $\omega$ is known, for example for finite difference discretisation of the Poisson equation. However, in many cases sophisticated eigenvalue analysis is required to determine the optimal $\omega$.
This exercise involves the manipulation and solution of the linear system resulting from
the finite difference solution to Poisson's equation in two dimensions. Let $A$ be a
sparse symmetric positive definite matrix of dimension $(N-1)^2 \times (N-1)^2$ created
using scipy.sparse (for a given $N$) by the function
buildA as follows:
import numpy as np
import scipy.sparse as sp
def buildA(N):
dx = 1 / N
nvar = (N - 1)**2;
e1 = np.ones((nvar), dtype=float);
e2 = np.copy(e1)
e2[::N-1] = 0
e3 = np.copy(e1)
e3[N-2::N-1] = 0
A = sp.spdiags(
(-e1, -e3, 4*e1, -e2, -e1),
(-(N-1), -1, 0, 1, N-1), nvar, nvar
)
A = A / dx**2;
return Aand let $\mathbf{f}_1$ and $\mathbf{f}_2$ be the vectors defined in
buildf1 and buildf2
def buildf1(N):
x = np.arange(0, 1, 1/N).reshape(N, 1)
y = x.T
f = np.dot(np.sin(np.pi*x), np.sin(np.pi*y))
return f[1:,1:].reshape(-1,1)def buildf2(N):
x = np.arange(0, 1, 1/N).reshape(N, 1)
y = x.T
f = np.dot(np.maximum(x,1-x), np.maximum(y,1-y))
return f[1:,1:].reshape(-1, 1)We will consider manipulation of the matrix $A$ and solution of the linear systems $A\mathbf{U}_i=\mathbf{f}_i$. The solution to this linear system corresponds to a finite difference solution to Poisson's equation $-\nabla^2 u = f$ on the unit square with zero Dirichlet boundary conditions where $f$ is either $\sin(\pi x) \sin (\pi y)$ or $\max(x,1-x) \max(y,1-y)$. PDEs of this type occur (usually with some additional reaction and or convection terms) very frequently in mathematical modelling of physiological processes, and even in image analysis.
import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg
import scipy.optimize
import matplotlib.pylab as plt
def jacobi(A, b, x0=None, tol=1e-5, max_iter=1000):
if x0 is None:
x0 = np.zeros_like(b)
x = np.copy(x0)
b_norm = np.linalg.norm(b)
# jacobi method: M = D
M = A.diagonal().reshape(-1, 1)
invM = 1/M
# main relaxation iteration
for i in range(max_iter):
r = b - A @ x
error = np.linalg.norm(r) / b_norm
if error < tol:
break
x += invM * r
return x, i
num = 20
iterations = np.empty((num, 2), dtype=int)
iterations[:] = np.nan
Ns = np.logspace(0.5, 1.5, num=num, dtype=int)
for j, buildf in enumerate((buildf1, buildf2)):
for i, N in enumerate(Ns):
A = buildA(N)
f = buildf(N)
max_iter = 10*N
x, iters = jacobi(A, f, max_iter=max_iter)
if iters < max_iter:
iterations[i, j] = iters
plt.plot(Ns, iterations)
plt.xlabel('N')
plt.ylabel('iterations')
plt.show()scipy.optimize.minimize_scalar.def SOR(A, b, omega, x0=None, tol=1e-5, max_iter=300):
if x0 is None:
x0 = np.zeros_like(b)
x = np.copy(x0)
b_norm = np.linalg.norm(b)
# SOR method
D = sp.spdiags((A.diagonal()), (0), *A.shape)
L = sp.tril(A, k=-1)
M = (1/omega) * D + L
# main relaxation iteration
for i in range(max_iter):
r = b - A @ x
error = np.linalg.norm(r) / b_norm
if error < tol:
break
x += sp.linalg.spsolve_triangular(M, r)
return x, i
N = 64
A = buildA(N)
f = buildf2(N)
def SOR_iterations(omega):
x, i = SOR(A, f, omega, max_iter=10, tol=1e-32)
return np.linalg.norm(A @ x - f)
res = scipy.optimize.minimize_scalar(SOR_iterations, bracket=[0.1, 1.0, 1.99], tol=1e-2)
print('ideal omega is', res.x, 'versus analytic value of', 2 / (1 + np.sin(np.pi/N)))