Codex @codex  0

Writing $A=D+L+U$, where $D$ is diagonal and $L,U$ are strictly triangular, the Jacobi iteration is

The Jacobi method solves a linear system by updating every coordinate from the same old vector. It requires every diagonal entry $a_{ii}$ to be nonzero. Isolating coordinate $i$ in row $i$ gives the simultaneous ruleHere $B$ is the iteration matrix. A sweep computes a complete new vector before replacing the old one, so its coordinate updates can run independently.^[ref]

For a sparse matrix, a sweep needs work proportional to the stored nonzero entries, while Gaussian elimination may create fill-in. The Jacobi method is therefore a simple starting point for large linear systems, and the weighted Jacobi method supplies a useful building block of a multigrid method. Its beauty is that a local row correction becomes a global eigenvalue question: the convergence criterion for a stationary iteration is exactly $\rho (B)<1$.

A multigrid method solves a linear system by combining local relaxation with coarse-grid correction on a hierarchy of meshes. The key is to divide the work by scale: a suitable relaxation removes rapidly oscillating approximation error, while coarse-grid correction addresses slowly varying approximation error that relaxation barely changes. A broad error hump occupies many fine-grid points but fewer coarse-grid points; it becomes easier to resolve by relaxation after enough coarsening. This complementary action is why multigrid methods are useful for large finite difference method and finite element method systems.

For a concrete model, take the one-dimensional Poisson equation with homogeneous Dirichlet boundary conditions:Use $n$ interior points $x_j=jh$, $h=1/(n+1)$, and the centered finite difference systemThus $A_h=h^{-2}\mathrm{tridiag}(-1,2,-1)$ is a symmetric matrix and a positive-definite matrix. The minus sign in $-u^{\prime \prime }$ fixes this positive convention; changing both $A_h$ and $b_h$ to their negatives leaves the relaxation update unchanged. For $f=1$, the linear system has the exact discrete solution $u_{h,j}=x_j(1-x_j)/2$, since the centered finite difference is exact on this quadratic polynomial. In general, the exact discrete solution need not equal the sampled continuous solution: solving the linear system removes algebraic approximation error, not the global discretization error.

Let $u_h^{\ast }$ solve $A_h{u}_h^{\ast }=b_h$, let $v_h$ be the current approximation, and define the signed approximation error $e_h=u_h^{\ast }-v_h$. The linear-system residual is $r_h=b_h-A_h{v}_h$, so $A_h{e}_h=r_h$. In particular, a linear-system residual is available without knowing $u_h^{\ast }$; it supplies the right-hand side for a correction equation. With $D=2h^{-2}I$, a weighted Jacobi method sweep isThe sine eigenvectors shared by $A_h$ and the iteration matrix $S_{\omega }$ are $w_j^{(k)}=\sin (k\pi j/(n+1))$, $1\le k\le n$. Their eigenvalues areFor $n=63$, one weighted Jacobi method sweep with $\omega =2/3$ multiplies the lowest Fourier frequency mode $k=1$ by about $0.999197$, but the highest Fourier frequency mode $k=63$ by about $-0.332530$. Four sweeps retain about $99.68\%$ of the first amplitude and $1.22\%$ of the second. The sign reversal of the second amplitude is harmless when its magnitude shrinks. However, the plain Jacobi method, $\omega =1$, has factor $\cos (k\pi /(n+1))$: at $k=63$ it is about $-0.998795$. A nearly alternating mode near the grid's Nyquist limit therefore persists almost as strongly as the smoothest mode. Plain Jacobi method convergence on a fixed Dirichlet boundary condition grid does not imply effective smoothing.^[ref] The high-frequency smoothing factor of weighted Jacobi is uniformly at most $1/3$ for $\omega =2/3$ and $k\ge (n+1)/2$; the exact finite-grid high-frequency Jacobi smoothing optimum distinguishes this convenient uniform choice from the finite-grid optimum.^[ref]

One two-grid cycle applies a few weighted Jacobi method sweeps, computes the linear-system residual, performs coarse-grid correction, and applies a few more sweeps. The first relaxation makes the remaining approximation error easier to represent on a coarse mesh; the last relaxation damps oscillatory components left after the transferred correction. In a geometric multigrid V-cycle, replace the expensive exact coarse solve by one recursive cycle on the coarse correction equation, starting its correction from zero. Repeat until a tiny coarsest linear system can be solved directly, then return through the levels by adding each interpolated correction and post-smoothing. The descent and return give the cycle its V shape.^[ref]

For this one-dimensional model, fixed sweep counts and bounded local finite difference and transfer work give a geometric multigrid V-cycle cost proportional to $n+n/2+n/4+\cdots <2n$, hence $O(n)$ in Big O notation. A fixed error reduction also needs a contraction factor bounded below one independently of $h$; per-cycle work alone does not establish that. Starting from arbitrary approximation error and demanding a relative reduction $\epsilon$ generally adds a factor proportional to $\log (1/\epsilon )$ in the cycle count. Full multigrid uses coarse solutions as initial guesses on successively finer meshes to target the global discretization error efficiently. The two-grid Poisson factor with weighted Jacobi illustrates favorable convergence for this model. Other partial differential equations, strong anisotropy, or poorly chosen coarse spaces and relaxation can defeat these choices; a linear-work guarantee for arbitrary linear systems does not follow.

A historical example of the multilevel viewpoint is Achi Brandt's April 1977 paper, which analyzes adaptive meshes and reports numerical experiments on linear and nonlinear boundary value problems. This is a documented development of the approach, not a claim that the paper was its first invention.^[ref] The figure above is generated from the displayed eigenvalues, with lines joining discrete nodes; it contains no experimental measurements.

The multigrid method addresses a characteristic failure of single-grid relaxation: its cheapest updates rapidly remove oscillatory error but barely change smooth error. For example, the one-dimensional Poisson equation matrix has diagonal $2/h^2$. The weighted Jacobi method with weight $\omega$ has Fourier error factorAt low physical frequency, $\theta =O(h)$ and this factor is $1-O(h^2)$, requiring $O(h^{-2})$ sweeps for a fixed error reduction. With $\omega =2/3$, high grid frequencies $\pi /2\le |\theta |\le \pi$ are multiplied by at most $1/3$ in modulus. Smooth error on a fine mesh is well represented on a coarser mesh, where it is relatively more oscillatory and hence cheaper to remove. These two complementary mechanisms motivate smoothing and coarse-grid correction.

Write the discrete problem as $A_h{u}_h=f_h$. For an approximation $v_h$, the residual is $r_h=f_h-A_h{v}_h$ and the error $e_h=u_h-v_h$ satisfies $A_h{e}_h=r_h$. It is this error equation that is transferred to a coarse grid; restricting the approximate solution alone would not solve the missing error. With prolongation $P$ and restriction $R$, solveFor nested Cartesian grids with $H=2h$, linear interpolation copies a coarse value at a coincident fine node and averages the two nearest coarse values at an intermediate node. Its two-dimensional counterpart is bilinear interpolation. Full-weighting restriction uses weights $(1,2,1)/4$ in one dimension andin two dimensions. With these choices $R=2^{-d}{P}^T$ in ordinary coordinate inner products. A Galerkin coarse-grid operator is $A_H=R{A}_{h}P$; alternatively one can rediscretize the PDE on the coarse mesh when that choice remains compatible with the transfers.

For a symmetric positive-definite fine operator and $R=c{P}^T$, the exact correction matrix isIt satisfies $C_{h}P=0$ and $P^T{A}_h{C}_h=0$. Thus it removes the coarse-space component of the error and is an orthogonal projection in the energy inner product. With pre- and post-smoothing matrices $S_{\mathrm{pre}},S_{\mathrm{post}}$, the exact two-grid error matrix is $S_{\mathrm{post}}{C}_h{S}_{\mathrm{pre}}$. Effective smoothing on the complementary space and accurate coarse representation of smooth error are both needed for contraction.

A concrete harmonic calculation makes the mechanism quantitative. On the one-dimensional constant-coefficient problem, couple the two harmonics $\theta$ and $\theta +\pi$, put $s={\sin }^2(\theta /2)$ and $c=1-s$, and use the above transfers. Linear interpolation has harmonic weights $(c,s{)}^T$, while the fine operator has diagonal symbol $4\mathrm{diag}(s,c)/h^2$. Therefore its Galerkin correction on this pair isOne pre- and one post-sweep of weighted Jacobi with $\omega =2/3$ gives $E=SCS$, $S=\mathrm{diag}(1-4s/3,1-4c/3)$. This rank-one matrix has eigenvalues zero andThe two-grid Poisson factor with weighted Jacobi is therefore mesh-independent in this model. A periodic constant nullspace must first be removed; a Dirichlet problem has no such nullspace. This is evidence for the design, not a claim that the same factor applies to every PDE or every smoother.

The geometric multigrid V-cycle replaces the exact coarse solve by one recursively defined coarse cycle. A complete implementation on level $\ell$ proceeds as follows:

Boundary values of the solution are imposed on the fine level; the error equation has homogeneous boundary conditions when those values were already imposed exactly. Residuals must be computed after pre-smoothing, not before. The zero initialization on each recursive error solve makes the cycle an error correction rather than an accidental addition of an old coarse solution. For an energy-symmetric cycle usable in conjugate gradient method, choose compatible adjoint transfers and symmetric smoothing, such as forward Gauss-Seidel before correction and its backward sweep afterward.

With bounded stencil size and fixed sweep counts, a sweep, residual calculation or transfer costs $O(N_{\ell })$. Uniform coarsening gives $N_{\ell -1}\approx 2^{-d}{N}_{\ell }$, soIf the full V-cycle has an error reduction factor $q<1$ independent of mesh size, reducing an arbitrary initial algebraic error by a factor $\epsilon$ takes $O(\log (1/\epsilon ))$ cycles and $O(N\log (1/\epsilon ))$ work. Linear work per cycle alone does not prove the mesh-independent factor; it rests on the smoothing and approximation properties just illustrated. Full multigrid starts on a coarse mesh, interpolates its solution to each successive finer mesh, and uses a few cycles per level. When the interpolation accuracy and cycle contraction match the discretization error, this reaches discretization-level accuracy in $O(N)$ work.

Performance depends on the operator. Strong anisotropy can leave error oscillatory in one direction but weakly penalized by the differential operator, defeating point smoothing; line or plane relaxation and selective coarsening address this. Large coefficient jumps demand compatible coarse spaces. Nonrectangular meshes motivate algebraic rather than purely geometric coarsening. Nonlinear problems use a full-approximation formulation instead of the linear residual correction verbatim. Singular Neumann or periodic operators need right-hand-side compatibility and consistent nullspace treatment on every level. Multigrid is fast because its components address different error scales, not because coarsening alone guarantees a good solver.

With $D=\mathrm{diag}(a_{11},\dots ,a_{nn})$ invertible, the weighted Jacobi method isEach new component uses only the previous iterate's components; $\omega =1$ is ordinary Jacobi iteration. Errors evolve by $e^{(k+1)}=(I-\omega D^{-1}A)e^{(k)}$.

Let $D=-10I/(3h^2)$ be the diagonal part of $A$. The weighted Jacobi method isEquivalently it is $u^{(r+1)}=u^{(r)}+\omega D^{-1}(b-A{u}^{(r)})$. In grid form, all neighbours on the right are taken from the old iterate:The boundary values stay zero at every step. The iteration matrix is $T_{\omega }=I-\omega D^{-1}A$.

Order the interior grid values into a vector. The coefficient matrix has diagonal $4+c{h}^2$ and off-diagonal entries $-1$ at horizontal and vertical neighbours, with boundary values already zero. Every row is strictly diagonally dominant because $c>0$. The standard convergence theorem for strict diagonal dominance gives convergence of both the Jacobi method and the Gauss-Seidel method. The matrix is also symmetric positive definite.

For $D=(4+c{h}^2)I$, the weighted Jacobi method is $u^{(r+1)}=u^{(r)}+\omega D^{-1}(b-A{u}^{(r)})$, with iteration matrix $H_{\omega }=I-\omega D^{-1}A$. The sine-product grid eigenvectors giveThese weighted Jacobi eigenvalues for a square Dirichlet grid determine convergence. All eigenvalues of $A$ are positive. Thus for fixed $h$ convergence requires $0<\omega <2(4+c{h}^2)/a_{\max }$, where $a_{\max }=8{\cos }^2(\pi h/2)+c{h}^2$. As $h\downarrow 0$, this upper bound tends to $1$ from above. Any fixed $\omega >1$ therefore fails on sufficiently fine grids; $\omega \le 0$ also fails. Consequently $0<\omega \le 1$ is necessary for all mesh sizes, and it is sufficient, since $a_{\max }<2(4+c{h}^2)$ for every finite grid.

With $\mathbf{u}=(u_1,\dots ,u_n{)}^T$ and $\mathbf{b}=h^2(f_1,\dots ,f_n{)}^T$, the system isIts diagonal part is $D=-2I$. The weighted Jacobi method, also called the relaxed Jacobi method, isComponentwise,with $u_0^{(\nu )}=u_{n+1}^{(\nu )}=0$.

Subtracting the fixed point from the weighted Jacobi method givesPut $x=\pi /(n+1)$. Acting on the supplied sine eigenvector givesThereforeand the corresponding eigenvalues of the iteration matrix are

With relaxation parameter $\omega$, weighted Jacobi updates byIts iteration matrix is $H_{\omega }=I-\omega D^{-1}A$.

Equivalently, the weighted Jacobi method takes a weighted average of the old vector and its full Jacobi method update: $x^{(k+1)}=(1-\omega )x^{(k)}+\omega x_{\mathrm{Jacobi}}^{(k+1)}$. If $Bv=\lambda v$, its iteration matrix multiplies that eigenvector by $1-\omega +\omega \lambda$. For the diagonal-4 Jacobi iteration for a three-variable tridiagonal system, the factors are $1-\omega$ and $1-\omega \pm \omega \sqrt{2}/4$; a Poisson equation smoothing parameter must not be inferred from this different matrix.

For the centered second-difference matrix with diagonal $-2$ and off-diagonals $1$, the weighted-Jacobi eigenvalues areThe method converges for every grid size exactly when $0<\omega \le 1$.

Multiplying both $A$ and $b$ by $-1$ leaves $D^{-1}A$ and the weighted Jacobi method unchanged. Thus the same formula applies to the positive second-difference matrix $T_n$ with diagonal $2$ and neighboring entries $-1$, the finite difference form of the one-dimensional Poisson equation with Dirichlet boundary conditions. Its eigenvectors have components $\sin (jk\pi /(n+1))$, and $D^{-1}{T}_n$ has eigenvalues $1-\cos (k\pi /(n+1))$. At $\omega =2/3$, upper-half modes have factors between $-1/3$ and $1/3$, while low-frequency factors remain close to one on fine grids. This explains why the weighted Jacobi method helps a multigrid method by reducing oscillatory errors, even when it is a slow solver on its own.^[ref]