Неявные методы и BDF для жёстких задач

Мы знаем диагноз — жёсткость — и знаем лекарство — A-устойчивость. Осталось понять, какие методы её обеспечивают на практике, чем за неё платят и как выбрать солвер в реальном коде.

BDF (Backward Differentiation Formulas, формулы обратного дифференцирования) — семейство многошаговых неявных методов, аппроксимирующих производную по нескольким предыдущим точкам и текущей (ещё неизвестной); они жёстко-устойчивы и лежат в основе солверов scipy 'BDF' и (через близкие идеи) 'Radau'.

Зачем нужны неявные методы

Из предыдущего урока ясно: на жёсткой задаче явный метод вынужден держать шаг по самой быстрой моде, и это убивает производительность. Неявные методы вроде неявного Эйлера A-устойчивы — их шаг ограничен только требуемой точностью, а не устойчивостью. Можно идти крупным шагом по гладкой медленной динамике, не боясь взрыва от давно затухшей быстрой моды. Это превращает «5000 шагов» из первого урока в «несколько десятков». Ради этого выигрыша инженеры и терпят главное неудобство неявных методов — необходимость решать уравнение на каждом шаге.

В чём «неявность» и её цена

Сравните две формулы Эйлера:

Явный:    y_{n+1} = y_n + h * f(t_n,     y_n)        # правая часть готова
Неявный:  y_{n+1} = y_n + h * f(t_{n+1}, y_{n+1})    # y_{n+1} стоит в обеих частях

В явной формуле неизвестное y_{n+1} выражено напрямую — подставил и посчитал. В неявной оно входит и в левую, и в правую часть под знаком функции f. Это уравнение относительно y_{n+1}, которое в общем случае нелинейно и аналитически не решается. На каждом шаге приходится решать его численно — обычно методом Ньютона, а для систем — методом Ньютона с матрицей Якоби, то есть решать ещё и СЛАУ внутри каждой итерации. Вот цена A-устойчивости: один шаг неявного метода дороже шага явного в разы. Но шагов нужно несравнимо меньше, и на жёсткой задаче неявный метод выигрывает с огромным запасом.

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

Для линейного теста y' = λ*y неявность снимается алгеброй, и метод Ньютона не нужен. Подставим f = λ*y в неявный Эйлер: y_{n+1} = y_n + h*λ*y_{n+1}. Соберём y_{n+1} в одну сторону: y_{n+1}*(1 - h*λ) = y_n, откуда y_{n+1} = y_n / (1 - h*λ). Множитель за шаг — ровно тот самый R(z) = 1/(1 - z) из прошлого урока. Прогоним неявный Эйлер на жёстком y' = -50*y крупным шагом h = 0.05, при котором явный Эйлер взрывался, и сравним их бок о бок.

import math

def explicit_euler(h, steps, lam=-50.0, y0=1.0):
    y = y0
    out = [y]
    for _ in range(steps):
        y = y + h * (lam * y)
        out.append(y)
    return out

def implicit_euler(h, steps, lam=-50.0, y0=1.0):
    y = y0
    out = [y]
    factor = 1.0 / (1.0 - h * lam)   # для линейного теста неявность решается явно
    for _ in range(steps):
        y = y * factor
        out.append(y)
    return out

h = 0.05
exp_vals = explicit_euler(h, 5)
imp_vals = implicit_euler(h, 5)
print(f"шаг h={h}, множитель неявного = 1/(1-50h) = {1/(1 - h*-50):.4f}")
print("  n |   явный    |  неявный  |  точное exp(-50*t)")
for n in range(6):
    t = n * h
    exact = math.exp(-50 * t)
    print(f"  {n} | {exp_vals[n]:+9.4f} | {imp_vals[n]:8.4f} | {exact:.4f}")

Вывод:

шаг h=0.05, множитель неявного = 1/(1-50h) = 0.2857
  n |   явный    |  неявный  |  точное exp(-50*t)
  0 |   +1.0000 |   1.0000 | 1.0000
  1 |   -1.5000 |   0.2857 | 0.0821
  2 |   +2.2500 |   0.0816 | 0.0067
  3 |   -3.3750 |   0.0233 | 0.0006
  4 |   +5.0625 |   0.0067 | 0.0000
  5 |   -7.5938 |   0.0019 | 0.0000

Картина исчерпывающая. Явный Эйлер при h = 0.05 расходится — значения растут и скачут по знаку. Неявный Эйлер при том же шаге устойчиво затухает: множитель 0.2857 меньше единицы, и решение монотонно идёт к нулю, повторяя качественное поведение точного exp(-50*t). Точность при крупном шаге, конечно, грубая (неявный «садится» к нулю медленнее точного), но это вопрос точности, а не выживания — и решается он уменьшением шага или методом более высокого порядка. Главное: неявный метод не взрывается там, где явный безнадёжен.

Что такое BDF

Неявный Эйлер — это BDF первого порядка (BDF1). Идея BDF: аппроксимировать производную y' в точке t_{n+1} не вперёд, а назад — через несколько уже известных предыдущих значений плюс искомое y_{n+1}. Чем больше точек используется, тем выше порядок. Формулы для чтения:

BDF1 (неявный Эйлер):  (y_{n+1} - y_n) / h = f(t_{n+1}, y_{n+1})
BDF2:  (3*y_{n+1} - 4*y_n + y_{n-1}) / (2*h) = f(t_{n+1}, y_{n+1})
BDF3:  (11*y_{n+1} - 18*y_n + 9*y_{n-1} - 2*y_{n-2}) / (6*h) = f(t_{n+1}, y_{n+1})

Правая часть берётся в НОВОЙ точке t_{n+1} => метод неявный.

BDF1 и BDF2 A-устойчивы; BDF порядков 3–6 уже не A-устойчивы (барьер Далквиста не позволяет), но обладают более слабым, зато достаточным на практике свойством — жёсткой устойчивостью: они устойчивы в той части левой полуплоскости, которая важна для жёстких задач (большие отрицательные вещественные части), и лишь в узком секторе у мнимой оси теряют устойчивость. Этого достаточно для подавляющего большинства жёстких систем, поэтому BDF до 5-го порядка — рабочая лошадка промышленных солверов (легендарный код DASSL, scipy.integrate с method='BDF').

Как выбрать метод на практике

В реальном Python вы почти никогда не пишете солвер сами — берёте scipy.integrate.solve_ivp и выбираете метод параметром. Правило выбора простое:

from scipy.integrate import solve_ivp

# НЕЖЁСТКАЯ задача (плавная динамика, один масштаб времени):
sol = solve_ivp(f, [0, 10], y0, method='RK45')   # Dormand-Prince, по умолчанию

# ЖЁСТКАЯ задача (разные масштабы, явные методы взрываются):
sol = solve_ivp(f, [0, 10], y0, method='BDF')     # многошаговый неявный
sol = solve_ivp(f, [0, 10], y0, method='Radau')   # неявный Рунге-Кутты, очень устойчив
sol = solve_ivp(f, [0, 10], y0, method='LSODA')   # САМ переключается жёсткий/нежёсткий

# подсказать якобиан ускоряет неявные методы:
sol = solve_ivp(f, [0, 10], y0, method='BDF', jac=jacobian)

Логика такая. Если задача нежёсткая — RK45 (метод Дормана-Принса, адаптивный явный) даёт лучшее соотношение точности и скорости. Если жёсткая — нужен неявный: BDF (быстрый, многошаговый), Radau (очень устойчивый неявный Рунге-Кутты, хорош для высокой точности и сложных систем) или LSODA, который автоматически детектирует жёсткость по ходу интегрирования и сам переключается между Adams (нежёсткий режим) и BDF (жёсткий). Если вы не уверены, жёсткая задача или нет — начните с RK45; если он зависает, делает миллионы шагов или выдаёт предупреждение о малом шаге, переключайтесь на BDF или LSODA. Передача якобиана (jac=) ускоряет неявные методы, потому что избавляет от его численной оценки на каждой ньютоновской итерации.

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

Применять RK45 к жёсткой задаче и удивляться зависанию. Явный адаптивный метод не взорвётся (он уменьшит шаг), но сделает это ценой огромного числа шагов и будет считать вечность. Симптом — солвер «думает» минутами на простой с виду системе.
Забыть, что неявность стоит денег. BDF и Radau решают нелинейное уравнение на каждом шаге методом Ньютона. На нежёсткой задаче это лишние затраты — там RK45 быстрее. Неявный метод оправдан только жёсткостью.
Не передавать якобиан для больших систем. Без jac= scipy оценивает матрицу Якоби численно, тратя много вызовов f. Для систем из сотен уравнений аналитический якобиан ускоряет счёт в разы.
Считать BDF высоких порядков A-устойчивыми. A-устойчивы только BDF1 и BDF2. Порядки 3–6 лишь жёстко-устойчивы — этого хватает для типичных жёстких задач, но не для систем с собственными числами у самой мнимой оси (чисто колебательных).
Брать слишком высокий порядок «для точности». BDF выше 6-го порядка неустойчивы и в солверах не используются. Точность повышают уменьшением шага и контролем ошибки, а не безграничным ростом порядка.

Итоги

Неявные методы A-устойчивы — крупный шаг на жёсткой задаче без взрыва; цена — решение уравнения (метод Ньютона + СЛАУ) на каждом шаге.
Неявный Эйлер на жёстком y' = -50*y при h = 0.05 устойчив (множитель 0.2857), тогда как явный расходится.
BDF — многошаговые неявные формулы обратного дифференцирования; BDF1=неявный Эйлер, BDF1/2 A-устойчивы, BDF3–6 жёстко-устойчивы.
Выбор в scipy: нежёсткое → RK45 (Dormand-Prince); жёсткое → BDF / Radau / LSODA (авто-детект жёсткости).
Передача якобиана (jac=) ускоряет неявные методы; высокий порядок BDF не панацея — выше 6 нестабилен.