Почему Гаусс? (100 способов решить систему уравннений)
Обо всем по порядку
Начнем с простого. Пусть задана система линейных уравнений третьего порядка:
Метод Гаусса заключается в последовательном «уничтожении» слагаемых, находящихся ниже главной диагонали. На первом этапе первое уравнение умножается на
и вычитается из второго (и аналогично умножается на
и вычитается из третьего). То есть, после этого преобразования, получаем следующее:
Теперь второе уравнение умножается на
и вычитается из третьего:
Получили довольно простую систему, из которой легко находится
, затем
и
.
Внимательный читатель обязательно заметит: а что, если диагональные элементы равны нулю? Что же делать, если, например,
? Неужели знаменитый метод Гаусса на этом заканчивается?
Ничего страшного! Давайте найдем
и поменяем местами
-ую и первую строку (не ограничивая общности, предположим, что
). Заметим, что случая, когда все
быть не может, так как в этом случае определитель матрицы коэффициентов обращается в ноль и решить систему не предоставляется возможным. Итак, после перестановки 3-го уравнение на первую строку, выполняем действия, описанные ранее.
Поиском максимального по модулю элемента можно заниматься на каждой итерации, то есть на
-ой итерации искать
, затем менять
-ую и
-ую строчки. Алгоритм, в которм осуществляется поиск максимального элемента в столбце, называется методом Гаусса с выбором главного элемента в столбце.
Есть и другой способ. Можно на
-ой итерации искать
, затем менять уже не строчки, а столбцы. Но при этом важно запоминать индексы меняющихся столбцов в какой-нибудь массив
(чтобы потом восстановить точный порядок переменных).
Пример простого кода, реализующего данный алгоритм:
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(); } }
Почему Гаусс?
Существует и другой способ решения СЛАУ. Один из таких — метод Крамера. Он заключается в предварительном подсчете некоторого количества определителей, с помощью которых моментально вычисляются значения переменных. При исходной системе этот метод будет выглядеть следующим образом:
Но вспомним — что такое определитель?
По определению, определитель матрицы
есть сумма
где
— знак подстановки
Определитель содержит
слагаемых. Для решения системы необходимо посчитать
определителей. При достаточно больших
это очень затратно. Например, при
число операций становится
в то время как метод Гаусса с ассимптотикой
потребует всего лишь
операций.
Итерационные методы
В качестве алгоритмов решения СЛАУ подходят и так называемые итерационные методы. Они заключаются в последовательном вычислении приближений до тех пор, пока какое-то из них будет максимально близко к точному ответу.
Сначала выбирается какой-то произвольный вектор
— это нулевое приближение. По нему строится вектор
— это первое приближение. И так далее. Вычисления заканчиваются, когда
, где
— какое-то заданное наперед значение. Последнее неравенство означает, что наше «улучшение» решения с каждой итерацией получается почти незначительным.
Рассмотрим популярный метод Якоби, который является одним из простейших итерационных методов решения СЛАУ.
Для начала запишем систему в следующем виде:
Отделим
-ое слагаемое и выразим его через все остальное:
Теперь просто навесим «счетчики» на переменные и получим итерационный метод Якоби:
Заметим, что обязательным условием применения данного метода является отсутствие нулей по главной диагонали.
Реализация метода Якоби на Java:
В качестве
берется заранее вычисленное так называемое машинное эпсилон.
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; } }
Модификацией метода Якоби является метод релаксации. Его главное отличие заключается в том, что с помощью заранее подобранного параметра количество итераций значительно снижается. Опишем в кратце главную идею метода.
Из исходной системы снова выразим
, но расставим немного иначе счетчики и добавим параметр
:
При
это все превращается в метод Якоби.
Итак, будем искать какое-нибудь «хорошее»
. Зададим какое-нибудь число
и будем беребирать значения
, для каждого из которых будем считать нормы
. Для наименьшей из этих норм запомним данное значение
, и с его помощью будем решать нашу систему.
Иллюстрация метода на языке Java:
здесь
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; } }
Вместо заключения
Существует еще масса алгоритмов для решения СЛАУ. Например, метод квадратного корня, в котором искомая система заменяется на две «простых» системы, решения которых вычисляется по элементарным формулам; метод прогонки, который используется для так специфических трехдиагональных матриц. Каждый сам выбирает, какой метод ему использовать для своей проблемы.

