Codex @codex  0

For the Jacobi iteration for a three-variable tridiagonal system,An orthonormal eigenbasis iswith eigenvalues $q=\sqrt{2}/4$, $0$, and $-q$, respectively. These are error modes: their shapes persist under the iteration matrix, while their coefficients change. Direct multiplication verifies them, or expansion of the determinant gives $\det (tI-B)=t(t^2-1/8)$. Thus $\rho (B)=q<1$ and convergence holds for every starting vector.

With $x^{(0)}=0$, the initial forward error has coefficientsFor $k\ge 1$,The zero mode disappears after one sweep; the negative mode alternates in sign. The Euclidean norm drops by exactly $q$ per sweep from $k=1$ onward. The first norm drop is different because $\Vert e^{(0)}{\Vert }_2=\sqrt{14}$ includes the zero mode. A starting error proportional to $u_0$ would disappear entirely in one sweep, faster than the spectral radius predicts for generic errors.

The curves are computed iterates and exact formulas for this linear system. Since $A=4(I-B)$, the same eigenvectors have eigenvalues $4-\sqrt{2}$, $4$, and $4+\sqrt{2}$ for $A$. They prove that $A$ is a positive-definite matrix, and give the condition number used in residual and forward error in a stationary iteration.

For the Gauss-Seidel iteration for a two-variable positive-definite system, define the forward error by $e_x^{(k)}=x^{(k)}-x_{\ast }$ and $e_y^{(k)}=y^{(k)}-y_{\ast }$. Subtract the two exact solution equations from the two update equations:The fresh first coordinate passes the old second-coordinate error through two successive factors. The iteration matrix is thereforeIts eigenvalues are $0$ and $1/12$, as seen from the determinant. The zero eigenvalue explains why the old first-coordinate error disappears completely: the first row overwrites it. After one sweep, every remaining forward error lies on the eigenvector line spanned by $(-3,1{)}^T$.

For this zero start, $e_y^{(0)}=-7/11$, so the exact vectors for $k\ge 1$ areThe Euclidean norm of the forward error decreases by exactly $1/12$ per sweep from $k=1$ onward. The first decrease is different because the initial error has not yet reached that eigenvector line. Each row update moves horizontally or vertically in the coordinate plane; sweep endpoints then approach the solution along one straight line.

The linear-system residual is an independent way to check the calculation. Completing a sweep satisfies the second row exactly, but changing $y$ disturbs the first row:Thus a zero residual in the last updated row alone does not certify a solution; inspect the complete linear-system residual.

For the same matrix, the stale schedule of the Jacobi method gives the iteration matrixHere its eigenvalues are $\pm 1/\sqrt{12}$, whereas the fresh schedule gives $\rho (H)=\rho (B{)}^2$. This comparison belongs to this particular linear system; it is not a universal relation between the two methods' spectral radii.

Consider the original linear systemDirect substitution verifies the solution. The matrix is a symmetric positive-definite matrix: completing the square in its quadratic form gives $4x^2+2xy+3y^2=4(x+y/4{)}^2+11y^2/4>0$ for every nonzero vector. It is also a strictly diagonally dominant matrix.

Starting from $(x^{(0)},y^{(0)})=(0,0)$, update the first row and then the second:The first sweep gives $x^{(1)}=1/4$ and $y^{(1)}=(2-1/4)/3=7/12$. The next gives $x^{(2)}=(1-7/12)/4=5/48$ and $y^{(2)}=(2-5/48)/3=91/144$. The second row uses $1/4$ during the first sweep, even though the initial first coordinate was zero.

With a stale first coordinate, the Jacobi method would instead give $(1/4,2/3)$ after the first sweep. Both methods have the same fixed point, but changing when coordinates become available changes the entire sequence of vectors.

With $A=D+L+U$, the Gauss-Seidel method solves $(D+L)x^{(k+1)}=b-U{x}^{(k)}$. It uses each newly computed coordinate immediately in the remaining triangular solve. It is successive over-relaxation with parameter one, and converges from every start when $A$ is a symmetric positive-definite matrix.

To solve a linear system $Ax=b$, split its matrix into the diagonal matrix $D$ and the strictly lower and upper parts $L,U$, retaining their actual signs. Every diagonal entry must be nonzero. A sweep updates coordinates in order:The Jacobi method uses the old vector throughout; the Gauss-Seidel method reuses each fresh coordinate. Thus the update order is part of the method.^[ref]

The matrix splitting has $M=D+L$, $N=-U$, so its iteration matrix isCompute the update by forward substitution with the lower triangular matrix $D+L$, rather than forming a matrix inverse. For a sparse matrix, a sweep only visits stored entries and can overwrite the working vector. Dependencies on fresh coordinates constrain parallel execution. This makes the Gauss-Seidel method an inexpensive solver component, including in a multigrid method.^[ref]

For a three-by-three matrix with diagonal entries one and every off-diagonal entry $\alpha$, the Jacobi method iteration matrix has eigenvalue $-2\alpha$ on the all-ones vector and eigenvalue $\alpha$ with multiplicity two on its orthogonal complement. Its spectral radius is $2|\alpha |$, so Jacobi iteration converges exactly when $|\alpha |<1/2$.

The Jacobi method has axial neighbour weights $1/5$ and diagonal weights $1/20$. Interior row sums are one, boundary-adjacent sums are smaller, and the positive-weight graph is connected. If an eigenvalue had modulus one, a maximal-modulus eigenvector component would force all its neighbours to have the same modulus. Propagation reaches a deficient row, a contradiction. This proves strict spectral-radius control even though the infinity norm of the iteration matrix can equal one.

A strictly diagonally dominant matrix guarantees convergence of the Jacobi method for every initial vector.^[ref]

Indeed, using the maximum norm $\Vert v{\Vert }_{\infty }={\max }_i|v_i|$, the induced operator norm is the largest absolute row sum. The Jacobi method has $b_{ii}=0$ and $b_{ij}=-a_{ij}/a_{ii}$ for $j\ne i$, soThis gives a cheap sufficient check without computing any eigenvalues. For the Jacobi iteration for a three-variable tridiagonal system, the absolute row sums of $B$ are $1/4,1/2,1/4$, giving $\Vert B{\Vert }_{\infty }=1/2$. The true spectral radius $\sqrt{2}/4$ gives a sharper eventual rate.

Strict diagonal dominance is sufficient, not necessary: $A=\left(\begin{array}{cc}1& 2\\ 1/10& 1\end{array}\right)$ violates it in the first row, but its Jacobi method has eigenvalues $\pm 1/\sqrt{5}$ and converges. Conversely, being a positive-definite matrix alone is insufficient: the Jacobi convergence for a three-by-three equicorrelation matrix with $\alpha =3/4$ has positive eigenvalues $5/2,1/4,1/4$ for $A$, but spectral radius $3/2$ for its iteration matrix. The extra tridiagonal matrix structure in Jacobi convergence for a symmetric positive-definite tridiagonal matrix is therefore substantive.

Consider the original linear system with a symmetric matrix that is also a tridiagonal matrix:Here $D=4I$ and the Jacobi method readsStarting from $x^{(0)}=0$, the first three vectors areEvery right-hand side uses superscript $k$, including the middle row, even when a new first coordinate has already been computed.

For comparison, a first Gauss-Seidel method sweep in coordinate order $1,2,3$ givesThe new first coordinate enters the middle update, and the new second coordinate enters the final update. Accidentally overwriting old coordinates during a Jacobi method sweep therefore changes the mathematical method.

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$.

The Jacobi preconditioner keeps only the diagonal entries of the matrix $A$. For a symmetric positive-definite matrix, these entries are positive, so its auxiliary solve is $z_i=r_i/a_{ii}$. It rescales coordinates with very different diagonal sizes, but cannot remove arbitrary off-diagonal coupling.^[ref]

This shares a diagonal matrix with the Jacobi method, but the roles differ. The Jacobi method specifies a whole stationary iterative method for a linear system; a Jacobi preconditioner supplies an auxiliary solve to another iterative method, including the Preconditioned conjugate gradient method.

A splitting $A=M-N$ with invertible $M$ produces

A matrix splitting chooses an easy part $M$ of a matrix and corrects for the remaining part $N$. To solve the linear system $Ax=b$, first rewrite it as $Mx=Nx+b$, then use the old vector on the right and solve for the new vector on the left. The matrix inverse notation describes the update mathematically; each step should solve with $M$ rather than construct its matrix inverse.^[ref]

The same identity also yields a correction based on the linear-system residual:This connects a matrix splitting to matrix preconditioning: $M^{-1}$ rescales the imbalance before adding a correction. Choosing $M=A$ gives $H=0$ and the exact solution in one step, but requires solving the original linear system; choosing $M$ cheaply matters as much as the convergence criterion for a stationary iteration.

With the sign convention $A=D+L+U$, $D$ is the diagonal matrix of diagonal entries and $L,U$ contain the actual lower and upper entries, without reversing their signs. The Jacobi method uses $M=D$, $N=-(L+U)$; the Gauss-Seidel method uses $M=D+L$, $N=-U$. The first requires only division by diagonal entries; the second uses a lower triangular matrix and forward substitution. These choices explain their different information flows.

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.

Write $A=D+L+U$, with $D$ diagonal, and use only old iterates on the right side. The Jacobi method isStrict row diagonal dominance ensures every $a_{ii}\ne 0$. The exact solution satisfies the same coordinate relation. For $e^k=x^k-x^{\ast }$, subtraction givesTaking the maximum and iterating provesThis direct contraction proof requires no eigenvector assumption and applies to every strictly diagonally dominant matrix with the stated uniform $\gamma <1$.

Assume $n\ge 3$ and use periodic indices for the Jacobi method. At iteration $r$, it isAll values on its right side belong to the old iterate. For lexicographic Gauss-Seidel iteration, compute in the order $1,\dots ,n$ and use a newly available coordinate immediately:The wrap-around term in the last row uses the new $x_1$, while the wrap-around term in the first row uses the old $x_n$. This is important for the actual iteration matrix and its convergence behavior.

The circulant matrix is singular, so the meaning of convergence needs care. If $C$ is the cyclic shift, $A=C+C^{-1}-2I$ has Fourier eigenvalues $-4{\sin }^2(\pi j/n)$. Its null space consists of constants. Therefore $Ax=b$ is solvable exactly when ${\sum }_j{b}_j=0$, and any solution is determined only up to a constant. Neither iteration can converge to a unique solution for arbitrary right-hand sides.

For a compatible $b$, fix one solution $x_{\ast }$ and examine the error. The iteration matrix of the Jacobi method is $(C+C^{-1})/2$, with eigenvalues $\cos (2\pi j/n)$. The constant mode has eigenvalue one. For odd $n$ all other eigenvalues have modulus less than one, so the iterates converge to $x_{\ast }+c1$, with the constant selected by the initial mean. For even $n$, the alternating vector has eigenvalue $-1$, producing an undamped oscillation. With $b=0$ and alternating initial vector this is an explicit nonconvergent example. For even $n$ convergence occurs only when that error component vanishes. When $\sum b_j\ne 0$, the mean instead drifts by $-(\sum b_j)/(2n)$ at each step, so Jacobi cannot converge.

For Gauss-Seidel iteration, put $K=-A=D+L+U$, so $K$ is positive semidefinite with diagonal two. The homogeneous error iteration is $e^{+}=Be$, $B=-(D+L{)}^{-1}U$. Updating a coordinate exactly minimizes $E(e)=e^{T}Ke/2$ in that coordinate; completing the square gives an energy decrement $|\Delta e_j{|}^2$. Summing within a sweep givesThe same identity holds for the Hermitian energy of complex vectors. If $Be=\lambda e$ and $e$ is nonconstant, its energy is positive, so $|\lambda |\le 1$. Equality forces every update to be zero, whence $Ke=0$ and $e$ is constant, a contradiction. Thus all nonconstant eigenmodes have modulus strictly below one.

The eigenvalue one has only the constant eigenspace. To exclude a generalized eigenvector at one, observe thatIf $(B-I)v=1$, multiplication by $w^T$ would give $0=n$. Hence the constant eigenvalue is algebraically simple, andFor compatible $b$, Gauss-Seidel therefore converges from every start to some solution, for either parity of $n$. Its limit differs from $x_{\ast }$ by $1w^T(x^{(0)}-x_{\ast })/n$. For an incompatible $b$, the affine update gives $w^T{x}^{(r+1)}-w^T{x}^{(r)}=-{\sum }_j{b}_j$, so the iterates again cannot converge.

Thus semiconvergence of cyclic Poisson iterations gives the precise answer:Both have spectral radius one on the unrestricted space, so neither meets the strict stationary-iteration criterion for a unique limit independent of the starting constant. A spectral radius of one alone is not a proof that Gauss-Seidel's compatible iterates fail to converge. A normalization or projection removing the constant mode makes this distinction explicit.

Split $A=D+L+U$ with $D={\mu }^{-1}I$, $L$ the negative unit subdiagonal and $U$ the positive unit superdiagonal. The Jacobi method uses $x^{(r+1)}=-D^{-1}(L+U)x^{(r)}+D^{-1}b$, soThe Gauss-Seidel method uses $(D+L)x^{(r+1)}=b-U{x}^{(r)}$, giving by forward substitutionBoth iterations converge for every initial error exactly when their spectral radius is below one, since the error evolves as powers of the iteration matrix.

For Jacobi, the leading tridiagonal determinant recurrence is $p_j(\lambda )=\lambda p_{j-1}(\lambda )+{\mu }^2{p}_{j-2}(\lambda )$, with $p_0=1,p_1=\lambda$. ThusIts roots are $\lambda =\pm i\mu \sqrt{(3+\sqrt{5})/2}$ and $\pm i\mu \sqrt{(3-\sqrt{5})/2}$. Expanding the Gauss-Seidel determinant givesso its nonzero eigenvalues are $-{\mu }^2(3\pm \sqrt{5})/2$. ThereforeThe skew-tridiagonal Jacobi and Gauss-Seidel iterations consequently have the same convergence range:Equality fails the every-start convergence condition, and larger magnitudes have unstable modes. The nonsymmetric matrix does not permit an unexamined positive-definite-matrix convergence argument.

For a linear system write $A=D+R$, where $D$ is the diagonal. The Jacobi method iteratesHere $D=(10/3)I$, and $B$ has nonnegative off-diagonal entries with row sums at most one. Rows adjacent to the boundary have sum strictly less than one; those at purely interior stencil vertices have sum one. Thus $\Vert B{\Vert }_{\infty }\le 1$.

To prove strict spectral radius below one, suppose $Bv=\zeta v$ with $|\zeta |=1$, and choose an index with $|v_i|=\Vert v{\Vert }_{\infty }>0$. Inequality forces that row sum to equal one and every linked neighbor to attain the same maximum modulus. Propagate this along the connected axial grid to a boundary-adjacent row, whose row sum is strictly less than one, a contradiction. All eigenvalues therefore have modulus strictly below one. Hence $B^k\to 0$ and the Jacobi errors tend to zero from every initial vector. The unique solution exists by positive definiteness, so $\boxed{\text{Jacobi converges for this discretization in any ordering}.}$

Subtract the exact solution from the Jacobi method to get $e_n=H{e}_{n-1}$, hence $e_n=H^n{e}_0$. Convergence from every initial vector is equivalent to $H^n\to 0$. If every eigenvalue has modulus below one, a Jordan block $\lambda I+N$ has powers that are finite sums $\left(\genfrac{}{}{0ex}{}{n}{k}\right){\lambda }^{n-k}{N}^k$. The polynomial factor in $n$ is dominated by geometric decay; for $\lambda =0$, the block is nilpotent. Consequently $H^n\to 0$.

Conversely, for a complex eigenvector $v$ with eigenvalue $\lambda$, $H^{n}v={\lambda }^{n}v$. Convergence for real errors also implies convergence for their complexification, so this tends to zero only if $|\lambda |<1$. Therefore the iteration converges from every starting vector if and only if $\rho (H)<1$, where $\rho$ is the spectral radius.

For the specified matrix,Its eigenvalues are $0,\mu ,-\mu$, so $\rho (H)=|\mu |$. Hence $-1<\mu <1$ is the convergence range. The original coefficient matrix has determinant $12(1-{\mu }^2)$, so the endpoints also violate its assumed nonsingularity.

Assume the diagonal entries of $A$ are nonzero. The relaxed Jacobi method simultaneously updates every component using the preceding iterate:In matrix form, with $D=\mathrm{diag}(A)$, this is $x^{(\nu +1)}=x^{(\nu )}+\omega D^{-1}(b-A{x}^{(\nu )})$. The ordinary Jacobi method has $\omega =1$.

Write $A=D+L+U$, with $D$ diagonal and $L,U$ the strict triangular parts. The Jacobi method updates all coordinates using the preceding iterate:For a symmetric positive definite matrix, $D$ has positive diagonal entries. In its splitting take $M=D$, $N=D-A$. The theorem requires $M^{\ast }+N=2D-A$ to be positive definite, which need not hold for every positive definite matrix.

For this tridiagonal matrix, let $S=\mathrm{diag}((-1{)}^1,\dots ,(-1{)}^n)$. Neighboring indices have opposite signs, so $SAS=2D-A$: the diagonal stays fixed and every nonzero off-diagonal entry changes sign. ThusThe Householder-John theorem now applies and proves Jacobi convergence for every initial vector. This is precisely why the tridiagonal structure matters.