More about methods for solving systems of linear algebraic equations

It is of course very difficult to tell in detail about all the methods, but this topic seems to me interesting and extremely important, because everyone often comes across the task of finding a solution. In the first article Why Gauss? The Gauss method (including modifications) and some iterative methods were described. However, given the criticism of Sinn3r , I decided to describe other methods.

Start over

So, let us again need to solve a system of linear algebraic equations of order $$n . One of the most striking numerical methods is the method that uses$$L U decomposition of the matrix.

This intuitive method is as follows. Suppose we have a system of linear algebraic equations (written in vector form):

$$A x = b ,

Where $$A = ( a i j ) is a terrible and entanglednon-degeneratematrix. Imagine it as a product of two other matrices:$$L = ( l i j ) is the lower triangular matrix, and$$U = ( u i j ) is the upper triangular matrix. Then obviously the system takes the form:

$$L U x = b .

If to designate $$U x = y , then the solution of the source system is found from the solution of two other systems:

$$L y = b ,

$$U x = y .

The advantage of this method is that the solutions of the two auxiliary systems are found very simply (due to the fact that the matrices $$L and$$U - triangular). But is it easy to find these matrices? So, the problem of solving a system is reduced to the problem of constructing$$L U decompositions for the matrix.

Description of the algorithm that applies$$L U decomposition is described in sufficient detailhere. And while we look at another method.

Square root method

We impose additional conditions on the matrix $$A . We require that it be symmetric, that is,$$a i j = a j i or$$A t = A . Such a matrix can be represented as

$$A = S t D S ,

Where $$S is the upper triangular matrix,$$D is the diagonalrealmatrix, and its diagonal elements are equal to unity in absolute value. Similarly, the original system is reduced to two others:

$$S t D y = b ,

$$S x = y ,

which are also solved by elementary properties of matrices $$S and$$D .

The derivation of algorithmic formulas is rather cumbersome and remains outside the scope of the publication. If desired, they can be foundhere.

Sample Java code that implements the sweep method:

import java.io.FileNotFoundException;
import java.io.FileReader;
import java.io.PrintWriter;
import java.util.Locale;
import java.util.Scanner;
publicclassSquareRootMethod{
    publicstaticvoidmain(String[] args)throws FileNotFoundException {
        Scanner scanner = new Scanner(new FileReader("input.txt"));
        scanner.useLocale(Locale.US);
        int n = scanner.nextInt();
        double[][] a = newdouble[n + 1][n + 1];
        double[] b = newdouble[n + 1];
        for (int i = 1; i <= n; i++) {
            for (int j = 1; j <= n; j++) {
                a[i][j] = scanner.nextDouble();
            }
            b[i] = scanner.nextDouble();
        }
        double[] x = newdouble[n + 1];
        double[] d = newdouble[n + 1];
        double[][] s = newdouble[n + 1][n + 1];
        // for i = 1
        d[1] = Math.signum(a[1][1]);
        s[1][1] = Math.sqrt(Math.abs(a[1][1]));
        for (int j = 2; j <= n; j++) {
            s[1][j] = a[1][j] / (s[1][1] * d[1]);
        }
        // for i > 1//searching S and D matrixfor (int i = 2; i <= n; i++) {
            double sum = 0;
            for (int k = 1; k <= i - 1; k++) {
                sum += s[k][i] * s[k][i] * d[k];
            }
            d[i] = Math.signum(a[i][i] - sum);
            s[i][i] = Math.sqrt(Math.abs(a[i][i] - sum));
            double l = 1 / (s[i][i] * d[i]);
            for (int j = i + 1; j <= n; j++) {
                double SDSsum = 0;
                for (int k = 1; k <= i - 1; k++) {
                    SDSsum += s[k][i] * d[k] * s[k][j];
                }
                s[i][j] = l * (a[i][j] - SDSsum);
            }
        }
        //solve of the system (s^t * d)y = bdouble[] y = newdouble[n + 1];
        y[1] = b[1] / (s[1][1] * d[1]);
        for (int i = 2; i <= n; i++) {
            double sum = 0;
            for (int j = 1; j <= i - 1; j++) {
                sum += s[j][i] * d[j] * y[j];
            }
            y[i] = (b[i] - sum) / (s[i][i] * d[i]);
        }
        //solve of the system sx = y
        x[n] = y[n] / s[n][n];
        for (int i = n - 1; i >= 1; i--) {
            double sum = 0;
            for (int k = i + 1; k <= n; k++) {
                sum += s[i][k] * x[k];
            }
            x[i] = (y[i] - sum) / s[i][i];
        }
        //output
        PrintWriter pw = new PrintWriter("output.txt");
        for (int i = 1; i <= n; i++) {
            pw.printf(Locale.US, "x%d = %f\n", i, x[i]);
        }
        pw.flush();
        pw.close();
    }
}

Sweep method

Using this method, only specific systems can be solved, having no more than three unknowns in each row. That is, with the system

$$A x = b

matrix $$A is tridiagonal:

$$A = ( C 1 - B 1 0 ... 0 0 - A 2 C 2 - B 2 ... 0 0 ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ 0 0 ... - A n - 1 C n - 1 - B n - 1 0 0 ... 0 - A n C n ) .

Immediately, we note that there is a connection between neighboring solutions:

$$

$$- some unknown numbers. If we find them and any one variable, we can find all the others.

Derivation of formulas

$$

$$

present here . Well, in the end

$$

Note that in the search formulas $$ there is a division by number $$which may turn out to be zero to keep track of. But in fact the following assertion holds , the proof of which is here : the sweep algorithm is correct and stable if the conditions are fulfilled:

$$

$$

Consider a Java program code that solves the system.

$$

at $$.

package SystemOfLinearAlgebraicEquations;
import java.io.FileNotFoundException;
import java.io.FileReader;
import java.io.PrintWriter;
import java.util.Locale;
import java.util.Scanner;
publicclassTridiagonalMatrixAlgorithm{
	publicstaticvoidmain(String[] args)throws FileNotFoundException {
		finalint n = 21;
		double[] A = newdouble[n + 1];
		double[] B = newdouble[n + 1];
		double[] C = newdouble[n + 1];
		double[] F = newdouble[n + 1];
		double xi = 2;
		C[1] = 1;
		F[1] = 0;
		for(int i = 2; i <= n-1; i++) {
			A[i] = 1;
			C[i] = xi;
			B[i] = 1;
			F[i] =  - (double) i / n;
		}
		A[n] = 0;
		C[n] = 1;
		F[n] = 1;
		double[] alpha = newdouble[n + 1];
		double[] beta = newdouble[n + 1];
		// right stroke
		alpha[1] = B[1] / C[1];
		beta[1] = F[1] / C[1];
		for (int i = 2; i <= n - 1; i++) {
			double denominator = C[i] - A[i] * alpha[i-1];
			alpha[i] = B[i] / denominator;
			beta[i] = (F[i] + A[i] * beta[i - 1]) / denominator;
		}
		double[] x = newdouble[n + 1];
		// back stroke
		x[n] = (F[n] + A[n] * beta[n - 1]) / (C[n] - A[n] * alpha[n - 1]);
		for (int i = n - 1; i >= 1; i--) {
			x[i] = alpha[i] * x[i + 1] + beta[i];
		}
		// output
		PrintWriter pw = new PrintWriter("output.txt");
		for (int i = 1; i <= n; i = i + 5) {
			pw.printf(Locale.US, "x%d = %f\n", i, x[i]);
		}
		pw.flush();
		pw.close();
	}
}