1
$\begingroup$

I am trying to implement the below description from Ch. 11 of Heath's "Scientific Computing An Introductory Survey" of the Gauss-Seidel iterative method for solving a system of linear equations

\begin{align} x^{(k+1)} &= D^{-1}(b - Lx^{(k+1)} - Ux^{(k)})\\ &= (D+L)^{-1}(b - Ux^{(k)})\, \end{align}

However, I must be missing something because writing the equation exactly as is fails to produce the correct results.

using LinearAlgebra: lu
# From dicretization of Laplace equation
A = [
  4.0  -1.0  -1.0   0.0
 -1.0   4.0   0.0  -1.0
 -1.0   0.0   4.0  -1.0
  0.0  -1.0  -1.0   4.0]

# block diagonal matrix of A
D = [
  4.0  -1.0   0.0   0.0
 -1.0   4.0   0.0   0.0
  0.0   0.0   4.0  -1.0
  0.0   0.0  -1.0   4.0]

b = [0, 0, 1, 1] # rhs for laplace w/ Dirichlet boundary conditions
L, U = lu(A)

# gauss seidel using matrix terms
niters = 100
D_plus_L_inv = inv(D + L)
x = zeros(length(b))
for k in 1:niters
  x = D_plus_L_inv*(b - U*x)
end 

# clearly not equal
@show A \ b
# A \ b = [0.12499999999999999, 0.12499999999999999, #0.37499999999999994, 0.37499999999999994]

@show x
# x = [0.021795262674411238, 0.023610129063946845, 0.14893710589872386, 0.14211027447080438]
$\endgroup$
2
  • $\begingroup$ I provide an answer asked on stack overflow here $\endgroup$ Nov 18, 2023 at 14:06
  • 1
    $\begingroup$ That is the same answer that you obtained here and the reason why we suggest to avoid cross-posting. $\endgroup$
    – nicoguaro
    Nov 18, 2023 at 18:30

2 Answers 2

4
$\begingroup$

The $L$ and $U$ of Gauss-Seidel are different from the $L$ and $U$ that come from the LU factorization. For Gauss-Seidel, $L$ and $U$ are what you get if you zero out the upper or lower part, respectively, so that $A=L+D+U$. If you print A-(L+D+U) with your current code, you'll see that the're not equal.

For a pointwise Gauss-Seidel in MATLAB, you want

D = diag(diag(A))
L = tril(A, -1)
U = triu(A, 1)

For a blocwise Gauss-Seidel, I don't know of analogous MATLAB functions, but you can manually copy the block lower, diagonal, and upper parts with a loop. (Or just hard code them like you did D.)

$\endgroup$
2
$\begingroup$

In the Gauss iterative method (Gauss-seidel) we decompose the linear system $Ax=b$ into $x^{k+1} = (L+U)^{-1}(b-Ux^{k})$ as exposed in the reference indicated (Scientific Computing An Introductory Survey Heath, page 469), but as well noted in the previous answer, there was a point in which you used the $LU$ photoration to generate the matrices $L$ and $U$ from the matrix $ A$ and applying them to the Gauss iterative method. These matrices generated by the direct $LU$ method are different from the $L$, $U$ and $D$ of the Gauss iterative method, where.

  • $L_{ij} = A_{ij}$ if $i>j$, and $L_{ij} = 0$ if $i \leq j$.
  • $U_{ij} = A_{ij}$ if $i<j$, and $U_{ij} = 0$ if $i \geq j$.
  • $D$ a diagonal matrix with $D_{ii} = A_{ii}$ for $i=j$ and null otherwise.

Based on the above, I made the following implementation in the $javaScript$ language, which you can test in your favorite browser.

"use strict";

// This function performs the sum of two matrices
function sumMatrix(A, B) {
    if (A.length !== B.length) {
        console.log('Problem! The inserted matrices are not the same size.');
        process.exit(0);
    } else {
        let matrixOut = [];
        for (let i = 0; i < A.length; i++) {
            matrixOut[i] = [];
            for (let j = 0; j < A[i].length; j++) {
                matrixOut[i][j] = A[i][j] + B[i][j];
            }
        }
        return matrixOut;
    }
}

// This function performs the sum of two vectors
function sumVector(A, B) {
    if (A.length !== B.length) {
        console.log('Problem! The inserted matrices are not the same size.');
        process.exit(0);
    } else {
        let matrixOut = [];
        for (let i = 0; i < A.length; i++) {
            matrixOut[i] = A[i] + B[i];
        }
        return matrixOut;
    }
}

// This function performs the subtraction of two vectors.
function subtractionVector(A, B) {
    if (A.length !== B.length) {
        console.log('Problem! The inserted matrices are not the same size.');
        process.exit(0);
    } else {
        let matrixOut = [];
        for (let i = 0; i < A.length; i++) {
            matrixOut[i] = A[i] - B[i];
        }
        return matrixOut;
    }
}

// Multiplies an nxm matrix by a mx1 vector
function multiplyMatrix(A, B) {
    if (A[0].length !== B.length) {
        console.log("Problem! The inserted matrices have different orders.");
        process.exit(0);
    } else {
        let matrixOut = [];
        for (let i = 0; i < A.length; i++) {
            matrixOut[i] = 0;
            for (let j = 0; j < A[i].length; j++) {
                matrixOut[i] += A[i][j] * B[j];
            }
        }
        return matrixOut;
    }
}

// Norm the vector
function norm(vector) {
    let sumSquare = 0;
    for (let i = 0; i < vector.length; i++) {
        sumSquare += vector[i] ** 2;
    }
    return Math.sqrt(sumSquare);
}

// Inverse matrix
function inverseMatrix(A) {
    const n = A.length;
    for (let k = 0; k < n; k++) {
        const con = A[k][k];
        A[k][k] = 1;
        for (let j = 0; j < n; j++) {
            A[k][j] = A[k][j] / con;
        }
        for (let i = 0; i < n; i++) {
            if (i !== k) {
                const con = A[i][k];
                A[i][k] = 0.0;
                for (let j = 0; j < n; j++) {
                    A[i][j] = A[i][j] - A[k][j] * con;
                }
            }
        }
    }
    return A;
}

// Matrix diagonal
function diagonalmatrix(A) {
    // Diagonal matrix
    let D = [];
    for (let i = 0; i < A.length; i++) {
        D[i] = [];
        for (let j = 0; j < A[i].length; j++) {
            if (i === j) {
                D[i][j] = A[i][j];
            } else {
                D[i][j] = 0;
            }
        }
    }
    return D;
}

// Lower triangular
function lowerMatrix(A) {
    // Lower triangular
    let L = [];
    for (let i = 0; i < A.length; i++) {
        L[i] = [];
        for (let j = 0; j < A[i].length; j++) {
            if (i > j) {
                L[i][j] = A[i][j];
            } else {
                L[i][j] = 0;
            }
        }
    }
    return L;
}

// Upper triangular
function upperMatrix(A) {
    // Upper triangular
    let U = [];
    for (let i = 0; i < A.length; i++) {
        U[i] = [];
        for (let j = 0; j < A[i].length; j++) {
            if (i < j) {
                U[i][j] = A[i][j];
            } else {
                U[i][j] = 0;
            }
        }
    }
    return U;
}

// Gauss method
function methodGauss(A, b, x0, tol) {
    let L = lowerMatrix(A);
    let D = diagonalmatrix(A);
    let U = upperMatrix(A);
    let counter = 0, x = x0, x_previous = x0;
    do {
        x_previous = x;
        x = multiplyMatrix(inverseMatrix(sumMatrix(D, L)), subtractionVector(b, multiplyMatrix(U, x)));
        counter++;
    } while (norm(subtractionVector(x, x_previous)) > tol);
    console.log(`\nNumber of iterations: ${counter}`);
    for (let i = 0; i < x.length; i++) {
        x[i] = x[i].toFixed(3);
        x[i] = parseFloat(x[i]);
    }
    return x;
}

// Enter the data
let A = [
    [4, -1, -1, 0],
    [-1, 4, 0, -1],
    [-1, 0, 4, -1],
    [0, -1, -1, 4]
];

let b = [0, 0, 1, 1];
let x0 = [1, 1, 1, 1]; // first kick

// Function call and parameter passing
let tolerance = 1.0e-8;
console.log('Solution: ', methodGauss(A, b, x0, tolerance));

In the Figure below we have the solution for the inserted data set, and also according to the indicated bibliography (Scientific Computing An Introductory Survey Heath. Example 11.5 Laplace Equation page 460/461). Hope this helps!.

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.