Обыкновенные дифференциальные уравнения

1687 год. Ньютон публикует Principia. $F=ma$ - это не формула, это ОДУ: $m\ddot{x} = F(x,\dot{x},t)$. Три века спустя та же запись - основа симуляции климата, нейронных сетей и орбиты МКС. Численное решение ОДУ - это мост между уравнениями физики и предсказаниями о реальном мире.

• **Neural ODE (torchdiffeq):** ResNet как дискретизация ОДУ, backprop через adjoint method - O(1) памяти вместо O(N шагов)
• **SpaceX Dragon:** траектория вычисляется решением ОДУ движения ракеты с 6 степенями свободы каждые 10 мс
• **Climate models:** coupled ОДУ для атмосферы + океана + льда (>100 000 уравнений), метод BDF для жёстких взаимодействий
• **Physics-Informed Neural Networks (PINNs):** ОДУ как loss term при обучении нейросети - модель обучается удовлетворять дифференциальному уравнению

Метод Эйлера: явный, неявный, симплектический

**1687 год. Ньютон публикует Principia. F=ma - это не формула, это ОДУ: $m\ddot{x} = F(x,\dot{x},t)$.** Три века спустя та же запись используется для симуляции климата, нейронных сетей и орбиты МКС. Задача: дана производная $y' = f(t, y)$ и начальное условие $y(t_0) = y_0$. Метод Эйлера - прямая дискретизация определения производной. Явный Эйлер: $y_{n+1} = y_n + h f(t_n, y_n)$. Шаг $h$ - размер временного интервала. Глобальная ошибка $O(h)$: уменьшение шага вдвое уменьшает ошибку вдвое.

**Область устойчивости явного Эйлера** - это часть комплексной плоскости, где $|1 + h\lambda| \leq 1$. Для вещественного отрицательного $\lambda$ это требует $h < 2/|\lambda|$. Если в системе есть собственное значение $\lambda = -10^6$ (жёсткая система), шаг должен быть меньше $2 \times 10^{-6}$ - даже если решение меняется медленно.

**Симплектический Эйлер** - модификация для гамильтоновых систем (физические симуляции): сначала обновляем импульс, затем координату через новый импульс. Глобальная ошибка $O(h)$ как у обычного Эйлера, но энергия сохраняется в среднем - без долгосрочного дрейфа. Именно этот метод используется в игровых движках.

Рунге-Кутта 4-го порядка: точность против жёсткости

**SpaceX Dragon: траектория вычисляется решением ОДУ с 6 степенями свободы каждые 10 мс. При отказе решателя - отклонение 1 метр через минуту. При правильной схеме - 1 сантиметр.** RK4 (метод Рунге-Кутты 4-го порядка) - стандарт для нежёстких задач. Идея: вместо одной оценки производной в начале шага (Эйлер), взять взвешенное среднее из 4 оценок в разных точках интервала $[t_n, t_n+h]$.

**RK4: четыре оценки производной.** $k_1 = f(t_n, y_n)$, $k_2 = f(t_n + h/2, y_n + hk_1/2)$, $k_3 = f(t_n + h/2, y_n + hk_2/2)$, $k_4 = f(t_n + h, y_n + hk_3)$. Итог: $y_{n+1} = y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)$. Глобальная ошибка $O(h^4)$ - при уменьшении шага вдвое точность растёт в 16 раз. За один шаг - 4 вычисления $f$.

Жёсткие системы: когда RK4 бесполезен

**Euler method - это не плохой метод. Для нежёстких систем с маленьким шагом он работает. Плохо то что 'маленький шаг' зависит от задачи. Для жёстких систем 'маленький' - это $10^{-12}$.** Жёсткость (stiffness) - отношение наибольшего к наименьшему собственному значению матрицы Якоби $\partial f / \partial y$. Химические реакции горения: быстрые H-OH реакции с $\lambda \sim 10^{10}$ и медленные изменения температуры с $\lambda \sim 1$. Индекс жёсткости $10^{10}$. RK45 потребует $10^{10}$ шагов на 1 секунду симуляции.

**BDF (Backward Differentiation Formulas)** - семейство неявных методов для жёстких систем. BDF-k использует $k$ предыдущих значений для аппроксимации производной. BDF-2: $\frac{3}{2}y_{n+1} - 2y_n + \frac{1}{2}y_{n-1} = h f(t_{n+1}, y_{n+1})$. Область A-устойчивости охватывает почти всю левую полуплоскость - шаг ограничен только точностью, не устойчивостью. scipy использует BDF через `method='BDF'` или Radau IIA (`method='Radau'`).

**scipy.integrate.solve_ivp с method='RK45' (по умолчанию) плохо работает для жёстких систем.** Если solve_ivp долго не завершается или требует миллионы шагов - система жёсткая. Переключитесь на `method='BDF'` или `method='Radau'`. Radau - метод 5-го порядка с превосходной A-устойчивостью, предпочтительнее BDF для точных расчётов.

Адаптивный шаг и Neural ODE

**Neural ODE (Chen et al. 2018) - это ResNet, где количество слоёв стремится к бесконечности. Backpropagation через Neural ODE - это решение сопряжённого ОДУ назад во времени.** Адаптивный контроль шага - ключевой инструмент: в гладких участках шаг большой, в резких - маленький. Метод Дормана-Принса (RK45): два решения разных порядков на одних вычислениях. Встроенная пара RK4/RK5: оценка локальной ошибки без дополнительных вычислений $f$.

**PI-регулятор для управления шагом.** Базовое правило: $h_{new} = h \cdot (\varepsilon_{tol} / \varepsilon_{local})^{1/(p+1)}$. Проблема - осцилляции шага. PI-контроллер стабилизирует: $h_{new} = h \cdot (\varepsilon_{tol}/e_n)^{k_I} \cdot (e_{n-1}/e_n)^{k_P}$, где $k_I \approx 0.3/(p+1)$, $k_P \approx 0.4/(p+1)$. RTOL - относительная (масштабируется с $|y|$), ATOL - абсолютная (защищает от деления на нуль при $y \approx 0$).

**Выбор RTOL/ATOL.** Типичные значения: rtol=1e-6, atol=1e-9. Слишком малые - миллионы шагов. Слишком большие - неточное решение. Для Neural ODE тренировки: rtol=1e-3, atol=1e-4 (скорость важнее точности). Для физических симуляций: rtol=1e-8, atol=1e-10.

ОДУ: ключевые идеи

• Явный Эйлер $y_{n+1} = y_n + h f(t_n, y_n)$: порядок $O(h)$, ограничение устойчивости $h < 2/|\lambda_{max}|$
• RK4: 4 оценки производной, порядок $O(h^4)$ - при уменьшении шага вдвое точность растёт в 16 раз; стандарт для нежёстких задач
• Жёсткие системы: индекс жёсткости = отношение max/min собственных значений; implicit методы BDF/Radau не имеют ограничения устойчивости на шаг
• Адаптивный шаг (Dormand-Prince RK45): встроенная оценка ошибки, dense output, PI-контроллер шага; основа Neural ODE через adjoint backprop

Связанные темы

ОДУ соединяют непрерывную математику с дискретными вычислениями и пронизывают физику, ML и численный анализ:

• Собственные значения и SVD — Жёсткость системы ОДУ = отношение собственных значений матрицы Якоби; implicit методы требуют LU-разложения на каждом шаге
• Цифровая обработка сигналов: дискретизация — Дискретизация по времени в ОДУ и теорема Найквиста - одна идея в разных контекстах
• Deep Learning: Neural ODE — Neural ODE = ResNet с бесконечной глубиной; backprop через adjoint method
• Сложность алгоритмов — Порядок точности O(h^p) - прямое применение Big-O анализа к методам численного решения

Вопросы для размышления

• Климатическая модель содержит 100 000 уравнений с временными масштабами от секунд (конвекция) до десятилетий (изменение климата). Какой метод и почему - RK45 или BDF? Как изменится выбор при переходе от прогноза погоды (дни) к симуляции климата (годы)?
• ResNet с 1000 слоями - это дискретизация ОДУ с шагом $h = 1/1000$. Neural ODE - это предел при $h \to 0$. Что теряется и что приобретается при переходе к Neural ODE? Как adjoint метод решает проблему памяти при backpropagation?
• Symplectic Euler сохраняет энергию в среднем, хотя имеет лишь первый порядок точности. RK4 более точен, но не сохраняет энергию. Для симуляции планетарных орбит на миллиарды лет - какой метод лучше и почему точность и сохранение энергии - разные требования?

Связанные уроки

• sci-04
• alg-01-big-o
• dsp-05
• dl-05
• nm-07