Математика / Начала анализа

Метод Ричардсона: стационарная итерация для Ax=b

Метод Ричардсона прибавляет к текущему приближению невязку, умноженную на параметр τ. Для симметричной положительно определённой матрицы разумный τ превращает эту простую формулу в сходящуюся стационарную итерацию.

Опубликовано: Обновлено:

Формула

$$x^{(k+1)}=x^{(k)}+\tau(b-Ax^{(k)}),\qquad r^{(k)}=b-Ax^{(k)}$$
Схема Ричардсон: коррекция по двум невязкам
Ричардсон: коррекция по двум невязкамРичардсон: коррекция по двум невязкамСтарт (0; 0), невязка (2; 8)τ = 1/3 → первый вектор (2/3; 8/3)Новая невязка: (2/3; −8/3)Второй вектор: (8/9; 16/9)

Сравнивайте приближение с точным решением и проверяйте невязку исходной системы.

Обозначения

$τ$
длина шага стационарной коррекции, обратная единица A
$r^(k)$
невязка b−Ax^(k), как у b
$x^(k)$
текущее приближение решения, единица x
$λ_min, λ_max$
крайние собственные значения симметричной положительно определённой матрицы, единица A

Условия применения

Для симметричной положительно определённой A достаточно выбрать 0<τ<2/λ_max. Если известны λ_min и λ_max, классический постоянный выбор τ=2/(λ_min+λ_max) минимизирует на этом интервале наихудший модуль множителя ошибки. Нужны согласованные единицы: τA должно быть безразмерным. Контролируют норму невязки и заранее фиксируют допустимую относительную точность.

Ограничения

Для нессиметричной или неопределённой матрицы простое правило через λ_min и λ_max неприменимо без отдельного анализа. Слишком малый τ даёт медленное движение, слишком большой может вызвать расходимость. Постоянный τ не использует информацию о каждом новом остатке и часто уступает современным методам Крылова. Не следует путать эту итерацию с экстраполяцией Ричардсона для производных: там меняют шаг сетки, здесь исправляют вектор решения.

Подробное объяснение

Невязка r=b−Ax показывает, какие уравнения системы сейчас нарушены. Метод Ричардсона добавляет к вектору x поправку τr. Если A=I и τ=1, одна поправка сразу даёт точное решение; для общей A та же идея требует подобрать масштаб, потому что разные направления могут изменяться по-разному.

Пусть x* — точное решение. Тогда ошибка после шага равна e^(k+1)=(I−τA)e^(k). Для симметричной положительно определённой матрицы собственные направления A не смешиваются, а в направлении собственного значения λ ошибка умножается на 1−τλ. Поэтому условие |1−τλ|<1 для всех λ приводит к 0<τ<2/λ_max.

Выбор τ=2/(λ_min+λ_max) уравнивает по модулю худшие множители на краях спектра. Его фактор сжатия равен (κ−1)/(κ+1), где κ=λ_max/λ_min. Эта формула объясняет медленность плохо обусловленной задачи: если κ велико, множитель близок к единице даже при наилучшем постоянном τ.

Метод особенно полезен как диагностическая и учебная схема. Он связывает невязку, спектральный масштаб и скорость, не скрывая их внутри сложного решателя. Но в реальном расчёте границы спектра могут быть неизвестны, а матрица — не симметричной положительно определённой; тогда нужно применять проверенную для данного класса альтернативу, а не переносить формулу без условий.

Как пользоваться формулой

  1. Убедитесь, что A симметрична и положительно определена либо что для неё есть отдельное доказательство сходимости выбранной итерации.
  2. Оцените λ_max и выберите τ в интервале от нуля до 2/λ_max; при известных двух границах используйте их симметричный выбор.
  3. Вычислите невязку r=b−Ax и прибавьте к текущему x вектор τr, сохраняя одинаковые единицы величин.
  4. После шага пересчитайте норму невязки; по её истории проверяйте сжатие и прекращайте итерацию по заданному допуску.

Историческая справка

В численных расчётах XX–XXI веков имя Льюиса Фрая Ричардсона связано с итерационными идеями в численных расчётах. Современный метод для Ax=b записывается на языке матриц, спектра и норм. Netlib Templates определяет стационарные итерации формой x^(k)=Bx^(k−1)+c; схема Ричардсона реализует этот принцип при B=I−τA и делает параметр шага частью алгоритма.

В XX веке вычислительная практика прошла переход от последовательных приближений к анализу больших разреженных систем, возникающих из дискретных моделей. Фиксированный параметр естественно сравнивают с методами, которые адаптируют направление или длину шага по ходу расчёта. Одноимённая экстраполяция в численном дифференцировании решает другую задачу: она комбинирует результаты на нескольких шагах сетки. Современное употребление имени различает эти две конструкции по объекту вычисления и формуле.

Историческая линия формулы

Принятое имя Richardson iteration относится к коррекции решения линейной системы по невязке. В отличие от экстраполяции Ричардсона, здесь не комбинируют ответы на разных шагах сетки. Стационарный параметр τ определяет матрицу I−τA и тем самым изменение ошибки на каждом проходе.

Пример

Пусть A=diag(2,4), b=(2,8), x^(0)=(0,0). Здесь λ_min=2, λ_max=4, поэтому τ=2/(2+4)=1/3. Невязка r^(0)=(2,8), значит x^(1)=(2/3,8/3). Далее Ax^(1)=(4/3,32/3), r^(1)=(2/3,−8/3), и x^(2)=(8/9,16/9). Точное решение равно (1,2). Ошибка первой компоненты сохраняет знак, а ошибка второй меняет знак: множители равны 1/3 для λ=2 и −1/3 для λ=4. Норма невязки после первого шага не нулевая, поэтому даже диагональная матрица не даёт одношагового решения при общем постоянном τ. Все числа безразмерны; при физических единицах τ должен компенсировать единицу коэффициентов A.

Частая ошибка

Не ставьте τ=1 автоматически: это оправдано только для особого масштаба A. Не измеряйте сходимость одним компонентом x — разные спектральные направления могут вести себя по-разному. Нельзя применять SPD-формулу к нессиметричной матрице только потому, что её диагональ положительна. Наконец, название не делает эту страницу дублем экстраполяции Ричардсона: здесь τ умножает невязку линейной системы, а не комбинирует значения на шагах h и h/2.

Практика

Задачи с решением

Один шаг

Условие. A=2I, b=(6,−2), x^(0)=0, τ=1/2. Найдите x^(1).

Решение. r^(0)=b. Добавляем (1/2)b=(3,−1) к нулевому вектору.

Ответ. (3, −1), точное решение за один шаг.

Допустимый шаг

Условие. Для SPD-матрицы λ_max=5. Какой интервал допустим для τ?

Решение. Требуется 0<τ<2/λ_max=2/5.

Ответ. 0<τ<0,4.

Дополнительные источники

  • https://www.netlib.org/templates/templates.html
  • https://fncbook.com/linear-systems/
  • https://fncbook.com/structure/

Связанные формулы

Математика

Прямой ход метода Гаусса

$R_i\leftarrow R_i-\frac{a_{ik}}{a_{kk}}R_k$

Прямой ход метода Гаусса зануляет коэффициенты под ведущими элементами. В результате система приводится к ступенчатому виду, из которого решение находят обратной подстановкой.

Математика

Обратная подстановка в методе Гаусса

$x_i=\frac{b'_i-\sum_{j=i+1}^{n}u_{ij}x_j}{u_{ii}}$

Обратная подстановка находит неизвестные после прямого хода метода Гаусса. Она идет снизу вверх по ступенчатой системе: сначала последняя ведущая переменная, затем предыдущие.

Математика

Метод Гаусса-Жордана

$\left[A\mid b\right]\sim\left[I\mid x\right]$

Метод Гаусса-Жордана продолжает метод Гаусса до приведенного ступенчатого вида. Если система имеет единственное решение, расширенная матрица превращается в [I|x], и ответ читается сразу.

Математика

Число обусловленности для задачи МНК

$\kappa_2(A^\top A)=\frac{\sigma_{\max}^2}{\sigma_{\min}^2}=\kappa_2(A)^2.$

При переходе к нормальным уравнениям число обусловленности фактически возводится в квадрат, поэтому ошибки округления и шум в данных могут заметно усилиться.