Три методи, три області
Для тестового рівняння 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 (скорочено на ступені 1)
- RK2: R(z) = 1 + z + z²/2 (скорочено на ступені 2)
- RK4: R(z) = 1 + z + z²/2 + z³/6 + z⁴/24 (скорочено на ступені 4)
Область стійкості — це {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 зв'язаними ОДЗ (проблема перехоплення Hamming Navy), збурення розміру ε в кожному вимірі може виробити загальне зміщення ε√28 ≈ 5,3ε від справжнього розв'язку. Трубку повинно розумітися в термінах L2 норми вся розмірів, не тільки максимального зміщення в будь-якому одному вимірі.
Стійкість у високих вимірах: система, де кожен компонент затухає незалежно (кожне власне значення має від'ємну дійсну частину), все ж може показати великі комбіновані зміщення, тому що помилки компонентів додаються в L2 нормі. 28-вимірна трубка — це не просто 28 незалежних 1-вимірних трубок — геометрія їх з'єднує.
Від геометрії до дизайну
Геометричні поглади розділів 18-20 збираються в набір принципів дизайну для числової симуляції:
Вибір розміру кроку: h повинен розмістити h·λ всередину області стійкості для кожного власного значення. Для жорстких систем, неявні методи видаляють обмеження стійкості, залишаючи тільки вимоги точності.
Накопичення помилок у високих вимірах: глобальна помилка — вектор у n-вимірному просторі. Її норма зростає як √n разів помилка на компонент. Високовимірні симуляції потребують більш жорстких вимог точності на крок.
Зворотний зв'язок як стабілізатор: якщо симуляція включає зворотний зв'язок (обчислений вихід впливає на подальші входи, як у системі керування), конвергентний зворотний зв'язок гасить помилки. Симуляція може толерувати неточні входи для величин у межах циклу зворотного зв'язку.
Нестійкість як сигнал: для задач з розбіжними полями напрямку, нестійкість може бути використана: напрямок розбіжності несе інформацію про помилку початкової умови, дозволяючи коригуючу настройку.