5

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:

no converge

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?

Adam
  • 223
  • 1
  • 6

3 Answers3

3

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)

Binxu Wang 王彬旭
  • 2,664
  • 1
    The condition you mention is only a sufficient condition, and this may not be apparent from the way you phrased it. The Jacobi methods converges for any initial approximation if and only if the spectral radius of $D^{-1}(L+U)$ is less than 1, which holds for a class of matrices much wider than strictly diagonally dominant matrices. If you don't impose convergence for all initial approximations, you get an even wider set of matrices. – PierreCarre May 24 '22 at 09:44
2

The iteration matrix $-D^{-1}(L+U)$ has eigenvalues $\pm i \frac{\sqrt{5}}{2}$ and $0$. This means that the iteration function $G(x)=D^{-1}b-D^{-1}(L+U)x$ is not contractive in any norm and the fixed point method (Jacobi's method is just the fixed point method for this particular choice of $G$) cannot be convergent for an arbitrary initial approximation.

PierreCarre
  • 20,974
  • 1
  • 18
  • 34
2

If you really want to use Jacobi Iteration to play with this equation, you can rewrite it as $$\begin{bmatrix}2&-1&1\\2&2&2\\-1&-1&2\end{bmatrix}^T\begin{bmatrix}2&-1&1\\2&2&2\\-1&-1&2\end{bmatrix}\begin{bmatrix}x\\y\\z\end{bmatrix}=\begin{bmatrix}2&-1&1\\2&2&2\\-1&-1&2\end{bmatrix}^T\begin{bmatrix}-1\\4\\-5\end{bmatrix},$$ or as $$\begin{bmatrix}9&3&4\\3&6&1\\4&1&9\end{bmatrix}\begin{bmatrix}x\\y\\z\end{bmatrix}=\begin{bmatrix}11\\14\\-3\end{bmatrix},$$ wich has a diagonally dominant matrix.

You can find related results searching for "\(Ax=b\) jacobi method" on SearchOnMath.

José Claudinei Ferreira
  • 150