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 $$$LU$decomposition of the matrix.

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

$$

$$Ax = b,$$

Where $$$A = (a_ {ij})$- terrible and confusing non-degenerate matrix. Imagine it as a product of two other matrices: $$$L = (l_ {ij})$- lower triangular matrix, and $$$U = (u_ {ij})$- upper triangular matrix. Then obviously the system takes the form:

$$

$$LUx = b.$$

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

$$

$$Ly = b,$$

$$

$$Ux = 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 $$$LU$decomposition for the matrix.

Description of the algorithm that applies $$$LU$- decomposition is described in some detail here . 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_ {ij} = a_ {ji}$or $$$A ^ t = A$. Such a matrix can be represented as

$$

$$A = S ^ tDS,$$

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

$$

$$S ^ tDy = b,$$

$$

$$Sx = 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 found here .

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; public class SquareRootMethod { public static void main(String[] args) throws FileNotFoundException { Scanner scanner = new Scanner(new FileReader("input.txt")); scanner.useLocale(Locale.US); int n = scanner.nextInt(); double[][] a = new double[n + 1][n + 1]; double[] b = new double[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 = new double[n + 1]; double[] d = new double[n + 1]; double[][] s = new double[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 matrix for (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 = b double[] y = new double[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

$$

$$Ax = b$$

matrix $$$A$is tridiagonal:

$$

$$A = \ begin {pmatrix} C_1 & -B_1 & 0 & \ dots & 0 & 0 \\ -A_2 & C_2 & -B_2 & \ dots & 0 & 0 \\ vdots & \ vdots & \ vdots & \ ddots & \ vdots & \ vdots \\ 0 & 0 & \ dots & -A_ {n-1} & C_ {n-1} & -B_ {n-1} \\ 0 & 0 & \ dots & 0 & -A_n & C_n \\ \ end {pmatrix}.$$

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

$$

$$x_ {i-1} = \ alpha_ {i-1} x_i + \ beta_ {i-1},$$

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

Derivation of formulas

$$

$$\ alpha_1 = \ dfrac {B_1} {C_1}, \ \ alpha_i = \ dfrac {B_i} {C_i - A_i \ alpha_ {i-1}}, \ i = \ overline {2, n},$$

$$

$$\ beta_1 = \ dfrac {b_1} {C_1}, \ \ beta_i = \ dfrac {b_i + A_i \ beta_ {i-1}} {C_i - A_i \ alpha_ {i-1}}, \ i = \ overline { 2, n}$$

present here . Well, in the end

$$

$$x_n = \ dfrac {b_n + A_n \ beta_ {n-1}} {C_n - A_n \ alpha_ {n-1}}, \ x_i = \ alpha_i x_ {i + 1} + \ beta_i, \ i = n- 1, n-2, \ dots, 1$$

Note that in the search formulas $$$\ alpha_i, \ beta_i$there is a division by number $$$C_i - A_i \ alpha_ {i-1},$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:

$$

$$C_1, C_n \ neq 0, A_i, B_i \ neq 0 \ (i = \ overline {2, n}),$$

$$

$$| C_i | \ geq | A_i | + | B_i |.$$

Consider a Java program code that solves the system.

$$

$$\ left \ {\ begin {aligned} & x_1 = 0, \\ & x_ {i-1} + \ xi x_ {i + 1} = \ dfrac {i} {n}, \ (i = \ overline {2, n-1}), \\ & x_n = 1, \\ \ end {aligned} \ right.$$

at $$$\ xi = 2, n = 21$.

 package SystemOfLinearAlgebraicEquations; import java.io.FileNotFoundException; import java.io.FileReader; import java.io.PrintWriter; import java.util.Locale; import java.util.Scanner; public class TridiagonalMatrixAlgorithm { public static void main(String[] args) throws FileNotFoundException { final int n = 21; double[] A = new double[n + 1]; double[] B = new double[n + 1]; double[] C = new double[n + 1]; double[] F = new double[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 = new double[n + 1]; double[] beta = new double[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 = new double[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(); } }