Три метода, три области
Для тестового уравнения dy/dx = λy три явных метода ОДУ имеют следующие области стабильности в комплексной плоскости hλ:
Метод Эйлера (первый порядок): область стабильности — диск |1 + hλ| ≤ 1, окружность радиуса 1 с центром в (-1, 0). Отрицательное вещественное hλ должно лежать в [-2, 0].
Метод Рунге-Кутты 2 (метод средней точки) (второй порядок): область стабильности |1 + hλ + (hλ)²/2| ≤ 1. Больше, чем диск Эйлера, но все ещё ограничена.
Метод Рунге-Кутты 4 (четвёртый порядок): область стабильности удовлетворяет |1 + hλ + (hλ)²/2 + (hλ)³/6 + (hλ)⁴/24| ≤ 1. Отрицательное вещественное hλ расширяется примерно до -2.785. Область существенно больше, чем у Эйлера.
Обратный метод Эйлера (неявный): область стабильности — вся комплексная плоскость кроме диска |1 - hλ|⁻¹ > 1, эквивалентно |1/(1-hλ)| ≤ 1. Для λ в левой полуплоскости (Re(λ) < 0) это безусловно стабильно — никаких ограничений на h от условия стабильности.
Функция усиления
Для любого метода Рунге-Кутты коэффициент усиления на шаг R(hλ) является полиномиальной аппроксимацией e^(hλ):
- Эйлер: R(z) = 1 + z (усечено на первой степени)
- RK2: R(z) = 1 + z + z²/2 (усечено на второй степени)
- RK4: R(z) = 1 + z + z²/2 + z³/6 + z⁴/24 (усечено на четвёртой степени)
Область стабильности — {z ∈ ℂ : |R(z)| ≤ 1}. Усиление истинного решения: |e^z| = e^(Re(z)). Для Re(z) < 0 (стабильное ОДУ) истинное решение затухает. Численный метод стабилен, если |R(z)| ≤ 1 — соответствуя поведению затухания.
Чисто мнимые собственные значения: колебательные системы
Многие физические системы имеют чисто мнимые собственные значения: λ = iω (колебания без затухания). Система масса-пружина, небесная механика, динамика маятника.
Для λ = iω: hλ = ihω лежит на мнимой оси.
Стабильность Эйлера на мнимой оси: |1 + ihω|² = 1 + (hω)² > 1 для любого h > 0. Эйлер нестабилен для любого размера шага на чисто мнимых собственных значениях. Вычисленное «колебание» растёт без границ.
Стабильность RK4 на мнимой оси: область стабильности расширяется примерно до |hω| ≤ 2.83 на мнимой оси. При достаточно малом h метод RK4 обрабатывает незатухающие колебания. Эйлер не может.
Эта геометрия объясняет, почему Эйлер не работает на консервативных системах (масса-пружина, орбиты, волновые уравнения) даже при малом h, в то время как RK4 справляется хорошо.
Геометрия жёстких задач
Жёсткая система ОДУ имеет собственные значения с очень различными величинами. Коэффициент жесткости: κ = max|Re(λᵢ)| / min|Re(λᵢ)| >> 1.
Почему жесткость дорога для явных решателей:
Стабильность требует h·max|λᵢ| ≤ C (где C зависит от метода). Наиболее отрицательное собственное значение устанавливает границу.
Точность для медленной динамики требует h·min|λᵢ| ≥ ε (правильное разрешение самого медленного режима).
Если κ велико, эти два требования заставляют использовать очень малый h: достаточно малый для стабильности быстрого режима, достаточно большой для выборки медленного режима. Количество шагов масштабируется с κ.
Геометрическая картина в спектре собственных значений: собственные значения якобиана ∂f/∂y образуют набор точек в комплексной плоскости. Область стабильности явного решателя должна содержать все точки h·λᵢ. Если собственные значения охватывают диапазон от -1 до -1000, область стабильности должна охватывать диапазон 1000 по вещественной оси — требуя очень малый h.
Неявные решатели: область стабильности обратного метода Эйлера охватывает всю левую полуплоскость. Все собственные значения с Re(λ) < 0 автоматически находятся внутри области стабильности независимо от h. Ограничение на h исходит только от точности, а не от стабильности.
Коэффициент жесткости и стоимость
Рассмотрим сеть химических реакций с быстрыми реакциями (временная шкала 10⁻⁶ с) и медленными реакциями (временная шкала 1 с).
Коэффициент жесткости: κ = 10⁶ / 1 = 10⁶.
С RK4 (предел стабильности h·|λ_max| ≤ 2.785): h_max = 2.785 / 10⁶ ≈ 2.8 × 10⁻⁶ с.
Для интегрирования 10 с времени реакции: шагов = 10 / (2.8 × 10⁻⁶) ≈ 3.6 × 10⁶.
С обратным методом Эйлера (безусловно стабильный): h может быть выбран для точности медленных реакций. h = 10⁻² с (100 выборок на масштабе 1 с). Шагов = 10 / 10⁻² = 1000.
Коэффициент стоимости: явный 3.6 миллиона шагов против неявного 1000 шагов — множитель 3600. Каждый неявный шаг требует решения линейной системы (стоимость на шаг выше), но общая стоимость намного ниже для очень жёстких задач.
Почему n-мерные трубки — не то, что вы думаете
В 2D «трубка» радиуса ε вокруг кривой C — это набор точек на расстоянии ε от C. Поперечное сечение — окружность радиуса ε. Объём трубки растёт пропорционально её длине.
В n измерениях геометрия трубки коренным образом меняется из-за явления из главы 9:
Парадокс n-мерного угла: в n-мерном пространстве почти весь объём n-мерного гиперкуба находится в углах — не в центральной области. По мере увеличения n доля объёма на расстоянии ε от центра стремится к нулю для любого фиксированного ε.
Применительно к трубкам решений ОДУ:
В 2D: если истинное решение проходит через центр трубки, большинство близлежащих точек находятся близко к кривой. Малые возмущения держат вас близко к истинному решению.
В высоких измерениях: большинство точек внутри ограничивающего бокса трубки фактически находятся далеко от кривой истинного решения. «Объём» трубки доминируется углами — регионами, которые далеко от центра в нескольких измерениях одновременно.
Следствие для моделирования: с 28 связанными ОДУ (проблема американского флота с перехватом), возмущение размера ε в каждом измерении может произвести полное смещение ε√28 ≈ 5.3ε от истинного решения. Трубку нужно понимать в терминах нормы L2 по всем измерениям, а не только максимального смещения в любом одном измерении.
Стабильность в высоких размерностях: система, в которой каждый компонент затухает независимо (каждое собственное значение имеет отрицательную вещественную часть), может всё ещё показывать большие объединённые смещения, потому что ошибки компонентов суммируются в норме L2. 28-мерная трубка — это не просто 28 независимых 1-мерных трубок — геометрия их связывает.
От геометрии к проектированию
Геометрические идеи глав 18-20 объединяются в набор принципов проектирования для численного моделирования:
Выбор размера шага: h должен помещать h·λ внутри области стабильности для каждого собственного значения. Для жёстких систем неявные методы снимают ограничение стабильности, оставляя только требования точности.
Накопление ошибки в высоких размерностях: глобальная ошибка — это вектор в n-мерном пространстве. Его норма растёт как √n раз ошибка на компонент. Моделирование в высоких размерностях требует более жёстких требований к точности на шаг.
Обратная связь как стабилизатор: если моделирование включает обратную связь (вычисленный выход влияет на последующие входы, как в системе управления), сходящаяся обратная связь подавляет ошибки. Моделирование может допускать неточные входы для величин внутри цикла обратной связи.
Нестабильность как сигнал: для задач с расходящимися направлениями поля нестабильность может быть использована: направление расходимости несёт информацию об ошибке начального условия, позволяя корректирующей настройке.