To solve the system
$$2x_1-\hphantom2x_2+\hphantom2x_3=-1\\2x_1+2x_2+2x_3=\hphantom-4\\-x_1-x_2+2x_3=-5$$
with Jacobi iteration, we let
$$A=2I_3,\qquad L+U=\begin{bmatrix}0&-1&1\\2&0&2\\-1&-1&0\end{bmatrix},\qquad b=\begin{bmatrix}-1\\4\\-5\end{bmatrix}$$
so that $(A+L+U)\cdot x=b$ is our system. Since $A^{-1}=\frac12I_3$, the Jacobi iteration is
$$\text{iter}(x)=A^{-1}(b-(L+U)x)=\begin{bmatrix}\frac{x_2}2-\frac{x_3}2-\frac12\\-x_1-x_3+2\\\frac{x_1}2+\frac{x_2}2-\frac52\end{bmatrix}$$
Certainly the exact solution $(1,2,-1)$ is a fixed point of $\text{iter}$, but when I try to use it it never converges. Here is a plot of the iteration spiraling away from the solution:
The code to make that in Mathematica is
With[{A = 2 IdentityMatrix@3,
LpU = {{0, -1, 1}, {2, 0, 2}, {-1, -1, 0}}, b = {-1, 4, -5}},
With[{iter = Inverse@A.(b - LpU.#) &, init = {10, -10, -1}},
Graphics3D[{Arrow@NestList[iter, init, 20], Point@init},
Boxed -> False]]]
Have I done something wrong?
Best Answer
Your implementation is not wrong! But Jacobi algorithm is guaranteed to converge only for strictly diagonally dominant system of linear equations. wiki
Your system is not strictly dominant. So it's not too surprising that it diverges...
Though compare this question I answered a few days ago, where Gauss Seidel method is applied to a non-diagonal dominant system, but it converged. (Can non diagonally dominant system of linear equations be solved by jacobi or guass seidel method)