Численные методы решения СЛАУ: Гаусс, LU, Якоби и Зейдель

Расчёт температуры в котле на сетке 1000×1000 даёт миллион линейных уравнений с миллионом неизвестных. Прямой метод Гаусса из линейной алгебры справится? Формально да — около $\frac{2}{3}{n}^3\approx 7\cdot 10^{17}$ арифметических операций. На машине с триллионом операций в секунду это примерно неделя, плюс миллион на миллион хранимых коэффициентов — терабайты памяти. А элементы матрицы при этом почти нулевые: каждая точка сетки связана лишь с четырьмя соседями. Лобовая атака здесь расточительна — нужны другие стратегии.

Почему лобовой Гаусс — не всегда выход#

Метод Гаусса — золотой стандарт для плотных систем, где каждый коэффициент значим: он надёжен (с выбором главного элемента) и предсказуем. Уравнения упругой балки, электросети, разреженной сетки из задачи теплопроводности устроены иначе: в строке длиной миллион отличны от нуля пять элементов.

Время — только полбеды. Плотная матрица размера $10^6$ требует хранить $10^{12}$ чисел — терабайты, и это при том, что в задаче теплопроводности отличны от нуля примерно пять чисел в строке. Вся конструкция держится на одной мысли: и время, и память тратятся на нули. Отсюда два направления атаки — умные прямые методы, которые учитывают структуру заполнения, и итерационные, которые нулей вообще не касаются.

Выбор главного элемента: спасение от катастрофы#

Решим руками маленькую систему — и устроим ей маленькую ручную катастрофу:

Гаусс без затей исключает $x$ из второй строки: ведущий элемент $0,0001$, множитель $m=10000$. Вторая строка превращается в $-9999y=-9998$. Если машина хранит лишь три значащие цифры, оба числа округляются до $-1,00\cdot 10^4$, и выходит $y=1,00$. Обратный ход: $0,0001x=1-1,00=0$, откуда $x=0$. Не «немного неточно» — совсем другой ответ: ноль вместо единицы, стопроцентная ошибка.

1. 1
   Смотрим на первый столбец: $|1|>|0,0001|$ — меняем строки местами, ведущим становится элемент $1$
2. 2
   Исключаем $x$: множитель $m=0,0001$; вторая строка: $(1-0,0001)y=1-0,0002$, то есть $0,9999y=0,9998$
3. 3
   Округления до трёх цифр почти ничего не портят: $y\approx 1,00$, обратный ход даёт $x=2-1,00=1,00$
4. 4
   Ошибка порядка $0,01\%$ — три верные цифры вместо откровенного мусора

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

LU-разложение: одна факторизация на много правых частей#

Гаусс делает две вещи одновременно: приводит матрицу к треугольному виду и правую часть — вместе с ней. Если их разделить, получится LU-разложение (что это даёт при многих правых частях и как считать L и U — отдельная словарная статья): $A=LU$, где $L$ — нижнетреугольная с единицами на диагонали, $U$ — верхнетреугольная. Решение $Ax=b$ распадается на два треугольных:

Зачем разделять? Затем, что факторизация — дорогой этап порядка $n^3$, а оба треугольных хода — дешёвые, порядка $n^2$. Если правых частей много (та же матрица жёсткости при разных нагрузках; нормальные уравнения МНК при переборе моделей), матрицу факторизуют один раз, а дальше каждая новая правая часть обходится в копейки. Мини-пример: $A=\left(\begin{array}{cc}2& 1\\ 4& 3\end{array}\right)=\left(\begin{array}{cc}1& 0\\ 2& 1\end{array}\right)\left(\begin{array}{cc}2& 1\\ 0& 1\end{array}\right)$. Для $b=(5;11)$: из $Ly=b$ получается $y=(5;1)$, из $Ux=y$ — $x=(2;1)$. Подстановка в исходную систему подтверждает.

Число обусловленности: когда виновата не программа#

Вот система, которая сведёт с ума любой алгоритм:

С правой частью $b=(2;2)$ решение — $x=(2;0)$. Сдвинем второй компонент на $0,0001$: при $b^{\prime }=(2;2,0001)$ решение ускакивает в $x^{\prime }=(1;1)$. Изменение данных на пять стотысячных — изменения решения в разы. Это свойство самой задачи, и называется оно плохой обусловленностью. Мерой служит число обусловленности $\kappa (A)=\Vert A\Vert \cdot \Vert A^{-1}\Vert$: во сколько раз относительная ошибка решения может превысить относительную ошибку данных. Здесь $\kappa \approx 4\cdot 10^4$ — теряются около четырёх цифр точности, что бы вы ни программировали.

Полезная модель: $\kappa =10^k$ — считайте, что $k$ младших цифр результата ненадёжны. Классический пациент — матрица Гильберта с элементами $\frac{1}{i+j-1}$: уже при $n=10$ обусловленность порядка $10^{13}$, и даже точный алгоритм в двойной точности сохранит три-четыре верные цифры. Откуда такие матрицы берутся? Из постановок задач: регрессия со степенями $x$ даёт матрицу Вандермонда, обработка сигналов — матрицы Тёплица. Прежде чем обвинять программу, посчитайте $\kappa$ — или хотя бы чуть шевелите правую часть и смотрите, как пляшет ответ: если от пыли в пятом знаке решение меняется в разы, вы имеете дело с обусловленностью, а не с багом. Спасает тут не алгоритм, а пересмотр постановки: больше данных, регуляризация, другая параметризация. Формальное определение через нормы матриц даёт статья число обусловленности.

Итерации Якоби и Зейделя: как это работает#

Идея итерационных методов: угадать начальное приближение и уточнять его по простому правилу. Разрешим каждое уравнение относительно «своей» диагональной переменной. Для системы $10x+y=12$, $x+10y=21$ (решение $(1;2)$) получаем:

Стартуем с нуля. Первая итерация: $x^{(1)}=1,2$, $y^{(1)}=2,1$. Вторая: $x^{(2)}=(12-2,1)/10=0,99$, $y^{(2)}=(21-1,2)/10=1,98$. Значения прыгают вокруг решения, но амплитуда падает — итерации сходятся. Метод Зейделя отличается одной деталью: свежие компоненты подставляются немедленно, не дожидаясь конца итерации:

Сходятся ли итерации? Точный критерий — спектральный радиус матрицы перехода меньше единицы; проверять его в лоб тяжело, поэтому пользуются достаточным условием: строгое диагональное преобладание $|a_{ii}|>{\sum }_{j\ne i}|a_{ij}|$ для каждой строки. В нашем примере $10>1$ — выполняется с запасом, и оба метода сходятся.

Разгадка ускорения — в матрице перехода. У Якоби её спектральный радиус здесь равен $0,1$: каждая итерация сокращает ошибку примерно в десять раз. У Зейделя радиус равен квадрату якобиевского — $0,01$, и ошибка падает в сто раз за итерацию. Числа из таблицы подтверждают: после третьей итерации Якоби ошибся на $2\cdot 10^{-3}$, Зейдель — на $2\cdot 10^{-5}$. Общий закон: чем больше диагональный элемент на фоне остальных в строке, тем меньше радиус и быстрее сходимость.

Практическое замечание: диагональное преобладание не экзотика, а норма жизни для сеточных задач. Уравнение теплопроводности связывает температуру каждого узла с четырьмя соседями, и её диагональный коэффициент по модулю не меньше суммы влияний соседей. Проверка занимает одну строку кода: обойти строки, сравнить модули. Если преобладания нет, но матрица симметричная и положительно определённая, работает метод сопряжённых градиентов — прямой наследник идей этого урока.

Прямая или итерационная стратегия: что выбирать#

• Маленькая плотная система — Гаусс с главным элементом, без вариантов
• Много правых частей при одной матрице — LU-разложение: факторизация один раз, дальше по $n^2$ на правую часть
• Огромная разреженная система — Зейдель и его старшие братья (сопряжённые градиенты, многосеточные)
• Плохая обусловленность — враг любой стратегии: проверяйте $\kappa$, а не только код

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