Жёсткие (stiff) уравнения: что это

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

Жёсткое (stiff) уравнение — дифференциальное уравнение или система, в которой одновременно присутствуют процессы с резко разными характерными временами: быстрая компонента затухает почти мгновенно, а медленная живёт долго; из-за этого явные методы вынуждены идти шагом, который диктует уже исчезнувшая быстрая мода.

Зачем это знать

Представьте, что вы моделируете химическую реакцию. Один промежуточный продукт распадается за микросекунды, а основное вещество расходуется за часы. Или электрическую цепь, где паразитная ёмкость заряжается за наносекунды, а полезный сигнал меняется за миллисекунды. Или остывание детали, у которой тонкая поверхностная плёнка реагирует мгновенно, а массивное ядро — медленно. Во всех этих случаях вас интересует медленная, «полезная» динамика, и хочется считать её крупным шагом. Но наивный явный метод этого не позволит: он будет спотыкаться о давно затухший быстрый процесс. Если не знать про жёсткость, вы либо получите бессмысленный результат («числа улетели в бесконечность»), либо потратите часы машинного времени там, где грамотный метод справился бы за секунду. Жёсткость — это не баг вашего кода, это свойство самой задачи, и распознавать его нужно заранее.

Разные масштабы времени

Возьмём простейшую линейную систему из двух уравнений:

y1' = -1000 * y1
y2' = -1 * y2

Эти уравнения не связаны друг с другом, поэтому решаются независимо: y1(t) = y1(0) * exp(-1000*t) и y2(t) = y2(0) * exp(-t). Первая компонента затухает с характерным временем 1/1000 = 0.001 — фактически она исчезает за несколько тысячных секунды. Вторая компонента живёт с характерным временем 1 секунда. Отношение этих масштабов — 1000. Это число и называют коэффициентом жёсткости. Чем оно больше, тем жёстче задача.

В матричной форме система записывается как y' = A*y, где A — диагональная матрица с числами -1000 и -1 на диагонали. Числа на диагонали (а в общем случае — собственные числа матрицы A) и задают скорости затухания мод. Когда среди собственных чисел есть и очень большое по модулю отрицательное (-1000), и умеренное (-1), мы имеем дело с жёсткой системой. Подчеркнём: дело именно в разбросе модулей, а не в их абсолютной величине.

Парадокс мёртвой моды

Самое неинтуитивное здесь вот что. Через сотую долю секунды быстрая компонента y1 практически равна нулю — она «умерла». Казалось бы, дальше можно считать только медленную часть крупным шагом, ведь быстрого процесса больше нет. Но явный метод так не умеет. Его устойчивость определяется самым большим по модулю собственным числом всей системы — то есть тем самым -1000, даже когда соответствующая ему компонента давно исчезла из решения. Метод вынужден держать шаг достаточно мелким, чтобы «угнаться» за призраком быстрой моды. Отсюда и название явления: задача требует абсурдно мелкого шага не из-за того, что в решении происходит что-то интересное, а из-за невидимого, уже затухшего процесса. Это ловушка, в которую попадают все явные схемы.

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

Почему явный метод чувствует мёртвую моду? Разберём на явном Эйлере для скалярного уравнения y' = -50*y. Шаг Эйлера: y_{n+1} = y_n + h * f(y_n) = y_n + h*(-50*y_n) = (1 - 50*h) * y_n. На каждом шаге значение умножается на коэффициент 1 - 50*h. Если |1 - 50*h| < 1, последовательность затухает (это правильно: точное решение тоже затухает). Если же |1 - 50*h| > 1, значения растут по модулю — численное решение взрывается, хотя точное стремится к нулю.

Найдём границу. Неравенство |1 - 50*h| < 1 означает -1 < 1 - 50*h < 1, то есть 0 < 50*h < 2, откуда h < 2/50 = 0.04. Итог: при h < 0.04 явный Эйлер устойчив, при h > 0.04 — расходится. Заметьте, что само точное решение exp(-50*t) ведёт себя идеально гладко при любом t; ограничение на шаг — целиком артефакт численной схемы. В системе с собственными числами -1000 и -1 быстрая мода даст условие h < 2/1000 = 0.002, и именно она продиктует шаг для всей системы, даже когда давно равна нулю.

def euler_step_factor(h, lam):
    # множитель явного Эйлера для y' = lam*y
    return 1.0 + h * lam

lam = -50.0
for h in [0.05, 0.04, 0.03, 0.01]:
    r = euler_step_factor(h, lam)
    status = "устойчив" if abs(r) < 1 else "ВЗРЫВ"
    print(f"h={h:.2f}: множитель 1-50h = {r:+.3f}, |множитель|={abs(r):.3f} -> {status}")

print("граница устойчивости: h = 2/50 =", 2/50)

Вывод:

h=0.05: множитель 1-50h = -1.500, |множитель|=1.500 -> ВЗРЫВ
h=0.04: множитель 1-50h = -1.000, |множитель|=1.000 -> ВЗРЫВ
h=0.03: множитель 1-50h = -0.500, |множитель|=0.500 -> устойчив
h=0.01: множитель 1-50h = +0.500, |множитель|=0.500 -> устойчив

Теперь прогоним сам метод на нескольких шагах, чтобы увидеть взрыв и затухание своими глазами. Возьмём y(0) = 1 и проследим первые шаги при h = 0.05 (взрыв) и h = 0.01 (устойчиво).

def run_euler(h, steps):
    lam = -50.0
    y = 1.0
    values = [y]
    for _ in range(steps):
        y = y + h * (lam * y)   # явный Эйлер
        values.append(y)
    return values

print("h=0.05 (h > 0.04, ожидаем взрыв):")
for i, v in enumerate(run_euler(0.05, 6)):
    print(f"  шаг {i}: y = {v:+.4f}")

print("h=0.01 (h < 0.04, ожидаем затухание):")
for i, v in enumerate(run_euler(0.01, 6)):
    print(f"  шаг {i}: y = {v:+.4f}")

Вывод:

h=0.05 (h > 0.04, ожидаем взрыв):
  шаг 0: y = +1.0000
  шаг 1: y = -1.5000
  шаг 2: y = +2.2500
  шаг 3: y = -3.3750
  шаг 4: y = +5.0625
  шаг 5: y = -7.5938
  шаг 6: y = +11.3906
h=0.01 (h < 0.04, ожидаем затухание):
  шаг 0: y = +1.0000
  шаг 1: y = +0.5000
  шаг 2: y = +0.2500
  шаг 3: y = +0.1250
  шаг 4: y = +0.0625
  шаг 5: y = +0.0312
  шаг 6: y = +0.0156

При h = 0.05 значения не просто растут, а ещё и меняют знак на каждом шаге — это классическая картина численной неустойчивости: осцилляции с растущей амплитудой. Точное решение exp(-50*t) при этом монотонно и спокойно стремится к нулю. При h = 0.01 численное решение ведёт себя разумно: каждый шаг уменьшает значение вдвое (множитель 0.5), что качественно верно.

Жёсткость — это про систему, а не про одно число

В скалярном случае y' = -50*y нет двух разных масштабов, и формально это ещё не жёсткая задача — просто пример, на котором удобно увидеть ограничение шага. Настоящая жёсткость появляется, когда в одной системе сосуществуют быстрые и медленные моды. Тогда возникает конфликт: интервал интегрирования (например, до t = 10) определяется медленной модой, а максимально допустимый шаг — быстрой. Число шагов получается как отношение этих величин, и для системы с числами -1000 и -1 на интервале до t = 10 явному Эйлеру понадобится порядка 10 / 0.002 = 5000 шагов, хотя «полезная» динамика плавная. Увеличишь жёсткость до миллиона — понадобятся миллионы шагов. Это и есть практическая боль, ради решения которой придумали неявные методы.

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

Считать, что взрыв — ошибка в коде. Если решение «улетает в бесконечность», первая мысль — опечатка. Но часто это корректно реализованный явный метод на жёсткой задаче с слишком крупным шагом. Лечится не отладкой, а сменой метода или уменьшением шага.
Путать жёсткость с большим собственным числом. Жёсткость — это разброс масштабов. Система только с числом -1000 (без медленной моды) не жёсткая в практическом смысле: вам и так нужен мелкий шаг, потому что динамика быстрая. Проблема возникает именно от соседства -1000 и -1.
Думать, что мелкий шаг всегда спасает. Технически да, но ценой огромного числа шагов и накопления ошибок округления. Уменьшать шаг до 0.002 на интервале до 10 — это 5000 итераций ради гладкой кривой. Неэкономно.
Ожидать, что затухшая мода «отпустит» метод. Интуиция подсказывает: раз быстрая компонента стала нулём, ограничение снимется. Нет. Устойчивость явного метода завязана на собственное число, а не на текущее значение компоненты.

Итоги

Жёсткость — это сосуществование процессов с резко разными характерными временами в одной задаче.
Признак: большой разброс модулей собственных чисел (например, -1000 и -1, коэффициент жёсткости ~1000).
Явные методы вынуждены держать шаг по самой быстрой моде, даже когда она уже затухла — отсюда абсурдно мелкий шаг.
Для y' = -50*y явный Эйлер устойчив только при h < 2/50 = 0.04; иначе решение растёт и осциллирует.
Уменьшение шага «работает», но неэкономно — правильный ответ для жёстких задач это неявные методы.