Как стать автором
Поиск
Написать публикацию
Обновить

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

Время на прочтение6 мин
Количество просмотров13K
Алгоритмы*Математика*
Из песочницы
Что вы будете делать, если вас попросят решить простенькую систему с тремя неизвестными? У каждого сформировался свой собственный и наиболее удобный лично для него подход. Существует масса способов решить систему линейных алгебраических уравнений. Но почему предпочтение отдается именно методу Гаусса?

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


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

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


Метод Гаусса заключается в последовательном «уничтожении» слагаемых, находящихся ниже главной диагонали. На первом этапе первое уравнение умножается на $ \dfrac{a_{21}}{a_{11}} $ и вычитается из второго (и аналогично умножается на $ \dfrac{a_{31}}{a_{11}} $ и вычитается из третьего). То есть, после этого преобразования, получаем следующее:

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


Теперь второе уравнение умножается на $ \dfrac{a_{32}'}{a_{22}'} $ и вычитается из третьего:

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


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

Внимательный читатель обязательно заметит: а что, если диагональные элементы равны нулю? Что же делать, если, например, $a_{11} = 0$? Неужели знаменитый метод Гаусса на этом заканчивается?

Ничего страшного! Давайте найдем $\max|a_{1j}|$ и поменяем местами $j$-ую и первую строку (не ограничивая общности, предположим, что $\max |a_{1j}| = a_{13}$). Заметим, что случая, когда все $a_{1j}=0$ быть не может, так как в этом случае определитель матрицы коэффициентов обращается в ноль и решить систему не предоставляется возможным. Итак, после перестановки 3-го уравнение на первую строку, выполняем действия, описанные ранее.

Поиском максимального по модулю элемента можно заниматься на каждой итерации, то есть на $k$-ой итерации искать $\max |a_{kj}|$, затем менять $j$-ую и $k$-ую строчки. Алгоритм, в которм осуществляется поиск максимального элемента в столбце, называется методом Гаусса с выбором главного элемента в столбце.

Есть и другой способ. Можно на $k$-ой итерации искать $\max |a_{ik}|$, затем менять уже не строчки, а столбцы. Но при этом важно запоминать индексы меняющихся столбцов в какой-нибудь массив $\alpha$ (чтобы потом восстановить точный порядок переменных).

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

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

	}
}

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


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

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


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


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



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

По определению, определитель матрицы $A = (a_{ij})$ есть сумма

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

где $N(i_1,\dots, i_n)$ — знак подстановки $i_1, \dots, i_n.$

Определитель содержит $n!$ слагаемых. Для решения системы необходимо посчитать $n+1$ определителей. При достаточно больших $n$ это очень затратно. Например, при $n = 12$ число операций становится $12!(12+1) = 6227020800,$ в то время как метод Гаусса с ассимптотикой $n^3$ потребует всего лишь $12^3 = 1728$ операций.

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


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

Сначала выбирается какой-то произвольный вектор $x^0$ — это нулевое приближение. По нему строится вектор $x^1$ — это первое приближение. И так далее. Вычисления заканчиваются, когда $||x^k - x^{k+1}|| < \varepsilon$, где $\varepsilon$ — какое-то заданное наперед значение. Последнее неравенство означает, что наше «улучшение» решения с каждой итерацией получается почти незначительным.

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

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

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


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

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


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

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



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

Реализация метода Якоби на Java:
В качестве $\varepsilon$ берется заранее вычисленное так называемое машинное эпсилон.

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

}

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

Из исходной системы снова выразим $x$, но расставим немного иначе счетчики и добавим параметр $\omega$:

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


При $\omega=1$ это все превращается в метод Якоби.

Итак, будем искать какое-нибудь «хорошее» $\omega$. Зададим какое-нибудь число $s$ и будем беребирать значения $\omega \in (0,2)$, для каждого из которых будем считать нормы $||x^{k+1}-x^k||, \ k = \overline{1,s}$. Для наименьшей из этих норм запомним данное значение $\omega$, и с его помощью будем решать нашу систему.

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

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

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


Существует еще масса алгоритмов для решения СЛАУ. Например, метод квадратного корня, в котором искомая система заменяется на две «простых» системы, решения которых вычисляется по элементарным формулам; метод прогонки, который используется для так специфических трехдиагональных матриц. Каждый сам выбирает, какой метод ему использовать для своей проблемы.
Теги:
Хабы:
+11
Комментарии5
-1
Карма
0
Рейтинг
@inchestnov

Пользователь

Отправить сообщение

Публикации

Показать лучшие за всё время
Изменить настройки темы

Истории

Ближайшие события

Серия занятий «Тренировки по алгоритмам 5.0» от Яндекса
Дата1 марта – 19 апреля
Время19:00
Место
Онлайн
Подробнее в календаре
Дизайнеры, выбирайте себе компанию по вайбам на Хабр Карьере
Дата15 – 28 апреля
Место
Онлайн
Подробнее в календаре
Конференция по анализу данных и технологиям ИИ «Data Fusion»
Дата17 – 18 апреля
Время11:00 – 20:00
Место
Москва
Подробнее в календаре
Бесплатный вебинар по 3D-моделированию катаны из Cyberpunk 2077
Дата18 апреля
Время20:00
Место
Онлайн
Подробнее в календаре
Weekend Offer в AliExpress
Дата20 – 21 апреля
Время10:00 – 20:00
Место
Онлайн
Подробнее в календаре
Вебинар «Истории Dodo Brands и Offprice: как компании мигрировали 1С в облако»
Дата23 апреля
Время11:00
Место
Онлайн
Подробнее в календаре
Вебинар «Информационная архитектура компании»
Дата23 апреля
Время19:00
Место
Онлайн
Подробнее в календаре
Вебинар «Специфика синхронных и асинхронных инсертов в ClickHouse»
Дата23 апреля
Время20:00
Место
Онлайн
Подробнее в календаре
Firebird Conf: конференция для разработчиков и администраторов СУБД Firebird
Дата6 июня
Время09:00 – 20:00
Место
Москва
Подробнее в календаре

Ваш аккаунт

Разделы

Информация

Услуги

Техническая поддержка
© 2006–2024, Habr