Численные методы
Адаптивные методы и управление шагом
Моделирование орбит спутников: в перигелии (близко к планете) сила притяжения меняется быстро, нужен малый шаг. В апогелии, медленно, большой шаг безопасен. Фиксированный шаг, либо 10^6 лишних вычислений, либо недостаточная точность у планеты. Адаптивные методы решают этот баланс автоматически: scipy.integrate.solve_ivp сам выбирает шаги, достаточно задать rtol и atol.
- **Астродинамика**: scipy.integrate используют для расчётов траекторий, адаптивный RK45 справляется там, где равномерная сетка провалится
- **МКЭ в инженерии**: ANSYS, ABAQUS используют hp-адаптацию для концентраторов напряжений у отверстий и трещин, без адаптации пришлось бы строить мелкую сетку везде
- **Климатическое моделирование**: adaptive timestepping в CESM/MOM6, шаг по времени адаптируется к числу Куранта, предотвращая нестабильность
Предварительные знания
Адаптивный шаг для ODE: embedded methods
Нейронные ODE (Neural ODE, NeurIPS 2018) используют адаптивный решатель жёстких уравнений: DOPRI-5 c шагом от 10⁻⁶ до 1 в одном прямом проходе. Решение ODE может иметь области с разным «темпом изменения»: в одних местах решение почти постоянно (большой шаг безопасен), в других, быстро осциллирует или имеет перегиб (нужен малый шаг). Фиксированный шаг, всегда компромисс: либо слишком медленно, либо слишком неточно. Адаптивный шаг решает оба недостатка.
Идея: вычислить два приближения разного порядка за почти одно число вычислений функции f. Разность = оценка локальной ошибки. **RK45 (Dormand-Prince, dopri5)**: - Решение порядка 4: y̅ₙ₊₁ (использует 6 вычислений f) - Решение порядка 5: yₙ₊₁ (использует те же 6 + 1 = 7 вычислений f) - Локальная ошибка: err = ||yₙ₊₁ - y̅ₙ₊₁|| Эта пара называется «embedded», метод 5-го порядка «встроен» в метод 4-го. scipy.integrate.solve_ivp использует именно dopri5 по умолчанию.
**Правило обновления шага**: при ошибке err и допуске tol, оптимальный следующий шаг h_new = h · (tol/err)^(1/(p+1)), где p, порядок метода. На практике добавляют safety factor 0.9 и ограничивают рост: h_new ∈ [h/5, 5h]. Шаг отклоняется (rejected step) если err > tol, принимается и повторяется с уменьшенным h.
Зачем в embedded методах RK45 нужны два решения разного порядка?
Адаптивное измельчение сетки для PDE
Для PDE адаптивность сложнее: нужно не просто изменить шаг по времени, но и перестроить пространственную сетку. Элементы с большой ошибкой измельчаются, с малой, укрупняются. Это особенно важно для задач с локальными особенностями: ударные волны, фронты реакций, трещины в механике.
**h-refinement** (h-адаптация): добавление новых узлов (уменьшение h) в элементах с большой ошибкой. Стандартный МКЭ, алгоритм Rivara для триангуляций. **p-refinement** (p-адаптация): повышение порядка полинома на элементе. Особенно эффективно для гладких решений, экспоненциальная сходимость. **hp-refinement**: сочетание обоих. Там, где решение гладкое, повышать p (получить спектральную сходимость), где особенность, дробить h. Теоретически: экспоненциальная сходимость в числе степеней свободы даже для задач с сингулярностями. **r-refinement** (relocation): перемещение узлов без изменения их числа. Equidistribution principle: переместить узлы так, чтобы ошибка была одинакова на каждом элементе.
В чём преимущество hp-refinement над чистым h-refinement?
Оценка ошибки: Richardson extrapolation и a posteriori estimators
Для управления адаптацией нужна оценка ошибки. Две стратегии: **a priori** (до вычисления, через теоретические оценки) и **a posteriori** (после вычисления, используя само найденное решение). A posteriori оценки намного мощнее для адаптации.
**Residual-based**: оценка через невязку уравнения: err_K ≈ h²||R(u_h)||_{L2(K)} + h^{3/2}||jump(∇u_h)||_{L2(∂K)}. Включает объёмную невязку и скачок потока на гранях. **Recovery-based (ZZ)**: сравнить производную u_h с «суперсходящейся» производной σ* = Gh(∇u_h), полученной L2-проекцией. Практически более точно и проще реализовать. **Goal-oriented**: оценка ошибки в конкретном функционале J(u) (например, давление в точке, подъёмная сила). Использует сопряжённую задачу, получить точный ответ на интересующий вопрос.
Метод Ричардсона (Richardson extrapolation) требует решения задачи дважды, с шагами h и h/2. Что даёт нам эта пара решений?
Управление сходимостью: tolerances и баланс работы
Практика: как правильно задавать допуски и интерпретировать их? rtol и atol в ODE-решателях, tol в нелинейных решателях, mesh tolerance в МКЭ, всё это разные понятия, и путаница здесь стоит дорого.
Оптимальная адаптивная стратегия: **equidistribution**, ошибка одинакова на каждом элементе/шаге. Для ODE: оптимальный шаг hᵢ ∝ (tol / C·|f^{(p+1)}|)^{1/(p+1)}, больше там, где решение гладкое, меньше, там, где быстро меняется. Алгоритм r-refinement перемещает узлы до достижения equidistribution. Критерий: E_K = const для всех K, где E_K, оценка ошибки на K. **Баланс точности и стоимости**: цель, минимизировать число вычислений f при заданной глобальной ошибке. Теоретический минимум достигается именно при equidistribution.
Почему нельзя использовать только rtol=1e-8 (без atol) для ODE, решение которого проходит через ноль?
Ключевые идеи
- **Embedded RK (RK45/dopri5)**: два решения разного порядка за почти одно число вычислений f → оценка ошибки → автоматическое управление шагом
- **h/p/hp-refinement**: разные стратегии улучшения PDE-сетки, hp даёт экспоненциальную сходимость даже при сингулярностях
- **Richardson extrapolation**: два решения с h и h/2 → ошибка и решение порядка p+1, «бесплатный» порядок точности
- **rtol + atol**: правило допуска |err| ≤ atol + rtol·|y|, никогда не задавать только rtol для решений, проходящих через ноль
Связанные темы
Адаптивность пронизывает все численные методы:
- Методы решения ODE — RK45 и LSODA, адаптивные ODE-решатели, реализующие идеи embedded methods на практике
- МКЭ и PDE — A posteriori оценки, основа адаптивного МКЭ, реализованного в deal.II, FEniCS, DUNE
- Жёсткие системы — Адаптивный шаг особенно важен для жёстких уравнений: LSODA автоматически переключается между нежёстким и жёстким решателем
Вопросы для размышления
- Если данный адаптивный ODE-решатель вдруг делает тысячи крошечных шагов в одной точке, что это означает физически, и как это диагностировать?
- При каких условиях h-refinement лучше p-refinement для PDE, и наоборот?
- Как бы выбрали atol для системы ОДУ, где разные компоненты имеют разные физические масштабы (например, концентрации 10^{-9} М и температуры 300 К)?