Обо всем по порядку
Начнем с простого. Пусть задана система линейных уравнений третьего порядка:
Метод Гаусса заключается в последовательном «уничтожении» слагаемых, находящихся ниже главной диагонали. На первом этапе первое уравнение умножается на и вычитается из второго (и аналогично умножается на и вычитается из третьего). То есть, после этого преобразования, получаем следующее:
Теперь второе уравнение умножается на и вычитается из третьего:
Получили довольно простую систему, из которой легко находится , затем и .
Внимательный читатель обязательно заметит: а что, если диагональные элементы равны нулю? Что же делать, если, например, ? Неужели знаменитый метод Гаусса на этом заканчивается?
Ничего страшного! Давайте найдем и поменяем местами -ую и первую строку (не ограничивая общности, предположим, что ). Заметим, что случая, когда все быть не может, так как в этом случае определитель матрицы коэффициентов обращается в ноль и решить систему не предоставляется возможным. Итак, после перестановки 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;
}
}
Вместо заключения
Существует еще масса алгоритмов для решения СЛАУ. Например, метод квадратного корня, в котором искомая система заменяется на две «простых» системы, решения которых вычисляется по элементарным формулам; метод прогонки, который используется для так специфических трехдиагональных матриц. Каждый сам выбирает, какой метод ему использовать для своей проблемы.
Комментарии (3)
ProstoUser
25.07.2018 19:46+11Давным давно в
далекой галактикеодном научном институте пришлось мне программировать решение системы линейных уравнений как небольшую и (как тогда казалось) тривиальную часть другой задачи.
Но почему-то метод Крамера, который я быстренько написал, давал совершенно идиотские результаты — подстановка решения в систему не давало нулей справа. Оказалось — всего лишь матрица была слабо обусловленной.
Мой научный руководитель предложил посмотреть библиотеку математических методов и поискать там. Так вот, в библиотек оказалось несколько десятков уже запрограммированных методов. Каждый — наилучшим образом подходящий для решения систем с какими-то определенными свойствами.
Это я к тому, что задача решения систем линейных уравнений не так проста, как кажется на первый взгляд. И проблем там на пару порядков больше, чем упомянуто в этой статье. И выбор правильного метода — это не «кому какой нравится», а сначала анализ своей системы уравнений, и уже на основе этого анализа подбор метода, который даст хороший результат именно в этом случае. Выбор же неправильного метода приведет к тому, что решение будет просто неправильным.
В общем, это не про программирование, а, скорее, про математику. В общем, тема не так, чтобы раскрыта :-)
Ildarovich
26.07.2018 15:37Итерационные алгоритмы практичнее, они эффективны по памяти и работают на больших и очень больших (разреженных) матрицах. Гаусс не выдерживает большого количества переменных. Крамер тем более. Управлять погрешностями тоже непросто. Сложнее распараллеливается. Так что правильнее говорить: «почему НЕ Гаусс?».
aamonster
Наверное, минуснул бы, если бы мог.
Ответ на вопрос в заголовке не требует статьи. "Потому что он самый простой и понятный".
А вот упомянув другие методы — вы вообще не сказали, зачем они нужны, раз есть такой чудесный метод Гаусса. Не упомянули его недостатки, из-за которых эти методы приходится искать (проблема накопления погрешности при плохо обусловленной матрице и иногда даже при нормальной; сложность выбора ведущего элемента). Не сказали, в каких условиях какой метод нужен. Много букв ни о чём.