Область устойчивости численных методов

Чтобы сравнивать численные методы честно, придумали единый язык — линейный тест и фактор усиления. На нём видно, кто из методов «терпеливый», а кто разваливается при первом же крупном шаге.

Область устойчивости метода — множество комплексных значений z = h*λ, при которых численное решение линейного теста y' = λ*y не растёт по модулю; иными словами, множество таких z, для которых фактор усиления удовлетворяет условию |R(z)| ≤ 1.

Зачем нужен общий язык

В прошлом уроке мы вручную выводили условие h < 0.04 для явного Эйлера на конкретном уравнении. Но методов много, уравнений ещё больше, и каждый раз заново выводить границу шага неудобно. Нужен универсальный измеритель устойчивости, не зависящий ни от конкретного λ, ни от конкретного h по отдельности, а только от их произведения. Такой измеритель — фактор усиления R(z). Он позволяет нарисовать для каждого метода его «область устойчивости» на комплексной плоскости и потом, зная собственные числа задачи, мгновенно прикинуть допустимый шаг. Это рабочий инструмент: именно так в учебниках и документации сравнивают RK4, Эйлеры и BDF.

Линейный тест Далквиста

Договоримся проверять все методы на одном эталонном уравнении:

y' = λ*y,   λ — комплексное число
точное решение: y(t) = y(0) * exp(λ*t)

Почему именно оно? Любую линейную систему y' = A*y заменой переменных приводят к набору независимых уравнений вида y_i' = λ_i*y_i, где λ_i — собственные числа матрицы A. Поэтому поведение метода на простом скалярном тесте с комплексным λ полностью описывает его поведение на линейной системе. А нелинейные задачи вблизи положения равновесия линеаризуются — и опять сводятся к этому же тесту. Если Re(λ) < 0, точное решение затухает, и мы хотим, чтобы численное тоже затухало.

Фактор усиления R(z)

Применим метод к тесту за один шаг и выразим y_{n+1} через y_n. Для любого одношагового метода получится y_{n+1} = R(z) * y_n, где z = h*λ — безразмерная комбинация, а R(z) — рациональная функция, своя у каждого метода. Число R(z) показывает, во сколько раз меняется решение за шаг. Условие устойчивости простое: |R(z)| ≤ 1. Если оно выполнено, численное решение не растёт; если |R(z)| > 1 — взрывается. Выведем R для трёх ключевых методов:

Явный Эйлер:    y_{n+1} = y_n + h*λ*y_n        => R(z) = 1 + z
Неявный Эйлер:  y_{n+1} = y_n + h*λ*y_{n+1}    => R(z) = 1 / (1 - z)
RK4:            R(z) = 1 + z + z^2/2 + z^3/6 + z^4/24

Формула RK4 — это первые пять членов разложения exp(z), что не случайно: точный множитель за шаг равен exp(z), и чем больше членов ряда метод воспроизводит, тем он точнее.

Как работает под капотом

Разберём каждый метод по его R(z) для вещественного отрицательного λ (тогда z = h*λ — отрицательное вещественное число).

Явный Эйлер, R = 1 + z. Условие |1 + z| < 1 на вещественной оси даёт -2 < z < 0, то есть h*|λ| < 2, или h < 2/|λ|. Область устойчивости на комплексной плоскости — круг радиуса 1 с центром в точке -1. Маленькая область — отсюда жёсткие проблемы.

Неявный Эйлер, R = 1/(1 - z). Условие |1/(1 - z)| < 1 равносильно |1 - z| > 1. Для любого z с Re(z) < 0 точка 1 - z лежит правее единицы по вещественной части, и её модуль больше 1. Значит, неявный Эйлер устойчив при ВСЕХ z с Re(z) < 0 — вся левая полуплоскость целиком. Это свойство называют A-устойчивостью, и именно оно делает метод пригодным для жёстких задач: шаг можно брать любым, взрыва не будет.

RK4. Его область устойчивости заметно больше, чем у явного Эйлера: на вещественной оси она тянется примерно до z = -2.785. Но это всё равно ограниченная область — RK4 не A-устойчив и на жёстких задачах тоже спотыкается, просто чуть позже.

Грубая карта областей устойчивости (вещественная ось z):

  z:  -3.0  -2.785  -2.0   -1.0    0.0
       |      |      |      |      |
Явн.Эйлер:           [############]        (устойчив на -2 .. 0)
RK4:        [###################]           (устойчив ~ -2.785 .. 0)
Неяв.Эйлер: ##########################...  (устойчив на всей левой полуоси)

Проверим численно фактор усиления явного Эйлера и найдём его границу. Берём вещественные z и считаем |1 + z|.

def R_explicit_euler(z):
    return 1.0 + z

print("Явный Эйлер: |R(z)| = |1 + z| для вещественных z")
for z in [-0.5, -1.0, -1.5, -2.0, -2.5]:
    r = R_explicit_euler(z)
    status = "устойчив" if abs(r) <= 1 else "растёт"
    print(f"  z={z:+.1f}: R={r:+.2f}, |R|={abs(r):.2f} -> {status}")

print("граница: |1+z|=1 при z=-2, то есть h = 2/|λ|")

Вывод:

Явный Эйлер: |R(z)| = |1 + z| для вещественных z
  z=-0.5: R=+0.50, |R|=0.50 -> устойчив
  z=-1.0: R=+0.00, |R|=0.00 -> устойчив
  z=-1.5: R=-0.50, |R|=0.50 -> устойчив
  z=-2.0: R=-1.00, |R|=1.00 -> устойчив
  z=-2.5: R=-1.50, |R|=1.50 -> растёт

Видно, как при z = -1 множитель ровно ноль (за один шаг решение «садится» в ноль — это особенность Эйлера в этой точке), а после z = -2 модуль превышает единицу и начинается рост. Теперь сравним все три метода в одной точке, где явный Эйлер уже сломался, а неявный держится. Возьмём z = -3 (например, λ = -30, h = 0.1). Для комплексных z пригодится модуль через abs() — Python считает его и для комплексных, и для вещественных чисел.

def R_explicit(z):
    return 1 + z

def R_implicit(z):
    return 1 / (1 - z)

def R_rk4(z):
    return 1 + z + z**2/2 + z**3/6 + z**4/24

for z in [-1.0, -2.0, -3.0]:
    re = abs(R_explicit(z))
    ri = abs(R_implicit(z))
    rk = abs(R_rk4(z))
    print(f"z={z:+.1f}: |R_явн|={re:.3f}  |R_RK4|={rk:.3f}  |R_неявн|={ri:.3f}")

Вывод:

z=-1.0: |R_явн|=0.000  |R_RK4|=0.375  |R_неявн|=0.500
z=-2.0: |R_явн|=1.000  |R_RK4|=0.333  |R_неявн|=0.333
z=-3.0: |R_явн|=2.000  |R_RK4|=0.375  |R_неявн|=0.250

При z = -3 явный Эйлер уже взрывается (|R| = 2), а и RK4, и неявный Эйлер ещё спокойно затухают. Но если двинуться к z = -4, RK4 даст |R| ≈ 1.1 и тоже сломается, тогда как неявный Эйлер продолжит давать |R| < 1 для сколь угодно больших по модулю z. В этом и состоит качественное превосходство неявных схем на жёстких задачах.

А-устойчивость словами

Метод называют A-устойчивым, если его область устойчивости содержит всю левую полуплоскость Re(z) < 0. Смысл прост: для любого затухающего точного решения (а Re(λ) < 0 — это и есть затухание) численное решение тоже затухает при любом размере шага. Никакой связи между шагом и собственными числами больше нет — вы выбираете h только из соображений точности, а не выживания. Барьер Далквиста гласит, что среди многошаговых методов A-устойчивыми могут быть лишь методы не выше второго порядка — поэтому в реальных солверах используют не чистую A-устойчивость, а более мягкое свойство жёсткой устойчивости, о котором речь в следующем уроке.

Частые ошибки

Сравнивать методы по порядку точности, забыв про устойчивость. RK4 точнее Эйлера, но на жёсткой задаче его маленькая область устойчивости сделает его бесполезным при крупном шаге. Точность и устойчивость — две независимые оси качества.
Считать z по отдельным λ и h. Устойчивость зависит только от произведения z = h*λ. Один и тот же метод устойчив или нет в зависимости от этой комбинации, а не от λ или h по отдельности.
Игнорировать комплексные λ. У систем с колебаниями собственные числа комплексные, и проверять нужно |R(z)| для комплексного z, а не только на вещественной оси. Модуль abs() в Python работает и там.
Путать |R(z)| < 1 и |R(z)| <= 1. На самой границе (|R| = 1) решение не растёт, но и не затухает — оно «висит». Для затухания нужно строгое неравенство; на границе метод нейтрально устойчив.

Итоги

Линейный тест y' = λ*y — универсальный полигон: к нему сводится любая линейная система через собственные числа.
За один шаг любой одношаговый метод даёт y_{n+1} = R(z)*y_n, где z = h*λ; устойчивость это |R(z)| ≤ 1.
Явный Эйлер: R = 1 + z, маленький круг устойчивости, h < 2/|λ| для вещественного λ.
RK4: R — частичная сумма ряда exp(z), область больше (до z ≈ -2.785), но всё ещё ограничена.
Неявный Эйлер: R = 1/(1 - z), устойчив для всех z с Re(z) < 0 — это A-устойчивость, ключ к жёстким задачам.