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!.
