Почему Гаусс? (100 способов решить систему уравннений)

Posted on by tms
Что вы будете делать, если вас попросят решить простенькую систему с тремя неизвестными? У каждого сформировался свой собственный и наиболее удобный лично для него подход. Существует масса способов решить систему линейных алгебраических уравнений. Но почему предпочтение отдается именно методу Гаусса?

Обо всем по порядку

Начнем с простого. Пусть задана система линейных уравнений третьего порядка:

$$display$$\left\{ \begin{aligned} a_{11}x_1 + a_{12}x_2 + a_{13}x_3 = b_1,\\ a_{21}x_1 + a_{22}x_2 + a_{23}x_3 = b_2, \\ a_{31}x_1 + a_{32}x_2 + a_{33}x_3 = b_3. \\ \end{aligned} \right.$$display$$

Метод Гаусса заключается в последовательном «уничтожении» слагаемых, находящихся ниже главной диагонали. На первом этапе первое уравнение умножается на

$inline$ \dfrac{a_{21}}{a_{11}} $inline$

и вычитается из второго (и аналогично умножается на

$inline$ \dfrac{a_{31}}{a_{11}} $inline$

и вычитается из третьего). То есть, после этого преобразования, получаем следующее:

$$display$$\left\{ \begin{aligned} a_{11}x_1 + a_{12}x_2 + a_{13}x_3 = b_1,\\ a_{22}’x_2 + a_{23}’x_3 = b_2′, \\ a_{32}’x_2 + a_{33}’x_3 = b_3′. \\ \end{aligned} \right.$$display$$

Теперь второе уравнение умножается на

$inline$ \dfrac{a_{32}’}{a_{22}’} $inline$

и вычитается из третьего:

$$display$$\left\{ \begin{aligned} a_{11}x_1 + a_{12}x_2 + a_{13}x_3 = b_1,\\ a_{22}’x_2 + a_{23}’x_3 = b_2′, \\ a_{33}»x_3 = b_3». \\ \end{aligned} \right.$$display$$

Получили довольно простую систему, из которой легко находится

$inline$x_3$inline$

, затем

$inline$x_2$inline$

и

$inline$x_1$inline$

.

Внимательный читатель обязательно заметит: а что, если диагональные элементы равны нулю? Что же делать, если, например,

$inline$a_{11} = 0$inline$

? Неужели знаменитый метод Гаусса на этом заканчивается?

Ничего страшного! Давайте найдем

$inline$\max|a_{1j}|$inline$

и поменяем местами

$inline$j$inline$

-ую и первую строку (не ограничивая общности, предположим, что

$inline$\max |a_{1j}| = a_{13}$inline$

). Заметим, что случая, когда все

$inline$a_{1j}=0$inline$

быть не может, так как в этом случае определитель матрицы коэффициентов обращается в ноль и решить систему не предоставляется возможным. Итак, после перестановки 3-го уравнение на первую строку, выполняем действия, описанные ранее.

Поиском максимального по модулю элемента можно заниматься на каждой итерации, то есть на

$inline$k$inline$

-ой итерации искать

$inline$\max |a_{kj}|$inline$

, затем менять

$inline$j$inline$

-ую и

$inline$k$inline$

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

Есть и другой способ. Можно на

$inline$k$inline$

-ой итерации искать

$inline$\max |a_{ik}|$inline$

, затем менять уже не строчки, а столбцы. Но при этом важно запоминать индексы меняющихся столбцов в какой-нибудь массив

$inline$\alpha$inline$

(чтобы потом восстановить точный порядок переменных).

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

import java.io.FileReader; import java.io.IOException; import java.io.PrintWriter; import java.util.ArrayList; import java.util.Collections; import java.util.Locale; import java.util.Scanner;  public class GuassianEliminationSearchMainElementsAtString { 	public static void main(String[] args) throws IOException { 		 		Scanner sc = new Scanner(new FileReader("input.txt")); 		sc.useLocale(Locale.US); 		 		int n = sc.nextInt(); 		double[][] a = new double[n + 1][n + 1]; 		double[] b = new double[n + 1]; 		 		// input 		for (int i = 1; i <= n; i++) {  			for (int j = 1; j <= n; j++) { 				a[i][j] = sc.nextDouble(); 			} 			b[i] = sc.nextDouble(); 		}  		int[] alpha = new int[n + 1]; // array of indices 		for (int i = 1; i <= n; i++) { 			alpha[i] = i; 		}  		for (int m = 1; m <= n; m++) { 			double max = Math.abs(a[m][m]); 			int count = m; 			for (int i = m + 1; i <= n; i++) { 				if (Math.abs(a[m][i]) > max) { // search max elements at the string 					max = Math.abs(a[m][i]); 					count = i; 				} 			}  			int tmp = alpha[m]; // swap strings 			alpha[m] = alpha[count]; 			alpha[count] = tmp;  			for (int i = m; i <= n; i++) { 				double tmp2 = a[i][m]; 				a[i][m] = a[i][count]; 				a[i][count] = tmp2; 			}  			for (int i = m + 1; i <= n; i++) { // guassian right stroke 				b[i] = b[i] - a[i][m] * b[m] / a[m][m]; 				for (int j = m + 1; j < n; j++) { 					a[i][j] = a[i][j] - a[i][m] * a[m][j] / a[m][m]; 				} 			}  		} // for m   		double[] x = new double[n+1]; 		for (int i = n; i >= 1; i--) { // guassian back stroke 			double sum = 0; 			for (int j = i + 1; j <= n; j++) { 				sum += a[i][j] * x[alpha[j]]; 			} 			x[alpha[i] - 1] = (b[i] - sum) / a[i][i]; 			 		}  		// output 		PrintWriter pw = new PrintWriter("output.txt"); 		for (int i = 0; i < n; i++) { 			pw.printf(Locale.US, "x%d = %.5f \n", i + 1, x[i]); 		}  		pw.flush(); 		pw.close();  	} }

Почему Гаусс?

Существует и другой способ решения СЛАУ. Один из таких — метод Крамера. Он заключается в предварительном подсчете некоторого количества определителей, с помощью которых моментально вычисляются значения переменных. При исходной системе этот метод будет выглядеть следующим образом:

$$display$$ \Delta = \begin{vmatrix} a_{11} & a_{12} & a_{13}\\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33}\\ \end{vmatrix}, \Delta_1 = \begin{vmatrix} b_1 & a_{12} & a_{13}\\ b_2 & a_{22} & a_{23} \\ b_3 & a_{32} & a_{33}\\ \end{vmatrix}, $$display$$

$$display$$ \Delta_2 = \begin{vmatrix} a_{11} & b_1 & a_{13}\\ a_{21} & b_2 & a_{23} \\ a_{31} & b_3 & a_{33}\\ \end{vmatrix}, \Delta_3 = \begin{vmatrix} a_{11} & a_{12} & b_1\\ a_{21} & a_{22} & b_2 \\ a_{31} & a_{32} & b_3\\ \end{vmatrix}, $$display$$

$$display$$ x_i = \dfrac{\Delta_i}{\Delta}.$$display$$

Но вспомним — что такое определитель?

По определению, определитель матрицы

$inline$A = (a_{ij})$inline$

есть сумма

$$display$$ \sum\limits_{1\leq i_1 < \dots < i_n \leq n} (-1)^{N(i_1, \dots, i_n)} a_{1i_1}\dots a_{ni_n}, $$display$$

где

$inline$N(i_1,\dots, i_n)$inline$

— знак подстановки

$inline$i_1, \dots, i_n.$inline$

Определитель содержит

$inline$n!$inline$

слагаемых. Для решения системы необходимо посчитать

$inline$n+1$inline$

определителей. При достаточно больших

$inline$n$inline$

это очень затратно. Например, при

$inline$n = 12$inline$

число операций становится

$inline$12!(12+1) = 6227020800,$inline$

в то время как метод Гаусса с ассимптотикой

$inline$n^3$inline$

потребует всего лишь

$inline$12^3 = 1728$inline$

операций.

Итерационные методы

В качестве алгоритмов решения СЛАУ подходят и так называемые итерационные методы. Они заключаются в последовательном вычислении приближений до тех пор, пока какое-то из них будет максимально близко к точному ответу.

Сначала выбирается какой-то произвольный вектор

$inline$x^0$inline$

— это нулевое приближение. По нему строится вектор

$inline$x^1$inline$

— это первое приближение. И так далее. Вычисления заканчиваются, когда

$inline$||x^k — x^{k+1}|| < \varepsilon$inline$

, где

$inline$\varepsilon$inline$

— какое-то заданное наперед значение. Последнее неравенство означает, что наше «улучшение» решения с каждой итерацией получается почти незначительным.

Рассмотрим популярный метод Якоби, который является одним из простейших итерационных методов решения СЛАУ.

Для начала запишем систему в следующем виде:

$$display$$ \sum\limits_{j\leq n} a_{ij}x_j = b_i, \ i = \overline{1,n}. $$display$$

Отделим

$inline$i$inline$

-ое слагаемое и выразим его через все остальное:

$$display$$ x_i = \dfrac{b_i — \sum\limits_{j\neq i} a_{ij}x_j}{a_{ii}}, \ i = \overline{1,n}.$$display$$

Теперь просто навесим «счетчики» на переменные и получим итерационный метод Якоби:

$$display$$ x_i^k = \dfrac{b_i — \sum\limits_{j\neq i} a_{ij}x_j^k}{a_{ii}}, \ i = \overline{1,n},\ k = 0,1,\dots.$$display$$

Заметим, что обязательным условием применения данного метода является отсутствие нулей по главной диагонали.

Реализация метода Якоби на Java:
В качестве

$inline$\varepsilon$inline$

берется заранее вычисленное так называемое машинное эпсилон.

import java.io.FileNotFoundException; import java.io.FileReader; import java.io.PrintWriter; import java.util.Locale; import java.util.Scanner;  public class JacobiMethod {  	public static void main(String[] args) throws FileNotFoundException {  		Scanner sc = new Scanner(new FileReader("input.txt")); 		sc.useLocale(Locale.US);  		int n = sc.nextInt();  		double[][] a = new double[n + 1][n + 1]; 		double[] b = new double[n + 1]; 		double[] x0 = new double[n + 1];  		for (int i = 1; i <= n; i++) { 			for (int j = 1; j <= n; j++) { 				a[i][j] = sc.nextDouble(); 			} 			b[i] = sc.nextDouble(); 			x0[i] = b[i] / a[i][i]; 		} 				 		double EPS = EPSCalc(); 		double[] x = new double[n+1]; 		double norm = Double.MAX_VALUE; 		int counter = 0; 		do{ 			for(int i = 1; i <= n; i++) { 				double sum = 0; 				for(int j = 1; j <= n; j++) { 					if(j == i) continue; 					sum += a[i][j] * x0[j]; 				} 				x[i] = (b[i] - sum) / a[i][i]; 			} 			norm = normCalc(x0, x, n); 			for(int i = 1; i <= n; i++) { 				x0[i] = x[i]; 			} 			counter++; 		} while(norm > EPS); 		 		PrintWriter pw = new PrintWriter("output.txt"); 		pw.println(counter + " iterations"); 		for (int i = 1; i <= n; i++) { 			pw.printf(Locale.US, "x%d = %f\n", i, x0[i]); 		} 		pw.flush(); 		pw.close();  	} 	 	static double normCalc(double []x1, double[] x2, int n) { 		double sum = 0; 		for(int i = 1; i <= n; i++) { 			sum += Math.abs(x1[i] - x2[i]); 		} 		return sum; 	} 	 	static double EPSCalc () { 		double eps = 1;         while (1 + eps > 1) {             eps /= 2;         }         return eps; 	}  }

Модификацией метода Якоби является метод релаксации. Его главное отличие заключается в том, что с помощью заранее подобранного параметра количество итераций значительно снижается. Опишем в кратце главную идею метода.

Из исходной системы снова выразим

$inline$x$inline$

, но расставим немного иначе счетчики и добавим параметр

$inline$\omega$inline$

:

$$display$$ x_i^k = \dfrac{\omega\left(b_i — \sum\limits_{j = 1}^{i-1}a_{ij}x_j^{k+1} — \sum\limits_{j = i+1}^n a_{ij}x_j^k\right)}{a_{ii}} + (1-\omega)x_i^k, \ i = \overline{1,n},\ k = 0,1,\dots.$$display$$

При

$inline$\omega=1$inline$

это все превращается в метод Якоби.

Итак, будем искать какое-нибудь «хорошее»

$inline$\omega$inline$

. Зададим какое-нибудь число

$inline$s$inline$

и будем беребирать значения

$inline$\omega \in (0,2)$inline$

, для каждого из которых будем считать нормы

$inline$||x^{k+1}-x^k||, \ k = \overline{1,s}$inline$

. Для наименьшей из этих норм запомним данное значение

$inline$\omega$inline$

, и с его помощью будем решать нашу систему.

Иллюстрация метода на языке Java:
здесь

$inline$s=5$inline$ 

import java.io.FileNotFoundException; import java.io.FileReader; import java.io.PrintWriter; import java.util.Locale; import java.util.Scanner;  public class SuccessiveOverRelaxation { 	 	public static void main(String[] args) throws FileNotFoundException { 		 		Scanner sc = new Scanner(new FileReader("input.txt")); 		sc.useLocale(Locale.US);  		int n = sc.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] = sc.nextDouble(); 			} 			b[i] = sc.nextDouble(); 		} 		 		double EPS = EPSCalc(); 		 		double w = bestRelaxationParameterCalc(a, b, n); 		         double[] x = new double[n + 1];                  int counter = 0;         double maxChange = Double.MAX_VALUE;         do {         	maxChange = 0;             for (int i = 1; i <= n; i++) {                  double firstSum = 0;                 for (int j = 1; j <= i - 1; j++) {                     firstSum += a[i][j] * x[j];                 }                  double secondSum = 0;                 for (int j = i + 1; j <= n; j++) {                     secondSum += a[i][j] * x[j];                 }                  double lastTerm = (1 - w) * x[i];                 double z = (b[i] - firstSum - secondSum);                 double temp = (w * z) / a[i][i] + lastTerm ;                 maxChange = Math.max(maxChange, Math.abs(x[i] - temp));                 x[i] = temp;             }             counter++;         } while(maxChange > EPS);                  PrintWriter pw = new PrintWriter("output.txt");         pw.println(w + " is the best relaxation parameter"); 		pw.println(counter + " iterations"); 		for (int i = 1; i <= n; i++) { 			pw.printf(Locale.US, "x%d = %f\n", i, x[i]); 		} 		pw.flush(); 		pw.close(); 		 	} 	 	static double bestRelaxationParameterCalc(double[][]a, double[]b, int n) { 		 	double bestW = 1, bestMaxChange = Double.MAX_VALUE;         for (double w = 0.05; w <= 2; w += 0.05) {              double[] x = new double[n + 1];                          double maxChange = 0;             for (int s = 0; s < 5; s++) {                 maxChange = 0;                 for (int i = 1; i <= n; i++) {                      double firstSum = 0;                     for (int j = 1; j <= i - 1; j++) {                         firstSum += a[i][j] * x[j];                     }                      double secondSum = 0;                     for (int j = i + 1; j <= n; j++) {                         secondSum += a[i][j] * x[j];                     }                      double lastTerm = (1 - w) * x[i];                     double z = (b[i] - firstSum - secondSum);                     double temp = (w * z) / a[i][i] + lastTerm;                     maxChange = Math.max(maxChange, Math.abs(x[i] - temp));                     x[i] = temp;                 }             }              if (maxChange < bestMaxChange) {                 bestMaxChange = maxChange;                 bestW = w;             }          }                  return bestW;          	} 	 	static double EPSCalc () { 		double eps = 1;         while (1 + eps > 1) {             eps /= 2;         }         return eps; 	} 	 }

Вместо заключения

Существует еще масса алгоритмов для решения СЛАУ. Например, метод квадратного корня, в котором искомая система заменяется на две «простых» системы, решения которых вычисляется по элементарным формулам; метод прогонки, который используется для так специфических трехдиагональных матриц. Каждый сам выбирает, какой метод ему использовать для своей проблемы.

FavoriteLoadingДобавить в избранное
Posted in Без рубрики

Добавить комментарий