Gauss-Seidel Simulator — Iterative Solution of Linear Systems
Solve a 3-variable linear system Ax=b with the Gauss-Seidel iteration. Sweep SOR omega from 0.5 to 1.9 to compare plain iteration with over-relaxation and watch the log10 residual decay in real time.
With defaults (x_0=0, omega=1, tol=1e-6) the iteration converges in about 11 steps. Setting omega between 1.1 and 1.2 reduces the count further, while omega above 1.5 introduces oscillations and slows convergence; omega greater than or equal to 2 diverges.
Results
—
Solution x
—
Solution y
—
Solution z
—
Iterations
—
Iteration trajectories (x, y, z)
Horizontal axis: iteration k. Vertical axis: component value. Red = x, green = y, blue = z. The dashed line marks the true solution (1, 1, 1). Each iteration is shown as a dot; omega changes both amplitude and convergence speed.
Residual ||x^(k+1)-x^(k)||_inf decay (log10)
Horizontal axis: iteration n. Vertical axis: log10 of the infinity-norm of the increment between successive iterates. Linear convergence appears as a straight line. The red dashed line is the tolerance; steeper slopes mean faster convergence near the optimal omega.
Matrix form (D diagonal, L strictly lower, U strictly upper):
$$(D + \omega L)\,x^{(k+1)} = \omega b - [\omega U + (\omega-1)D]\,x^{(k)}$$
$\omega=1$ recovers Gauss-Seidel, $1<\omega<2$ is over-relaxation and $0<\omega<1$ is under-relaxation. Kahan's theorem requires $\omega \in (0,2)$ for convergence. The optimal value is $\omega_{\text{opt}} = 2/(1+\sqrt{1-\rho_{GS}^{2}})$, where $\rho_{GS}$ is the spectral radius of the Gauss-Seidel iteration matrix.
What the Gauss-Seidel simulator does
🙋
With the defaults the iteration nails (1, 1, 1) in 11 sweeps. What is the advantage of this over just inverting the matrix once?
🎓
Good catch. For a 3-by-3 system, LU factorisation is essentially free, but real CAE problems have a million to a billion unknowns. Direct solvers grow as O(N^3) in both time and memory, which is hopeless for sparse systems coming from finite elements. Gauss-Seidel is the simplest iterative solver: it updates each component using the newest values of the others, so one sweep costs O(non-zeros). On stiffness matrices that are 99% zero it beats direct solvers by orders of magnitude.
🙋
When I push omega to 1.2 the iteration count drops to around 10. Is that what SOR means?
🎓
Exactly. SOR stands for Successive Over-Relaxation. The update is x^(k+1) = (1-omega) x^(k) + omega x_GS^(k+1): you advance slightly further than Gauss-Seidel suggests. For 1
🙋
If I set omega to 1.95 the values start oscillating wildly. Is there an upper bound?
🎓
Yes. Kahan's theorem says SOR can only converge for omega in (0, 2); the bound is sharp. Inside that interval the optimum depends on the matrix. For this tool's A = [[2,1,1],[1,3,1],[1,1,4]] the optimum lives somewhere around 1.1 to 1.3, but for the finite-difference Poisson matrix (CFD pressure correction) omega_opt is more like 1.7 to 1.9. In modern CAE codes Gauss-Seidel rarely runs alone; it appears as the smoother inside multigrid or as a preconditioner for CG and GMRES.
🙋
If I start from x_0 = 10 the iteration count creeps up to about 16. Does the starting guess really matter?
🎓
For diagonally dominant matrices Gauss-Seidel is linearly convergent: the error shrinks by a factor rho < 1 per sweep. The further the initial guess is from the true solution, the more sweeps you need to fall below the tolerance, roughly log(error/tol)/log(1/rho). Linear convergence is slower than Newton's quadratic rate, but Gauss-Seidel pays nothing for a Jacobian, is trivial to code and parallelises well. That is why a few Gauss-Seidel sweeps as a preconditioner inside CG is the workhorse of modern sparse solvers.
Physical model and key equations
The Gauss-Seidel method was discovered independently by Carl Friedrich Gauss (1823, in a private letter) and Philipp Ludwig von Seidel (1874, in print). It is an iterative solver for the linear system Ax=b that updates each component in place, reusing the latest values produced earlier in the same sweep.
The diagonal entries $a_{ii}$ must be non-zero. The tool's matrix A=[[2,1,1],[1,3,1],[1,1,4]] has diagonals 2, 3 and 4 satisfying $|a_{ii}| \geq \sum_{j\neq i} |a_{ij}|$ (2≥2, 3≥2, 4≥2), so it is weakly diagonally dominant; combined with symmetry and positive-definiteness this guarantees SOR convergence for every omega in (0, 2).
Successive over-relaxation (SOR) blends the Gauss-Seidel update with the previous iterate using the relaxation factor omega:
A typical convergence test compares the infinity-norm of the increment between successive iterates against a tolerance: $\|x^{(k+1)} - x^{(k)}\|_\infty < \varepsilon$, with $\varepsilon = 10^{-n}$ in this tool.
Real-world applications
Finite element solvers: Discretising structural, thermal and electromagnetic problems produces large sparse symmetric positive-definite matrices. Commercial codes (Abaqus, ANSYS, COMSOL) use preconditioned conjugate gradient (PCG) or GMRES with SOR and symmetric SOR (SSOR) preconditioners. SSOR is simple, memory-efficient and stable on dominant stiffness matrices, and remains practical up to a few million degrees of freedom.
Multigrid smoothers: CFD pressure-Poisson equations and heat-conduction multigrid solvers use a handful of Gauss-Seidel sweeps at every grid level as a smoother. Red-Black Gauss-Seidel splits the unknowns on a checkerboard so they update independently in parallel, and that pattern dominates GPU multigrid implementations, accelerating the pressure step of flow solvers by orders of magnitude.
Power flow analysis: The non-linear AC power-flow equations are historically solved with a few Gauss-Seidel iterations as a warm start before switching to Newton-Raphson. Gauss-Seidel is less accurate but far easier to implement and rarely diverges, so it survives in contingency screening at electric utilities to this day.
Image reconstruction: The iterative algebraic reconstruction techniques (ART, SART) used in CT and MRI are Gauss-Seidel iterations applied to a huge sparse projection matrix. Combining them with SOR (relaxation parameter typically 0.1 to 0.5) speeds reconstruction and lets clinicians lower the radiation dose for the same image quality.
Common misconceptions and pitfalls
The most common misconception is that Gauss-Seidel is always faster than Jacobi. On diagonally dominant matrices it is roughly twice as fast asymptotically, but on weakly dominant or non-dominant matrices Gauss-Seidel can be slower or even diverge. The tool's A is symmetric positive definite so both methods converge, but in CAE you always check matrix properties before choosing a solver.
A second pitfall is that the optimal omega can be computed. Theory gives omega_opt = 2/(1+sqrt(1-rho^2)), but computing the spectral radius is an eigenvalue problem that costs more than solving Ax=b. In practice one starts with omega in 1.7 to 1.9 (typical for Poisson-like problems), sweeps the value and keeps the one with the fewest iterations. The omega-sweep button reproduces exactly that workflow.
Finally, the test ||x^(k+1) - x^(k)||_inf < epsilon is not the same as accuracy. With small under-relaxation the increment is small even when you are still far from the true solution. A robust stopping rule combines a relative residual ||b - A x||_2 / ||b||_2 with the increment criterion, so production codes always implement both.
FAQ
The Gauss-Seidel method decomposes Ax=b component-wise and updates each x_i in place using the most recent values of the other components. The formula is x_i^(k+1) = (1/a_ii) (b_i - sum_{ji} a_ij x_j^(k)). With the default system A=[[2,1,1],[1,3,1],[1,1,4]], b=[4,5,6], which is diagonally dominant with exact solution (1,1,1), the tool reaches 1e-6 accuracy in about 11 iterations starting from (0,0,0) with omega=1.
The Jacobi method uses only the previous iterate to update every component simultaneously. Gauss-Seidel instead reuses the newest x_j^(k+1) inside the same sweep, which for diagonally dominant matrices typically halves the number of iterations needed. In this tool omega=1 is plain Gauss-Seidel and 1
omega=1 is plain Gauss-Seidel, 0
Sufficient conditions for Gauss-Seidel convergence are (1) strict diagonal dominance |a_ii| > sum_{j!=i} |a_ij| or (2) symmetric positive definiteness of A. Zero or tiny diagonal entries cause pivoting trouble. For SOR, Kahan's theorem guarantees divergence if omega is outside (0,2). The tool's A=[[2,1,1],[1,3,1],[1,1,4]] is weakly diagonally dominant and symmetric positive definite, so SOR converges for every omega in (0,2).