Let be a square system of n linear equations, where:
When and are known, and is unknown, we can use the Jacobi method to approximate . The vector denotes our initial guess for (often for ). We denote as the k-th approximation or iteration of , and is the next (or k+1) iteration of .
The computation of requires each element in except itself. Unlike the Gauss–Seidel method, we can't overwrite with , as that value will be needed by the rest of the computation. The minimum amount of storage is two vectors of size n.
Input:initial guess x(0) to the solution, (diagonal dominant) matrix A, right-hand side vector b, convergence criterion
Output:solution when convergence is reachedComments: pseudocode based on the element-based formula above
k = 0while convergence not reached dofori := 1 step until n doσ = 0forj := 1 step until n doifj ≠ ithenσ = σ + aijxj(k)endendxi(k+1) = (bi − σ) / aiiend
increment kend
The standard convergence condition (for any iterative method) is when the spectral radius of the iteration matrix is less than 1:
A sufficient (but not necessary) condition for the method to converge is that the matrix A is strictly or irreducibly diagonally dominant. Strict row diagonal dominance means that for each row, the absolute value of the diagonal term is greater than the sum of absolute values of other terms:
The Jacobi method sometimes converges even if these conditions are not satisfied.
Note that the Jacobi method does not converge for every symmetric positive-definite matrix. For example,
If we choose (0, 0, 0, 0) as the initial approximation, then the first approximate solution is given by
Using the approximations obtained, the iterative procedure is repeated until the desired accuracy has been reached. The following are the approximated solutions after five iterations.
0.6
2.27272
-1.1
1.875
1.04727
1.7159
-0.80522
0.88522
0.93263
2.05330
-1.0493
1.13088
1.01519
1.95369
-0.9681
0.97384
0.98899
2.0114
-1.0102
1.02135
The exact solution of the system is (1, 2, −1, 1).
importnumpyasnpITERATION_LIMIT=1000# initialize the matrixA=np.array([[10.,-1.,2.,0.],[-1.,11.,-1.,3.],[2.,-1.,10.,-1.],[0.0,3.,-1.,8.]])# initialize the RHS vectorb=np.array([6.,25.,-11.,15.])# prints the systemprint("System:")foriinrange(A.shape[0]):row=[f"{A[i,j]}*x{j+1}"forjinrange(A.shape[1])]print(f'{" + ".join(row)} = {b[i]}')print()x=np.zeros_like(b)forit_countinrange(ITERATION_LIMIT):ifit_count!=0:print(f"Iteration {it_count}: {x}")x_new=np.zeros_like(x)foriinrange(A.shape[0]):s1=np.dot(A[i,:i],x[:i])s2=np.dot(A[i,i+1:],x[i+1:])x_new[i]=(b[i]-s1-s2)/A[i,i]ifx_new[i]==x_new[i-1]:breakifnp.allclose(x,x_new,atol=1e-10,rtol=0.):breakx=x_newprint("Solution: ")print(x)error=np.dot(A,x)-bprint("Error:")print(error)