🌊 Диффуры

Как решать жёсткие дифференциальные уравнения

Q: Жёсткое уравнение — это всегда нелинейное?

Нет, наш пример $y' = -1000(y - \cos t)$ линейный. Жёсткость — про разброс временных масштабов, а не про нелинейность: она возникает в линейных системах с сильно разными собственными значениями и в нелинейных задачах, где линеаризация имеет такой же разброс. Нелинейность лишь добавляет проблем неявным методам: на каждом шаге приходится решать уравнение методом Ньютона.

Жёсткие ОДУ: почему явный Эйлер взрывается при h = 0,01, как работает безусловно устойчивый неявный и расчёт для y' = −1000(y − cos t) по шагам.

Сложность: Продвинутый 2 схемы действий

Симуляция химической реакции: константы скоростей различаются в тысячу раз. Вы запускаете численный метод с шагом $h = 0, 01$ — и к пятому шагу значения взлетают до $59050$ , хотя настоящее решение давно вышло на спокойный косинус с амплитудой единица. Уменьшаете шаг в десять раз — счёт растягивается на минуты. Это не глюк программы. Это жёсткость, и у неё есть точный диагноз с лечением.

Жёсткость — свойство задачи, а не метода. Решение содержит компоненты с сильно разными масштабами времени: быстрое затухает за доли секунды, медленное живёт секунды. Формально: собственные значения $λ$ линейной части различаются по модулю в десятки и сотни раз, $∣ λ_{быстр} ∣ ≫ ∣ λ_{медл} ∣$ . Беда в том, что явным методам шаг диктует самый быстрый процесс — даже когда он давно умер и в решении не участвует.

Простейшая модель — распадающаяся система из двух независимых мод: $u_{1}^{'} = - 1000 u_{1}$ и $u_{2}^{'} = - 0, 1 u_{2}$ с точным решением $u_{1} = a e^{- 1000 t}$ , $u_{2} = b e^{- 0, 1 t}$ . Первая мода мертва к $t = 0, 005$ , вторая живёт десятки секунд. Явному методу условие устойчивости запрещает шаг больше $0, 002$ — на интервале в десять секунд это пять тысяч шагов ради отслеживания экспоненты, которая превратилась в ноль на первом шаге. Медленной моде хватило бы шага $h = 1$ : точность разрешает, устойчивость запрещает.

Канонический пример: два масштаба в одном уравнении#

Возьмём задачу, где всё видно руками: $y^{'} = - 1000 (y - cos t)$ , $y (0) = 0$ . Это линейное уравнение первого порядка — приёмы решения стандартные. Ищем частное решение вида $A cos t + B sin t$ и подгоняем коэффициенты: при $sin t$ выходит $A = 1000 B$ , при $cos t$ — $B = 1000 (1 - A)$ , откуда

A = \frac{1 0 ^{6}}{1 0 ^{6} + 1} \approx 0, 999999, B = \frac{1000}{1 0 ^{6} + 1} \approx 0, 001

A и B — коэффициенты установившегося режима: почти чистый косинус с крошечной добавкой синуса

Честная проверка: подстановка «наивного» $y = cos t + C e^{- 1000 t}$ оставляет в уравнении остаток $- sin t \approx - 0, 48$ при $t = 0, 5$ — чистый косинус не решение, нужна поправка $B sin t$ . А вот полное семейство подходит:

y (t) = A cos t + B sin t + (y_{0} - A) e^{- 1000 t}

общее решение: подстановка в уравнение даёт тождество с точностью до 10⁻⁶ — проверено численно

Смотрите, что вышло: переходная добавка $(y_{0} - A) e^{- 1000 t}$ умирает за время порядка $5/1000 = 0, 005$ — уже при $t = 0, 005$ от неё остаётся $e^{- 5} \approx 0, 007$ доли. Дальше решение — медленный косинус. Два масштаба времени $0, 001$ и $1$ , отношение $1000 : 1$ — перед нами жёсткая задача в хрестоматийном виде.

Вывод формулы стоит повторить один раз руками: уравнение приводится к виду $y^{'} + 1000 y = 1000 cos t$ , умножение на интегрирующий множитель $e^{1000 t}$ и интегрирование дают общее решение за пять строк. Проверка подстановкой — вторая строка: производная $y^{'} = - A sin t + B cos t - 1000 (y_{0} - A) e^{- 1000 t}$ , подставляем в правую часть $- 1000 (y - cos t)$ и собираем слагаемые — всё сводится к тождеству. Навык проверки точных решений пригодится и при тестировании численных методов: без эталона не отличить устойчивость от точности.

Почему явный Эйлер взрывается#

Явный Эйлер шагает так: $y_{n + 1} = y_{n} + h f (t_{n}, y_{n})$ . Для линейной компоненты $y^{'} = λ y$ каждый шаг умножает её на множитель перехода $1 + hλ$ — и это вся теория устойчивости в одной формуле. Пока $∣1 + hλ ∣ \leq 1$ , компонента затухает или стоит; как только множитель по модулю больше единицы — экспоненциальный рост, взрыв. Для $λ = - 1000$ граница $h \leq 2/∣ λ ∣ = 0, 002$ :

Шаг h	Множитель 1 − 1000h	Что происходит
$0, 01$	$- 9$	рост в 9 раз за шаг: $0 \to 10 \to - 80 \to 730 \to - 6560 \to 59050$
$0, 002$	$- 1$	граница устойчивости: осцилляции не затухают
$0, 001$	$0$	переходный режим гасится за один шаг — устойчиво
$0, 0005$	$0, 5$	плавное затухание, но счёт в 20 раз дольше

множители явного Эйлера для λ = −1000: три режима — взрыв, граница, затухание

Вдумайтесь в первую строку: шаг $h = 0, 01$ по точности великолепен — решение меняется на масштабе $1$ , и такая сетка разрешает его с запасом. Но множитель $- 9$ каждую итерацию растит ошибку в 9 раз; после десяти шагов накопление составило бы $9^{10} \approx 3, 5$ миллиарда. Числа в таблице — не мультяшка: прямая прокрутка явного Эйлера с $y_{0} = 0$ даёт именно $10, - 80, 730, - 6560, 59050$ .

Пограничные строки таблицы тоже поучительны. При $h = 0, 002$ множитель равен $- 1$ : компонента не растёт, но и не умирает — счёт ходит по кругу с чередованием знака, и любой изгиб правой части раскачивает решение дальше. При $h = 0, 001$ множитель ноль — «мёртвый» режим: переходная компонента уничтожается одним шагом. Это бывает полезно, когда быструю моду нужно просто выбросить из решения, а не отобразить.

Вот суть жёсткости одним предложением: явный метод обязан шагом обслуживать самый быстрый процесс задачи, даже если тот физически уже закончился. Точность позволила бы $h = 0, 01$ , устойчивость запрещает всё, что больше $0, 002$ , — и разрыв этот тем страшнее, чем больше отношение масштабов.

Неявный Эйлер: безусловно устойчив#

Лекарство меняет всего одну букву — правая часть берётся в новой точке: $y_{n + 1} = y_{n} + h f (t_{n + 1}, y_{n + 1})$ . Неизвестное теперь сидит в обеих частях, и для нашего уравнения шаг превращается в решение линейного уравнения:

y_{n + 1} (1 + 1000 h) = y_{n} + 1000 h cos (t_{n + 1}) \Rightarrow y_{n + 1} = \frac{y _{n} + 1000 h cos ( t _{n + 1} )}{1 + 1000 h}

неявный шаг для y' = −1000(y − cos t); для линейных уравнений формула явная, для общих — метод Ньютона на каждом шаге

Множитель перехода стал $\frac{1}{1 - hλ}$ : при $λ = - 1000$ и $h = 0, 01$ это $1/11 \approx 0, 09$ — затухание, причём при любом $h$ . Пока $Re λ \leq 0$ , знаменатель по модулю не меньше единицы, и метод гасит все устойчивые компоненты — это называется A-устойчивостью. Плата честная: каждый шаг дороже (уравнение надо решить), но шаг можно брать по точности, а не по выживанию.

Как распознать жёсткость на практике#

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

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

Числовой пример: пять шагов неявного Эйлера#

Берём ту же задачу, $y_{0} = 0$ , $h = 0, 01$ , $1000 h = 10$ , и крутим формулу $y_{n + 1} = (y_{n} + 10 cos t_{n + 1}) /11$ . Рядом — точное решение:

t	Неявный Эйлер	Точное решение	\|Ошибка\|
$0, 00$	$0, 00000$	$0, 00000$	—
$0, 01$	$0, 90905$	$0, 99991$	$9, 1 \cdot 1 0^{- 2}$
$0, 02$	$0, 99155$	$0, 99982$	$8, 3 \cdot 1 0^{- 3}$
$0, 03$	$0, 99882$	$0, 99958$	$7, 6 \cdot 1 0^{- 4}$
$0, 04$	$0, 99917$	$0, 99924$	$7, 3 \cdot 1 0^{- 5}$
$0, 05$	$0, 99879$	$0, 99880$	$1, 1 \cdot 1 0^{- 5}$

ошибка затухает примерно в 11 раз за шаг — ровно множитель 1/(1+1000h); расчёт проверен скриптом

Первый шаг грубоват — ошибка $0, 09$ : метод «замыливает» резкий старт, это свойство первого порядка. Но дальше ошибка умирает геометрически, и к $t = 0, 05$ расхождение — одна стотысячная. Досчитайте до $t = 0, 5$ : неявный Эйлер даёт $0, 878057$ против точного $0, 878061$ . А явный Эйлер с тем же шагом к этому моменту давно превратил бы задачу в генератор чисел порядка $1 0^{12}$ — разница в семь порядков на одном и том же $h$ .

Таблицу легко проверить на калькуляторе — это лучший способ прочувствовать механизм. Возьмите $y_{1}$ : $y_{1} = 10 cos (0, 01) /11 = 9, 9995/11 = 0, 90905$ . Дальше каждый шаг — одно деление: точность мала (первый порядок), но множитель $1/11$ методично гасит всё, что метод натворил на старте. Именно так жёсткие задачи и решаются промышленно: неявный метод на крупном шаге плюс умный контроль ошибки вместо тысяч микрошагов явного.

Что брать в бою#

Оцените спектр: выпишите $λ$ линейной части или прикиньте отношения констант скоростей.
Прогоните явный метод с автоматическим шагом: если шаг упёрся в нижний предел — диагноз жёсткость подтверждён.
Переключитесь на неявный: в SciPy это solve_ivp с method='Radau' (неявный Рунге–Кутта) или 'BDF' (формулы разной кратности); LSODA переключается сам, в MATLAB стандарт — ode15s.
Проверьте на модельной задаче с известным решением: наш пример с $cos t$ для этого идеален.
Контроль: уменьшите допуск и убедитесь, что решение меняется в пятом знаке, а не в третьем.

На экзамене и в работе

Три прокола. Первый — путать точность и устойчивость: малый шаг — не панацея, он лишь переносит взрыв чуть дальше по времени, если множитель по модулю больше единицы. Второй — забыть про цену неявности: на каждом шаге решается уравнение, и для нелинейной задачи без метода Ньютона не обойтись. Третий — считать неявный Эйлер панацеей точности: он первого порядка и первый шаг делает грубо; его сила — устойчивость, а тонкую точность дают Radau или BDF высоких порядков.

Проверь себя+15 XP

Явный Эйлер, задача $y^{'} = - 1000 (y - cos t)$ , шаг $h = 0, 01$ . Во сколько раз меняется переходная компонента за один шаг?

Проверь себя+15 XP

Какой максимальный шаг допускает явный Эйлер для этой задачи по устойчивости?

Частые вопросы

Жёсткое уравнение — это всегда нелинейное?

Нет, наш пример $y^{'} = - 1000 (y - cos t)$ линейный. Жёсткость — про разброс временных масштабов, а не про нелинейность: она возникает в линейных системах с сильно разными собственными значениями и в нелинейных задачах, где линеаризация имеет такой же разброс. Нелинейность лишь добавляет проблем неявным методам: на каждом шаге приходится решать уравнение методом Ньютона.

Чем неявный Эйлер отличается от явного по затратам?

Один явный шаг — одно вычисление правой части. Один неявный шаг — решение уравнения относительно $y_{n + 1}$ : для линейной задачи формула, для нелинейной — 2–3 итерации Ньютона с якобианом. Зато шаг неявный метод позволяет взять в сотни раз больше, и итог быстрее на порядки: выгодность определяет устойчивость, а не цена одного шага.

Что такое A-устойчивость простыми словами?

Метод A-устойчив, если его область устойчивости покрывает всю левую полуплоскость: любая компонента с $Re λ \leq 0$ не растёт при каком угодно шаге. У неявного Эйлера множитель $1/ (1 - hλ)$ по модулю не превышает единицы ровно на этой полуплоскости. Практический смысл: шаг выбирается только из точности, устойчивость перестаёт быть ограничением.

Какой метод выбрать в Python для жёсткой системы?

Первый заход — solve_ivp с method='Radau': неявный Рунге–Кутта пятого порядка с оценкой ошибки. BDF хорош, когда правая часть дешёвая, а решение гладкое, — формулы разной кратности экономят вычисления. LSODA сам детектирует жёсткость и переключается между явной и неявной ветками — удобен, если жёсткость появляется лишь на части интервала. Явные Рунге–Кутты (RK45) для жёстких задач не годятся по определению.

Куда дальше: механика явного шага и его погрешности — в уроке о методе Эйлера, общая теория того, почему решения растут или умирают, — в устойчивости решений. Численная кухня целиком — в численном решении ОДУ, откуда берутся ошибки — в источниках погрешностей. Формулы под рукой — шпаргалка по дифференциальным уравнениям; тренировка — тренажёр по дифурам и тренажёр численных методов. Общий маршрут по теме — в гайде как решать дифференциальные уравнения и задачах с решениями.