Jacobi method
The Jacobi method is an iterative algorithm in numerical linear algebra used to approximate solutions to systems of linear equations of the form Ax = b, where A is an n \times n matrix with non-zero diagonal entries, x and b are vectors in \mathbb{R}^n, and the method refines an initial guess through successive updates without requiring matrix inversion.[1] It operates by decomposing A into its diagonal part D and the remaining off-diagonal part N = L + U (where L and U are the strict lower and upper triangular portions), yielding the iteration x^{(k+1)} = D^{-1}(b - N x^{(k)}), with convergence to the unique solution occurring under conditions such as the spectral radius of the iteration matrix D^{-1}N being less than 1.[2]
Named after the German mathematician Carl Gustav Jacob Jacobi, the method was introduced in 1845 in his paper "Über eine neue Auflösungsart der bei der Methode der kleinsten Quadrate vorkommenden linearen Gleichungen," published in Astronomische Nachrichten, where it was developed as a relaxation-type technique for solving normal equations in least squares problems arising from astronomical calculations.[3] Jacobi's approach built on earlier ideas from Carl Friedrich Gauss's 1823 work on iterative methods but emphasized simultaneous updates to all variables in each iteration, contrasting with sequential updates in later variants like the Gauss-Seidel method.[4] The method guarantees convergence for any initial guess if A is strictly diagonally dominant, meaning |a_{ii}| > \sum_{j \neq i} |a_{ij}| for all i, and it provides an error bound \|e^{(k)}\| \leq \rho(D^{-1}N)^k \|e^{(0)}\|, where e^{(k)} is the error at iteration k and \rho denotes the spectral radius.[1]
As the simplest stationary iterative method, the Jacobi method is highly parallelizable since component updates depend only on the previous iteration's values, making it suitable for distributed computing environments, though it often converges more slowly than the Gauss-Seidel method due to delayed incorporation of updates.[2] It serves as a foundational technique in numerical analysis, with applications in solving discretized partial differential equations (such as Poisson's equation via finite differences), eigenvalue problems, and large-scale systems in fields like chemical engineering and physics, where matrices exhibit diagonal dominance.[5] Despite its limitations for ill-conditioned or non-dominant systems, the method's elegance and ease of implementation have influenced subsequent iterative solvers, including successive over-relaxation (SOR) and preconditioned variants.[6]
Matrix-based iteration
The Jacobi method is an iterative technique for approximating the solution to a system of linear equations Ax = b, where A is an n \times n nonsingular matrix with nonvanishing diagonal entries, and x, b \in \mathbb{R}^n.[1][7]
The method begins by decomposing A into its diagonal matrix D, whose entries are the diagonal elements of A, and the remainder R = A - D.[1] This splitting yields A = D + R, and assuming D is invertible, the Jacobi iteration is given by
x^{(k+1)} = D^{-1} (b - R x^{(k)})
for k = 0, 1, 2, \dots.[1][7]
Rearranging the iteration formula leads to the fixed-point form x = Gx + c, where the iteration matrix G = -D^{-1} R and the constant vector c = D^{-1} b.[1] Often, R is expressed as the sum of the strict lower triangular part L and strict upper triangular part U of A, so R = L + U and G = -D^{-1}(L + U).[7]
An initial approximation x^{(0)} is selected, typically the zero vector or an informed estimate from the problem context, and subsequent iterates x^{(k+1)} are computed until a convergence criterion is met, such as \|x^{(k+1)} - x^{(k)}\| < \epsilon for a small tolerance \epsilon > 0.[1]
The error at each step, defined as e^{(k)} = x^{(k)} - x^* where x^* is the exact solution, propagates according to the matrix equation
e^{(k+1)} = G e^{(k)}.
This relation highlights how the spectral properties of G influence the reduction of error over iterations.[1] This matrix formulation provides a compact theoretical view equivalent to component-wise updates for each equation.[7]
Component-wise iteration
In the component-wise formulation of the Jacobi method, the solution vector is updated by solving each equation of the linear system Ax = b for a single unknown variable, using the most recent approximations for the other variables from the previous iteration.[1] This approach provides an explicit, scalar-level description of the iterations, which is particularly useful for understanding the computational steps involved in solving the system.[1]
The update rule for the i-th component of the solution vector at iteration k+1 is
x_i^{(k+1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{\substack{j=1 \\ j \neq i}}^n a_{ij} x_j^{(k)} \right),
where a_{ii} \neq 0 for all i = 1, \dots, n.[1] This formula isolates x_i from the i-th equation and substitutes the previous iteration's values for the off-diagonal terms.[1]
A key feature of this iteration is that all components of x^{(k+1)} are computed simultaneously using only the values from x^{(k)}, ensuring that no partial updates from the current iteration influence others within the same step.[1] This simultaneous nature distinguishes the Jacobi method from sequential update schemes and aligns with its matrix splitting interpretation, where the iteration can be compactly expressed as x^{(k+1)} = D^{-1} (b - (L + U) x^{(k)}), with D as the diagonal part of A.[1]
For effective application, the coefficient matrix A is typically assumed to be strictly diagonally dominant, meaning |a_{ii}| > \sum_{j \neq i} |a_{ij}| for each row i, or irreducible with weak diagonal dominance (i.e., |a_{ii}| \geq \sum_{j \neq i} |a_{ij}| for all i, with strict inequality for at least one row).[8] These properties help ensure the iterative process progresses toward the solution without detailed proof of convergence here.[8]
In numerical software, the component-wise updates lend themselves to vectorized implementations, where array operations compute the sums and divisions across all components efficiently, often in environments supporting parallel computation.[9]
The Iterative Process
Algorithm steps
The Jacobi method implements an iterative procedure to approximate the solution of the linear system Ax = b, where A is an n \times n matrix with non-zero diagonal entries and b is the right-hand side vector. The algorithm begins by decomposing A into its diagonal part D (containing the diagonal elements of A) and the remainder R = A - D (off-diagonal elements), enabling the iteration x^{(k+1)} = D^{-1}(b - R x^{(k)}). Equivalently, in component-wise form, each entry of the new iterate is updated as x_i^{(k+1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j \neq i} a_{ij} x_j^{(k)} \right) for i = 1, \dots, n. This update uses values from the previous iteration across all components simultaneously, ensuring that the computation of x^{(k+1)} relies solely on x^{(k)}.
The procedure requires selecting an initial approximation x^{(0)}, often the zero vector for simplicity, though a better-informed guess (e.g., based on domain knowledge) can reduce the number of iterations needed for convergence. The algorithm then proceeds in a loop: compute the next iterate x^{(k+1)}, check a convergence criterion such as the relative change \frac{\|x^{(k+1)} - x^{(k)}\|}{\|x^{(k)}\|} < \epsilon (where \epsilon is a user-specified tolerance, typically on the order of $10^{-6} to $10^{-10}), or monitor the residual \|Ax^{(k+1)} - b\|; if not met, increment the iteration counter and repeat until the criterion holds or a maximum number of iterations (e.g., 1000) is reached to prevent infinite loops.
If any diagonal element a_{ii} = 0, the method fails due to division by zero in the update formula, necessitating matrix pivoting, scaling, or switching to an alternative solver like the Gauss-Seidel method before applying Jacobi. In practice, the algorithm is implemented with safeguards, such as checking for zero diagonals during preprocessing.
Here is a standard pseudocode representation of the Jacobi iteration:
Input: Matrix A (n x n), vector b (n x 1), initial guess x^{(0)} (n x 1), tolerance ε > 0, maximum iterations max_iter
Decompose A into D (diagonal) and R = A - D
If any diagonal entry of D is zero, halt with error (method inapplicable)
k = 0
While k < max_iter:
For i = 1 to n:
x_i^{(k+1)} = (b_i - sum_{j=1 to n, j≠i} a_{ij} x_j^{(k)}) / a_{ii}
If ||x^{(k+1)} - x^{(k)}|| / ||x^{(k)}|| < ε:
Output x^{(k+1)} as approximate solution
Stop
Set x^{(k)} = x^{(k+1)}
k = k + 1
Output x^{(k)} (may not have converged)
Input: Matrix A (n x n), vector b (n x 1), initial guess x^{(0)} (n x 1), tolerance ε > 0, maximum iterations max_iter
Decompose A into D (diagonal) and R = A - D
If any diagonal entry of D is zero, halt with error (method inapplicable)
k = 0
While k < max_iter:
For i = 1 to n:
x_i^{(k+1)} = (b_i - sum_{j=1 to n, j≠i} a_{ij} x_j^{(k)}) / a_{ii}
If ||x^{(k+1)} - x^{(k)}|| / ||x^{(k)}|| < ε:
Output x^{(k+1)} as approximate solution
Stop
Set x^{(k)} = x^{(k+1)}
k = k + 1
Output x^{(k)} (may not have converged)
Implementation considerations
In implementing the Jacobi method for solving linear systems Ax = b, where A is an n \times n matrix, memory requirements include storage for the coefficient matrix A, the right-hand side vector b, and two iteration vectors: the current (old) approximation \mathbf{x}^{(k)} and the updated (new) approximation \mathbf{x}^{(k+1)}. This separation is necessary because all components of the new vector must be computed simultaneously using the old vector to maintain the method's simultaneous update property, preventing premature overwriting that could alter results mid-iteration.[10][11]
The computational complexity of each iteration is O(n^2) arithmetic operations for a dense matrix, as it requires a full matrix-vector multiplication equivalent to computing the off-diagonal contributions for all n components. The total cost depends on the number of iterations needed for convergence, which varies with the matrix properties but can be monitored via the residual norm \|r^{(k)}\|_2 = \|b - A\mathbf{x}^{(k)}\|_2 to terminate when below a tolerance, such as $10^{-6}.[12]
To enhance performance, row scaling can be applied by dividing each row of A and the corresponding entry of b by the diagonal element a_{ii}, normalizing the diagonals to 1 and simplifying the iteration matrix to I - D^{-1}A where D = I. This step simplifies the update formula by eliminating the need for division by a_{ii} in each iteration but does not induce diagonal dominance if the original matrix lacks it; convergence remains guaranteed only if the scaled matrix is strictly diagonally dominant, a property equivalent to that of the original matrix.
The Jacobi method offers strong potential for parallelization, as the update for each component x_i^{(k+1)} depends only on the previous iteration's values and is independent of other components, enabling distribution across processors or threads for the inner loop computations. In distributed-memory settings like MPI, synchronization via all-gather operations after each iteration ensures consistency, with per-iteration communication costs scaling as O(n \log p) for p processors, achieving near-linear speedups for large n (e.g., speedup of 7.55 for n=1000, p=10).[12]
Common pitfalls include division by near-zero or zero diagonal elements, which violates the method's assumption a_{ii} \neq 0 and can cause numerical instability or non-convergence; a preprocessing check for diagonal dominance or pivoting to ensure non-zero diagonals is essential. Debugging tips involve tracking the maximum change in the solution vector \max_i |x_i^{(k+1)} - x_i^{(k)}| alongside residuals to detect stagnation, and using double-precision arithmetic to mitigate floating-point errors in large systems.[13]
Convergence Theory
Conditions for convergence
The Jacobi method converges to the unique solution of the linear system Ax = b if the coefficient matrix A is strictly diagonally dominant, meaning |a_{ii}| > \sum_{j \neq i} |a_{ij}| for all rows i = 1, \dots, n.[14][15] This condition ensures that the spectral radius \rho(G) of the iteration matrix G = D^{-1}(L + U) satisfies \rho(G) < 1, where D, L, and U are the diagonal, strictly lower triangular, and strictly upper triangular parts of A, respectively.[16]
An alternative sufficient condition for convergence is that A is irreducibly diagonally dominant, where A is weakly diagonally dominant (|a_{ii}| \geq \sum_{j \neq i} |a_{ij}| for all i) with at least one strict inequality, and the graph of A is irreducible.[17][18]
Strict diagonal dominance is sufficient but not necessary for convergence; however, its absence can lead to divergence, as demonstrated by counterexamples such as certain non-diagonally dominant matrices where \rho(G) \geq 1, causing the iterates to fail to approach the solution.[16]
In general, the method converges for any initial guess x^{(0)} if and only if \rho(G) < 1; under this condition, \lim_{k \to \infty} x^{(k)} = x^*, the unique solution, with the error satisfying \|e^{(k)}\| \leq \|G\|^k \|e^{(0)}\| in a consistent matrix norm, where the asymptotic rate is determined by \rho(G).[19] These convergence conditions were formalized in the 1950s by Alexander Ostrowski and others through foundational work on the spectral properties of iteration matrices.[20]
Rate of convergence
The rate of convergence for the Jacobi method, assuming the prerequisite condition that the spectral radius ρ(G) of the iteration matrix G is less than 1, is asymptotically linear and governed by ρ(G), where the error norm satisfies ||e^{(k+1)}|| ≤ ρ(G) ||e^{(k)}|| for large k, with e^{(k)} denoting the error vector at iteration k. This implies that the number of iterations required to reduce the error by a factor of 10^{-m} is approximately m / -log_{10}(ρ(G)), making smaller values of ρ(G) essential for practical efficiency.[21][19]
In optimal scenarios, such as when the coefficient matrix A is well-conditioned and strictly diagonally dominant, ρ(G) approaches 0, enabling convergence in very few iterations—often fewer than 10 for systems where the off-diagonal elements are small relative to the diagonal. For symmetric positive definite matrices A, ρ(G) is bounded in terms of the condition number κ(A) = λ_max(A)/λ_min(A), with relations showing that ρ(G) ≤ 1 - c / κ(A) for some constant c depending on the matrix structure, thus linking slower convergence to ill-conditioning.[22][23]
For general nonsymmetric matrices, computing ρ(G) exactly requires finding the eigenvalues of G, which is costly for large systems; instead, practical estimation often employs the power method, iterating a random vector with G until the Rayleigh quotient stabilizes to approximate the dominant eigenvalue magnitude. This numerical approach is crucial for assessing convergence speed in applications, as theoretical bounds may be loose. To enhance the rate beyond the basic Jacobi iteration, polynomial acceleration methods like Chebyshev semi-iteration can be overlaid, minimizing the maximum deviation over a polynomial degree by exploiting eigenvalue distribution estimates, though implementation requires prior approximation of the spectrum.[24][25]
Extensions and Variants
Weighted Jacobi method
The weighted Jacobi method, also known as Jacobi over-relaxation (JOR), generalizes the standard Jacobi iteration by incorporating a relaxation parameter \omega to accelerate convergence for solving the linear system Ax = b, where A is an n \times n matrix and b \in \mathbb{R}^n.[26] The iteration proceeds as
x^{(k+1)} = (1 - \omega) x^{(k)} + \omega D^{-1} (b - R x^{(k)}),
where D is the diagonal part of A and R = A - D is the remainder, or equivalently,
x^{(k+1)} = x^{(k)} + \omega D^{-1} (b - A x^{(k)}).
[27][26] This formulation blends the previous iterate with a weighted update from the Jacobi step, allowing \omega > 1 for over-relaxation or $0 < \omega < 1 for under-relaxation to damp errors more effectively than the unweighted case (\omega = 1).[28]
For symmetric positive definite (SPD) matrices A, convergence is guaranteed when $0 < \omega < 2 / \lambda_{\max}(A_D), where \lambda_{\max}(A_D) is the largest eigenvalue of the scaled matrix A_D = D^{-1/2} A D^{-1/2}, ensuring the spectral radius \rho(G_\omega) < 1 with G_\omega = I - \omega D^{-1} A.[28][27] The method's iteration matrix G_\omega has real eigenvalues in [0, 1) under these conditions, with the error reduction rate determined by \rho(G_\omega).[26]
The optimal \omega minimizes \rho(G_\omega) and, for SPD A, is given by \omega^* = 2 / (\lambda_{\min}(A_D) + \lambda_{\max}(A_D)), yielding a convergence factor of (\kappa(A_D) - 1) / (\kappa(A_D) + 1), where \kappa(A_D) is the condition number of A_D.[28] Computing \omega^* requires estimating the extreme eigenvalues of A_D, often via power iterations, Lanczos methods, or bounds like the Gershgorin circle theorem applied to D^{-1} A.[28][26] This eigenvalue-based optimization significantly improves the asymptotic convergence rate compared to the standard Jacobi method, especially for matrices with clustered spectra.[27]
The weighted variant emerged in the mid-20th century amid developments in relaxation techniques for iterative solvers.[4]
Comparison with Gauss-Seidel method
The Gauss-Seidel method refines the Jacobi iteration by immediately incorporating the newly computed values into subsequent updates within the same iteration cycle.[26] The update rule for its i-th component is
x_i^{(k+1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j < i} a_{ij} x_j^{(k+1)} - \sum_{j > i} a_{ij} x_j^{(k)} \right),
where the sums use the latest available approximations for components j < i.[26]
Unlike the Jacobi method, which relies solely on the previous iteration's vector for all updates and thus supports fully parallel computation of components, the Gauss-Seidel approach is inherently sequential due to these forward dependencies.[26] This sequential structure often yields faster convergence, typically reducing the error by a factor corresponding to roughly twice the rate of Jacobi for compatible systems.[26]
For symmetric positive definite matrices, the spectral radius \rho(H) of the Gauss-Seidel iteration matrix H satisfies \rho(H) \leq [\rho(G)]^2, where G is the Jacobi iteration matrix \rho(G), establishing a theoretical basis for the improved asymptotic convergence rate of Gauss-Seidel.[29]
The Jacobi method's independence between component updates makes it advantageous for parallel architectures, while Gauss-Seidel suits sequential environments, with both requiring analogous conditions like strict diagonal dominance for guaranteed convergence.[30]
Illustrative Examples
2x2 system solution
To illustrate the Jacobi method, consider the following 2×2 linear system, which is strictly diagonally dominant:
\begin{cases}
4x + y = 5 \\
x + 3y = 4
\end{cases}
The exact solution is x = 1, y = 1. The matrix A = \begin{pmatrix} 4 & 1 \\ 1 & 3 \end{pmatrix} satisfies the strict diagonal dominance condition, as |4| > |1| and |3| > |1|, ensuring convergence of the Jacobi iteration.
The component-wise update formulas are:
x^{(k+1)} = \frac{5 - y^{(k)}}{4}, \quad y^{(k+1)} = \frac{4 - x^{(k)}}{3}
Starting with the initial guess \mathbf{x}^{(0)} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}, the first four iterations are computed as follows (values rounded to four decimal places for clarity):
| Iteration k | x^{(k)} | y^{(k)} |
|---|
| 0 | 0.0000 | 0.0000 |
| 1 | 1.2500 | 1.3333 |
| 2 | 0.9167 | 0.9167 |
| 3 | 1.0208 | 1.0278 |
| 4 | 0.9931 | 0.9931 |
| Iteration k | Residual \|\mathbf{A} \mathbf{x}^{(k)} - \mathbf{b}\|_\infty |
|---|
| 0 | 5.0000 |
| 1 | 1.3333 |
| 2 | 0.4165 |
| 3 | 0.1110 |
| 4 | 0.0278 |
3x3 system solution
To illustrate the application of the Jacobi method to a 3×3 system, consider the following diagonally dominant linear system:
\begin{cases}
3x_1 - 0.1x_2 - 0.2x_3 = 7.85 \\
0.1x_1 + 7x_2 - 0.3x_3 = -19.3 \\
0.3x_1 - 0.2x_2 + 10x_3 = 71.4
\end{cases}
[31] The exact solution is \mathbf{x}^* = [3, -2.5, 7]^T.[31] Starting with the initial guess \mathbf{x}^{(0)} = [0, 0, 0]^T, the element-wise updates produce the iterates shown in the table below up to the fourth iteration (rounded to four decimal places for clarity).
| Iteration k | x_1^{(k)} | x_2^{(k)} | x_3^{(k)} |
|---|
| 0 | 0.0000 | 0.0000 | 0.0000 |
| 1 | 2.6167 | -2.7571 | 7.1400 |
| 2 | 3.0008 | -2.4885 | 7.0064 |
| 3 | 3.0008 | -2.4997 | 7.0002 |
| 4 | 3.0000 | -2.5000 | 7.0000 |
Python code example
The Jacobi method can be efficiently implemented in Python using the NumPy library for matrix operations and vectorization, which avoids explicit loops for better performance on larger systems while preserving the iterative essence. The following function includes input validation to ensure the matrix has a non-zero diagonal, convergence monitoring via the Euclidean norm of successive updates, optional verbose logging of iterations, and support for the weighted variant through a relaxation parameter ω (where ω=1 yields the standard method). It returns the approximate solution, number of iterations, and a list of errors for analysis or plotting.
python
import numpy as np
import [matplotlib.pyplot as plt](/page/Matplotlib)
def jacobi(A, b, x0=None, tol=1e-5, max_iter=1000, omega=1.0, verbose=False):
"""
Solves Ax = b using the [Jacobi iterative method](/page/Iterative_method) (or weighted variant).
Parameters:
-----------
A : ndarray
n x n [coefficient matrix](/page/Coefficient_matrix).
b : ndarray
Right-hand side [vector](/page/Vector) of length n.
x0 : ndarray, optional
Initial guess (defaults to zero [vector](/page/Vector)).
tol : float
[Convergence](/page/Convergence) tolerance ([norm](/page/Norm) of update < tol).
max_iter : int
Maximum number of iterations.
omega : float
Relaxation parameter (1.0 for standard Jacobi; 0 < omega < 2 for weighted).
verbose : bool
If True, print iteration progress.
Returns:
--------
x : ndarray
Approximate solution.
iters : int
Number of iterations performed.
errors : list
List of successive errors for [convergence](/page/Convergence) analysis.
"""
n = len(b)
if A.shape != (n, n):
raise ValueError("A must be square and match b's length.")
diag = np.diag(A)
if np.any(diag == 0):
raise ValueError("A must have non-zero diagonal elements for invertibility.")
if x0 is None:
x0 = np.zeros(n)
x = x0.copy().astype(float)
R = A - np.diag(diag) # Off-diagonal matrix
D_inv = 1.0 / diag # Inverse diagonal (vectorized)
errors = []
for k in range(max_iter):
x_old = x.copy()
# Vectorized update: x = omega * D^{-1}(b - R x_old) + (1 - omega) x_old
x = omega * D_inv * (b - R @ x_old) + (1 - omega) * x_old
error = np.linalg.norm(x - x_old)
errors.append(error)
if verbose:
print(f"Iteration {k+1}: error = {error:.6f}")
if error < tol:
if verbose:
print(f"Converged after {k+1} iterations.")
return x, k + 1, errors
if verbose:
print(f"Reached maximum {max_iter} iterations.")
return x, max_iter, errors
import numpy as np
import [matplotlib.pyplot as plt](/page/Matplotlib)
def jacobi(A, b, x0=None, tol=1e-5, max_iter=1000, omega=1.0, verbose=False):
"""
Solves Ax = b using the [Jacobi iterative method](/page/Iterative_method) (or weighted variant).
Parameters:
-----------
A : ndarray
n x n [coefficient matrix](/page/Coefficient_matrix).
b : ndarray
Right-hand side [vector](/page/Vector) of length n.
x0 : ndarray, optional
Initial guess (defaults to zero [vector](/page/Vector)).
tol : float
[Convergence](/page/Convergence) tolerance ([norm](/page/Norm) of update < tol).
max_iter : int
Maximum number of iterations.
omega : float
Relaxation parameter (1.0 for standard Jacobi; 0 < omega < 2 for weighted).
verbose : bool
If True, print iteration progress.
Returns:
--------
x : ndarray
Approximate solution.
iters : int
Number of iterations performed.
errors : list
List of successive errors for [convergence](/page/Convergence) analysis.
"""
n = len(b)
if A.shape != (n, n):
raise ValueError("A must be square and match b's length.")
diag = np.diag(A)
if np.any(diag == 0):
raise ValueError("A must have non-zero diagonal elements for invertibility.")
if x0 is None:
x0 = np.zeros(n)
x = x0.copy().astype(float)
R = A - np.diag(diag) # Off-diagonal matrix
D_inv = 1.0 / diag # Inverse diagonal (vectorized)
errors = []
for k in range(max_iter):
x_old = x.copy()
# Vectorized update: x = omega * D^{-1}(b - R x_old) + (1 - omega) x_old
x = omega * D_inv * (b - R @ x_old) + (1 - omega) * x_old
error = np.linalg.norm(x - x_old)
errors.append(error)
if verbose:
print(f"Iteration {k+1}: error = {error:.6f}")
if error < tol:
if verbose:
print(f"Converged after {k+1} iterations.")
return x, k + 1, errors
if verbose:
print(f"Reached maximum {max_iter} iterations.")
return x, max_iter, errors
For demonstration, consider the following 3×3 diagonally dominant system, a common illustrative example in numerical linear algebra:
A = \begin{pmatrix}
4 & -1 & 0 \\
-1 & 4 & -1 \\
0 & -1 & 4
\end{pmatrix}, \quad
b = \begin{pmatrix}
3 \\
2 \\
3
\end{pmatrix}
The exact solution is x = (1, 1, 1)^T.
The code execution is as follows:
python
A = np.array([[4.0, -1.0, 0.0],
[-1.0, 4.0, -1.0],
[0.0, -1.0, 4.0]])
b = np.array([3.0, 2.0, 3.0])
x, iters, errors = jacobi(A, b, tol=1e-6, verbose=True)
print(f"\nApproximate solution: {x}")
print(f"Iterations: {iters}")
print(f"Final error: {errors[-1]:.2e}")
A = np.array([[4.0, -1.0, 0.0],
[-1.0, 4.0, -1.0],
[0.0, -1.0, 4.0]])
b = np.array([3.0, 2.0, 3.0])
x, iters, errors = jacobi(A, b, tol=1e-6, verbose=True)
print(f"\nApproximate solution: {x}")
print(f"Iterations: {iters}")
print(f"Final error: {errors[-1]:.2e}")
This yields verbose output showing decreasing errors (e.g., starting around 1.17 and dropping below 1e-6 in approximately 14 iterations), an approximate solution of [1.000000, 1.000000, 1.000000], 14 iterations, and a final error on the order of 1e-7, confirming convergence.
To visualize convergence, the errors can be plotted on a semilog scale:
python
plt.figure(figsize=(6, 4))
plt.semilogy(errors, 'b-o')
plt.xlabel('[Iteration](/page/Iteration)')
plt.ylabel('Relative [Error](/page/Error)')
plt.title('Convergence Plot for Jacobi Method')
plt.grid(True)
plt.show()
plt.figure(figsize=(6, 4))
plt.semilogy(errors, 'b-o')
plt.xlabel('[Iteration](/page/Iteration)')
plt.ylabel('Relative [Error](/page/Error)')
plt.title('Convergence Plot for Jacobi Method')
plt.grid(True)
plt.show()
The resulting plot shows a steady logarithmic decay in error, illustrating the method's linear convergence rate for this system.
The weighted extension (ω ≠ 1) can accelerate convergence for compatible matrices, with optimal ω often in (0, 2) determined by spectral analysis, though values near 1 are typical for standard Jacobi applications.[28] For efficiency, the vectorized form in the update step leverages NumPy's optimized linear algebra routines, outperforming element-wise loops for dimensions beyond small examples; the latter can be substituted for pedagogical clarity if needed.