Рассказать подробно про все методы конечно же очень трудно, но мне эта тема кажется интересной и чрезвычайно важной, поскольку с задачей нахождения решения все сталкиваются достаточно часто. В первой статье Почему Гаусс? был описан метод Гаусса (в том числе с модификацими) и некоторые итерационные методы. Однако, учитывая критику Sinn3r, я решил описать и другие методы.

Начнем сначала


Итак, пусть нам снова нужно решить систему линейных алгебраических уравнений порядка $n$. Одним из наиболее ярких численных методов является метод, который использует $LU$-разложение матрицы.

Это интуитивно понятный метод заключается в следующем. Пусть у нас есть система линейных алгебраических уравнений (записанная в вектором виде):

$ Ax = b,$

где $A = (a_{ij})$ — страшная и запутанная невырожденная матрица. Представим ее в виде произведения двух других матриц: $L = (l_{ij})$ — нижняя треугольная матрица, и $U=(u_{ij})$ — верхняя треугольная матрица. Тогда очевидно система принимает вид:

$ LUx = b.$

Если обозначить $Ux = y$, то решение исходной системы находится из решения двух других систем:

$ Ly = b,$

$Ux = y.$


Преимущество этого метода в том, что решения двух вспомогательных систем находится очень просто (в силу того, что матрицы $L$ и $U$ — треугольные). Но легко ли найти эти матрицы? Итак, задача решения системы сводится к задаче построения $LU$-разложения для матрицы.

Описание алгоритма, применяющего $LU$-разложение достаточно подробно описано здесь. А мы пока посмотрим на другой метод.

Метод квадратого корня


Наложим дополнитльные условия на матрицу $A$. Потребуем, чтобы она была симметрическая, то есть $a_{ij} = a_{ji}$ или $A^t = A$. Такую матрицу можно представить в виде

$A = S^tDS,$

где $S$ — верхняя треугольная матрица, $D$ — диагональная вещественная матрица, причем ее диагональные элементы по модулю равны единице. Аналогичным образом исходная система сводится к двум другим:

$S^tDy = b,$

$Sx = y,$

которые тоже решаются элементарно в силу свойств матриц $S$ и $D$.

Вывод алгоритмических формул довольно громоздкий и остается за рамками публикации. При желании, их можно найти здесь.

Пример кода на Java, реализующий метод прогонки:

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();
        
    }

}

Метод прогонки


С помощью этого метода можно решать только специфические системы, имеющие не более трех неизвестных в каждой строке. То есть при системе

$Ax = b$

матрица $A$ является трехдиагональной:

$A = \begin{pmatrix} С_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}. $



Сразу заметим, что имеется связь соседних решений:

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


$\alpha_k, \beta_k$ — какие-то неизвестные числа. Если мы найдем их и какую-то одну переменную, то сможем найти и все остальные.

Вывод формул

$\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}$

присутствует здесь. Ну и в итоге

$ 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$


Отметим, что в формулах поиска $\alpha_i, \beta_i$ присутствует деление на число $C_i - A_i\alpha_{i-1},$, которое может оказаться нулем, что нужно отслеживать. Но на самом деле имеет место следующее утверждение, доказательство которого есть здесь: алгоритм прогонки является корректным и устойчивым, если выполняются условия:

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

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


Рассмотрим код программы на Java, который решает систему

$\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.$

при $\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();

	}

}

Комментарии (4)


  1. third112
    30.07.2018 13:49

    Рассказать подробно про все методы конечно же очень трудно, но мне эта тема кажется интересной
    Ok. Действительно трудно. Действительно интересно.
    Нпр., метод Крамера. В вики сказано:
    метод, с точки зрения затрат времени на вычисления, считался непрактичным. Однако в 2010 году было показано, что метод Крамера может быть реализован со сложностью O ( n 3 ) {\displaystyle O(n^{3})} O(n^{3}), сравнимой со сложностью метода Гаусса

    Но ИМХО интересно и теоретически: если коэффициенты СЛАУ рациональные числа, то корни (если есть) будут рациональными, т.к. метод использует только 4 арифметических действия. Для программирования это м.б. очень полезно: у меня была задача, где нужно было сравнить корни 2х систем на равенство. Возникала ошибка округления, а с рациональными заработало. (Решал Гауссом).


  1. dannk
    30.07.2018 22:21
    +1

    Думаю стоило бы упомянуть, что LU — это тот же Гаусс, только сбоку, а метод квадратного корня обычно называют разложением Холецкого


    1. malkovsky
      31.07.2018 12:42

      Тут даже может возникнуть путаница. Есть метод Холецкого (Cholesky decomposition) — разложение вида A=LL^T, которое существует только для симметричных матриц, а есть матричный корень — разложение вида A=BB, которое обладает несколько другими свойствами. В статье вроде бы разложение Холесского описано


  1. malkovsky
    31.07.2018 12:52

    Кажется, у вас код про «Метод квадратного корня» (видимо разложение Холецкого) почему-то подписан как код к «методу прогонки».

    Вообще, разложение Холецкого, хоть это и не очевидно, — это аналог Гаусса для симметричных матриц. Вот тут на слайде 12-15 описана абстрактная рекурсивная реализация (ну еще вот тут на 20 слайде, суть та же). Если Холецкий применим, то он предпочтительней Гаусса.