Метод Гаусса
Метод Гаусса -- классический прямой алгоритм решения систем линейных алгебраических уравнений (СЛАУ) -- набора уравнений вида $Ax = b$, где $A$ -- матрица коэффициентов, $x$ -- вектор неизвестных, $b$ -- вектор правой части. В этом примере реализован метод Гаусса с частичным выбором ведущего элемента и нормировкой строки, что повышает численную устойчивость.
Задача
Найти вектор $x$, удовлетворяющий системе линейных алгебраических уравнений (СЛАУ):
$Ax = b$
В примере используется матрица Гильберта размера $n \times n$:
$A_{i,j} = \frac{1}{i + j + 1}, \quad i,j = 0, \ldots, n-1$
и вектор правой части $b_i = 1$. Матрица Гильберта -- классический пример плохо обусловленной матрицы: при росте $n$ число обусловленности растёт экспоненциально, что делает задачу хорошим тестом для численных методов.
Идея алгоритма
- Прямой ход -- последовательно для каждого столбца $k = 0, \ldots, n-1$:
- Найти строку с наибольшим по модулю элементом в столбце $k$ (частичный выбор) и поменять её местами с текущей строкой $k$.
- Проверить, не является ли ведущий элемент слишком малым (относительная проверка на вырожденность).
- Поделить всю строку $k$ на ведущий элемент $A_{k,k}$ (нормировка).
- Вычесть нормированную строку $k$ из всех строк ниже, обнуляя элементы в столбце $k$.
- Обратный ход -- начиная с последней строки, вычислить $x_i$ по формуле:
$x_i = \frac{b_i - \sum_{j=i+1}^{n-1} A_{i,j} \cdot x_j}{A_{i,i}}$
- Проверка -- вычислить невязку $r = A_{\text{orig}} \cdot x - b$ и вывести её нормы $|r|2$ и $|r|\infty$.
Разбор
Параметры и данные
Задаём размер системы $n = 256$ и порог $\varepsilon$ для определения вырожденности. Массив A хранит матрицу, массив B -- правую часть (в него же будет записано решение):
| параметры
цел n
вещ Eps
n := 256
Eps := 0.000000000001 | относительный порог по масштабу строки
| данные
вещ таб A[0:n-1,0:n-1]
вещ таб B[0:n-1]
Инициализация матрицы Гильберта
Заполняем матрицу по формуле $A_{i,j} = \frac{1}{i+j+1}$ и правую часть $b_i = 1$:
| инициализация: A(i,j) = 1/(i+j+1), B(i) = 1
нц для i от 0 до n-1
нц для j от 0 до n-1
A[i,j] := 1.0 / (i + j + 1.0)
кц
B[i] := 1.0
кц
Прямой ход: частичный выбор и нормировка
На каждом шаге $k$ ищем строку с максимальным $|A_{i,k}|$ ниже текущей и меняем строки местами. Затем проверяем ведущий элемент: если $|A_{k,k}| \le \varepsilon \cdot \max_j |A_{k,j}|$, значит строка почти нулевая -- прекращаем. Иначе нормируем строку (делим на $A_{k,k}$) и вычитаем из нижних строк:
| прямой ход Гаусса с частичным выбором; НОРМИРУЕМ всегда
нц для k от 0 до n-1
| выбор ведущей строки по максимуму |A[i,k]|
_imax := k
bestSq := A[k,k] * A[k,k]
нц для i от k+1 до n-1
curSq := A[i,k] * A[i,k]
если curSq > bestSq то
bestSq := curSq
_imax := i
все
кц
| swap k <-> _imax
если _imax <> k то
нц для j от k до n-1
tmp := A[k,j]
A[k,j] := A[_imax,j]
A[_imax,j] := tmp
кц
tmp := B[k]; B[k] := B[_imax]; B[_imax] := tmp
все
p := A[k,k]
| относительная проверка на вырожденность: |p| <= Eps * max_j |A[k,j]|
rowmax := 0.0
нц для j от k до n-1
sum := A[k,j]; если sum < 0.0 то sum := -sum все
если sum > rowmax то rowmax := sum все
кц
scale := Eps * rowmax
sum := p; если sum < 0.0 то sum := -sum все
если sum <= scale то
| почти нулевой пивот — прекращаем исключение
выход
все
| нормируем ведущую строку k
нц для j от k до n-1
A[k,j] := A[k,j] / p
кц
B[k] := B[k] / p
| исключаем ниже
нц для i от k+1 до n-1
tmp := A[i,k]
если tmp <> 0.0 то
нц для j от k до n-1
A[i,j] := A[i,j] - A[k,j] * tmp
кц
B[i] := B[i] - B[k] * tmp
все
кц
кц
Обратный ход
После прямого хода матрица приведена к верхнетреугольному виду с единицами на диагонали (благодаря нормировке). Вычисляем неизвестные от последней строки к первой. Результат записывается в массив B:
| обратный ход
нц для t от 0 до n-1
i := n - 1 - t
sum := B[i]
нц для j от i+1 до n-1
sum := sum - A[i,j] * B[j]
кц
если A[i,i] <> 0.0 то
B[i] := sum / A[i,i]
иначе
B[i] := 0.0
все
кц
Расчёт невязки
Для проверки заново вычисляем $r_i = \sum_j A_{i,j}^{\text{orig}} \cdot x_j - b_i$ (используя исходную формулу матрицы Гильберта) и выводим нормы $|r|2$ и $|r|\infty$:
| === расчёт невязки r = A_orig * X - 1 ===
вещ res_l2_sq, res_l2, res_inf, r, absv
res_l2_sq := 0.0
res_inf := 0.0
нц для i от 0 до n-1
sum := 0.0
нц для j от 0 до n-1
sum := sum + (1.0 / (i + j + 1.0)) * B[j] | A_orig * x
кц
r := sum - 1.0 | вычитаем b_orig
absv := r; если absv < 0.0 то absv := -absv все
res_l2_sq := res_l2_sq + r * r
если absv > res_inf то res_inf := absv все
кц
res_l2 := sqrt(res_l2_sq)
| выводим нормы невязки
вывод "||r||_2^2 = ", res_l2_sq, нс
вывод "||r||_2 = ", res_l2, нс
вывод "||r||_inf = ", res_inf, нс
Сложность алгоритма -- $O(n^3)$ по времени. Частичный выбор ведущего элемента и нормировка строки значительно повышают устойчивость метода к ошибкам округления.
Полная программа
алг main
нач
| Решение Ax=b методом Гаусса с частичным выбором, нормировкой строки
| и обратным ходом с делением на диагональ. Матрица Гильберта, b=1.
| Вектор решения записывается в B. В конце выводятся нормы невязки.
| параметры
цел n
вещ Eps
n := 256
Eps := 0.000000000001 | относительный порог по масштабу строки
| данные
вещ таб A[0:n-1,0:n-1]
вещ таб B[0:n-1]
| служебные
цел i, j, k, _imax, t
вещ p, tmp, bestSq, curSq, sum, rowmax, scale
| инициализация: A(i,j) = 1/(i+j+1), B(i) = 1
нц для i от 0 до n-1
нц для j от 0 до n-1
A[i,j] := 1.0 / (i + j + 1.0)
кц
B[i] := 1.0
кц
| прямой ход Гаусса с частичным выбором; НОРМИРУЕМ всегда
нц для k от 0 до n-1
| выбор ведущей строки по максимуму |A[i,k]|
_imax := k
bestSq := A[k,k] * A[k,k]
нц для i от k+1 до n-1
curSq := A[i,k] * A[i,k]
если curSq > bestSq то
bestSq := curSq
_imax := i
все
кц
| swap k <-> _imax
если _imax <> k то
нц для j от k до n-1
tmp := A[k,j]
A[k,j] := A[_imax,j]
A[_imax,j] := tmp
кц
tmp := B[k]; B[k] := B[_imax]; B[_imax] := tmp
все
p := A[k,k]
| относительная проверка на вырожденность: |p| <= Eps * max_j |A[k,j]|
rowmax := 0.0
нц для j от k до n-1
sum := A[k,j]; если sum < 0.0 то sum := -sum все
если sum > rowmax то rowmax := sum все
кц
scale := Eps * rowmax
sum := p; если sum < 0.0 то sum := -sum все
если sum <= scale то
| почти нулевой пивот — прекращаем исключение
выход
все
| нормируем ведущую строку k
нц для j от k до n-1
A[k,j] := A[k,j] / p
кц
B[k] := B[k] / p
| исключаем ниже
нц для i от k+1 до n-1
tmp := A[i,k]
если tmp <> 0.0 то
нц для j от k до n-1
A[i,j] := A[i,j] - A[k,j] * tmp
кц
B[i] := B[i] - B[k] * tmp
все
кц
кц
| обратный ход
нц для t от 0 до n-1
i := n - 1 - t
sum := B[i]
нц для j от i+1 до n-1
sum := sum - A[i,j] * B[j]
кц
если A[i,i] <> 0.0 то
B[i] := sum / A[i,i]
иначе
B[i] := 0.0
все
кц
| === расчёт невязки r = A_orig * X - 1 ===
вещ res_l2_sq, res_l2, res_inf, r, absv
res_l2_sq := 0.0
res_inf := 0.0
нц для i от 0 до n-1
sum := 0.0
нц для j от 0 до n-1
sum := sum + (1.0 / (i + j + 1.0)) * B[j] | A_orig * x
кц
r := sum - 1.0 | вычитаем b_orig
absv := r; если absv < 0.0 то absv := -absv все
res_l2_sq := res_l2_sq + r * r
если absv > res_inf то res_inf := absv все
кц
res_l2 := sqrt(res_l2_sq)
| выводим нормы невязки
вывод "||r||_2^2 = ", res_l2_sq, нс
вывод "||r||_2 = ", res_l2, нс
вывод "||r||_inf = ", res_inf, нс
кон