Еще про методы решения систем линейных алгебраических уравнений

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

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

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

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

$$

где $$ — страшная и запутанная невырожденная матрица. Представим ее в виде произведения двух других матриц: $$ — нижняя треугольная матрица, и $$ — верхняя треугольная матрица. Тогда очевидно система принимает вид:

$$

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

$$

$$

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

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

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

Наложим дополнитльные условия на матрицу $$. Потребуем, чтобы она была симметрическая, то есть $$ или $$. Такую матрицу можно представить в виде

$$

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

$$

$$

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

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

Пример кода на 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();              }  } 

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

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

$$

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

$$

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

$$

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

Вывод формул

$$

$$

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

$$

Отметим, что в формулах поиска $$ присутствует деление на число $$, которое может оказаться нулем, что нужно отслеживать. Но на самом деле имеет место следующее утверждение, доказательство которого есть здесь: алгоритм прогонки является корректным и устойчивым, если выполняются условия:

$$

$$

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

$$

при $$.

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