# Gauss-Seidel method

I am trying to solve a system of equations with Gauss Seidel method.I reached this point but there is some problem,and I can't find a result.This is the system.

x1-x2+8x3-x4=1
3x1-x2+2x3-11x4=4
11x1-x2+2x3-2x4=2
-x1+9x2-x3+2x4=1

GaussSeidelMat[a_?MatrixQ, b_?MatrixQ, x0_?MatrixQ, error_Real,
steps_Integer] := Block[
{l, u, x, abs},
x[0] = x0;
l = a SparseArray[{i_, j_} /; j <= i -> 1, {3, 3}];
u = a SparseArray[{i_, j_} /; j > i -> 1, {3, 3}];
Reap[Do[
x[i] = Inverse[l].(b - u.x[i - 1]);
abs = Norm[x[i] - x[i - 1]]/Norm[x[i]];
If[abs < error, Sow@x[i]; Break[]];
If[i == steps, Sow@x[steps]]
, {i, steps}]][[-1, -1, 1]]]

aa = {{1, -1, 8, -1}, {3, -1, 2, -11}, {11, -1, 2, -2}, {-1, 9, -1,
2}};
bb = {{1}, {4}, {2}, {1}};
pp = {{0}, {0}, {0}, {0}};
args = {N@aa, N@bb, N@pp, 0.05, 50};

GaussSeidelMat @@ args

• What is the problem that you reached? What happens when you run the code? Do you get errors? Does it run forever? Have you tried running just one step of your code (rather than 50) to see what it is doing at each step? Try these things first, read the error messages spit out by Mathematica, and see if you can track down the problem, then get back to us if you can't. Feb 21, 2020 at 0:59
• It looks like the dimensions of your arguments (aa, bb, and pp) don't match the arguments of the SparseArrays you are creating as l and u in the code. Fix that, and it seems to run just fine (although it seems like you need to do some normalization or something). Feb 21, 2020 at 1:01
• Side note: Do not use Inverse, better use LinearSolve (which will do the triangular solve for you). Inverse has $O(n^3)$ complexity while the triangle solve is only $O(n^2)$ (and much less if the triangular matrix is sparse) Feb 21, 2020 at 9:41

A necessary and sufficient condition for the convergence of Gauss Seidel iterative method is that the spectral radius of iterative matrix is less than 1. But the absolute value of the eigenvalue of the iterative matrix of this equation is greater than 1, so it is not convergent.

ClearAll["Global*"]
X[0] = {0, 0, 0, 0};(*Initial vector*)
b = {1, 4, 2, 1};
A = {{1, -1, 8, -1}, {3, -1, 2, -11}, {11, -1, 2, -2}, {-1, 9, -1, 2}};
DI = DiagonalMatrix@Diagonal[A];
L = LowerTriangularize[-A, -1];
U = UpperTriangularize[-A, 1];
B = IdentityMatrix[4] - Inverse[DI - L].A;
Abs[Eigenvalues[B]] // N(*The maximum absolute value of the eigenvalue of this matrix is greater than 1, so Gauss Seidel iterative method cannot converge*)


But in the following case, Gauss Seidel iteration can converge:

ClearAll["Global*"]
Solve[{8 x1 - 3 x2 + 2 x3 == 20,
4 x1 + 11 x2 - x3 == 33,
6 x1 + 3 x2 + 12 x3 == 36}, {x1, x2, x3}]

X[0] = {0, 0, 0};(*Initial vector*)
b = {20, 33, 36};
A = {{8, -3, 2}, {4, 11, -1}, {6, 3,
12}};

DI = DiagonalMatrix@Diagonal[A];
L = LowerTriangularize[-A, -1];
U = UpperTriangularize[-A, 1];
B = IdentityMatrix[3] - Inverse[DI - L].A;
N /@ Abs /@ Eigenvalues[B](*It can be seen that the absolute values of the eigenvalues of the iterative matrix B are all less than 1, so the Gauss Seidel iteration converges*)

f = Inverse[DI - L].b;
X[n_ /; 1 <= n] := X[n] = B.X[n - 1] + f
X[20] // N(*The exact solution is {3,2,1}*)


Your code can run with the following modifications, but the result is not convergent:

GaussSeidelMat[a_?MatrixQ, b_?MatrixQ, x0_?MatrixQ, error_Real,
steps_Integer] := Block[{l, u, x, abs}, x[0] = x0;
l = a SparseArray[{i_, j_} /; j <= i -> 1, {4, 4}];
u = a SparseArray[{i_, j_} /; j > i -> 1, {4, 4}];
Reap[Do[x[i] = Inverse[l].(b - u.x[i - 1]);
abs = Norm[x[i] - x[i - 1]]/Norm[x[i]];
If[abs < error, Sow@x[i]; Break[]];
If[i == steps, Sow@x[steps]], {i, 1, steps}]][[-1, -1, 1]]]

aa = {{1, -1, 8, -1}, {3, -1, 2, -11}, {11, -1, 2, -2}, {-1, 9, -1,
2}};
bb = {{1}, {4}, {2}, {1}};
pp = {{0}, {0}, {0}, {0}};
args = {N@aa, N@bb, N@pp, 0.05, 50};

GaussSeidelMat @@ args


So you should consider using other iterative methods to solve this problem.