Девлог №93 – Хвиля 72: фільтр Калмана, зворотне поширення, 3D-маятник, евтрофікація, акустична лінза та планетарні кільця

Нові симуляції

Фільтр Калмана: оптимальне оцінювання стану

Фільтр Калмана (1960) — один із найуспішніших алгоритмів прикладної математики. Він оптимально поєднує передбачення динамічної моделі із зашумленими показами датчиків, даючи незміщену оцінку стану з мінімальною дисперсією — і, що важливо, точно відслідковує, наскільки він невпевнений, через матрицю коваріації P.

4D-стан та операції з повними матрицями

Вектор стану — [x, y, v_x, v_y]: позиція та швидкість у 2D. Матриця коваріації P 4×4 відслідковує невизначеності та кореляції між усіма компонентами стану. Кожен крок передбачення обчислює:

x̂⁻ = F·x̂          (передбачення стану через матрицю переходу F)
P⁻  = F·P·Fᵀ + Q   (передбачення коваріації, Q = шум процесу)

Кожен крок оновлення (при отриманні спостереження 2D-позиції через H = [I₂|0₂]):

K  = P⁻·Hᵀ·(H·P⁻·Hᵀ + R)⁻¹  (коефіцієнт підсилення Калмана)
x̂  = x̂⁻ + K·(z − H·x̂⁻)       (оновлення стану)
P  = (I − K·H)·P⁻·(I − K·H)ᵀ + K·R·Kᵀ  (форма Джозефа)

Форма Джозефа для числової стабільності

Стандартне оновлення P = (I−KH)P⁻ страждає від асиметрії через похибки з плаваючою комою — P може втратити додатну визначеність після багатьох ітерацій. Форма Джозефа P = (I−KH)P⁻(I−KH)ᵀ + KRKᵀ алгебраїчно еквівалентна, але числово стабільна: вона гарантує, що P лишається симетричною додатно визначеною незалежно від округлення з плаваючою комою.

Візуалізація еліпса коваріації

2D невизначеність позиції візуалізується як еліпс, обчислений із власного розкладу підматриці позиції 2×2 матриці P. Напівосі — √(λ₁) і √(λ₂) (стандартні відхилення), а орієнтація слідує власним векторам. Еліпс візуально роздувається під час кроку передбачення й стискається, коли надходить спостереження — прекрасна демонстрація інформаційного циклу фільтра.

• 4 траєкторії: пряма лінія, кругова орбіта, випадкове блукання, ухильний маневр
• Порівняння RMSE: оцінка Калмана проти сирих вимірів проти лише передбачення моделі
• Повзунки шуму процесу Q та шуму вимірювання R
• Перемикання частоти спостережень, щоб побачити зростання коваріації під час прогалин

Зворотне поширення: потік градієнтів, зроблений видимим

Зворотне поширення — це просто ефективно застосоване ланцюгове правило до графа обчислень — але щоб зрозуміти його інтуїтивно, треба побачити, як градієнтні імпульси течуть назад крізь кожен шар. Ця симуляція робить цей потік явним.

Ручна реалізація ланцюгового правила

Архітектура мережі налаштовується: 2 входи → H прихованих юнітів → H прихованих юнітів → 2 (класифікація) чи 1 (регресія) виходи. Кожен прямий прохід обчислює активації пошарово. Кожен зворотний прохід обчислює похибки δ: вихідний шар отримує δ_L = ∂L/∂z_L, а далі поширюється: δ_l = (Wᵀ_{l+1} · δ_{l+1}) ⊙ σ'(z_l). Градієнти ваг: ∂L/∂W_l = δ_l · aᵀ_{l-1}.

Візуальна анімація градієнтів

Під час зворотного поширення кольорові імпульси рухаються справа наліво крізь діаграму мережі. Ширина ребра й інтенсивність кольору відображають |∂L/∂w| для кожної ваги. Великі градієнти світяться яскраво; майже нульові градієнти (зникаючий градієнт) показуються тьмяними лініями. З ReLU мертві нейрони (завжди на виході дають 0) позначаються явно — їхні вихідні градієнти назавжди нульові.

• Активації: ReLU, tanh, sigmoid; вихід: softmax+CE або лінійний+MSE
• Теплова карта межі рішень оновлюється кожні N кроків навчання
• Графік смуги втрат: втрати навчання з часом
• Покроковий режим одного зразка для вивчення алгоритму

3D-маятник: обмеження через множник Лагранжа

Сферичний маятник — маса на жорсткому стрижні, вільна коливатись у всіх напрямках — простий для фізичного опису, але незручний для числової симуляції. Обмеження полягає в тому, що відстань від точки підвісу до маси завжди точно дорівнює L. Замість переходу до сферичних координат (що вносить сингулярності gimbal-lock), симуляція працює в декартових координатах і забезпечує обмеження через множник Лагранжа.

Декартовий RK4 з проєкцією обмеження

Рівняння руху в декартових (x,y,z) координатах: m·r̈ = F_гравітація + λ·r̂, де λ — величина сили обмеження (натяг). На кожному кроці RK4, після вільного просування r та v, позиції проєктуються назад на сферу, а швидкості — на дотичну площину: λ = (−|v|² − r·g_vec) / L², потім застосовується поправка. Це гарантує точне виконання обмеження без координатних сингулярностей.

Три фізичні режими

• Плоска розетка: початкові умови в одній площині — при погляді згори утворює візерунок розетки, схожий на фігуру Лісажу, через нелінійний період
• Конічний маятник: маса обертається під фіксованим кутом θ з періодом T = 2π√(L cosθ / g) — коротшим за період простого маятника через зменшену ефективну гравітацію
• Маятник Фуко: додає силу Коріоліса F_c = −2mΩ × v з Ω = Ω_Земля·sin(широта)·ẑ, змушуючи площину коливань прецесувати зі швидкістю Ω·sin(широта)

• 3D-проєкція з обертанням перетягуванням миші
• Повзунок широти для режиму Фуко (0° = екватор, 90° = полюс)
• Довжина сліду та колірний градієнт за часом

Евтрофікація: точка перекидання озера

Евтрофікація — збагачення поживними речовинами, що веде до цвітіння водоростей та виснаження кисню — один із найважливіших прикладів катастрофічної точки перекидання в екологічній науці. Щойно озеро перетинає поріг, воно різко переходить у каламутний, бідний киснем стан, що може зберігатись десятиліттями навіть після зниження надходження фосфору.

П'ятивидова система ОДУ

Симуляція інтегрує п'ять пов'язаних ОДУ, що представляють P (розчинений фосфор), A (біомаса водоростей), O (розчинений кисень), F (біомаса риби) та D (детрит). Ключова нелінійність — член рециркуляції фосфору: за аноксійних умов (O нижче порогу) донні відклади вивільняють накопичений фосфор — позитивний зворотний зв'язок, що штовхає систему до альтернативного стійкого стану. Інтегратор RK4 надійно обробляє жорсткі переходи.

Гістерезис та фазовий графік

Фазовий графік простежує (навантаження фосфором, біомаса водоростей) з часом. Коли навантаження повільно зростає, система лишається у стані чистої води, доки різко не перескакує в стан цвітіння в точці верхньої сідло-вузлової біфуркації. Коли потім навантаження зменшується, система лишається у стані цвітіння аж до нижньої сідло-вузлової точки — це і є гістерезис. Площа, охоплена петлею гістерезису на фазовому графіку, представляє "борг відновлення" — скільки додаткових зусиль потрібно, щоб відновити озеро, порівняно із запобіганням цвітінню від самого початку.

• Візуалізація водойми: градієнт кольору води від прозоро-блакитного до каламутно-зеленого, поверхнева плівка цвітіння, іконки риб, темрява донних відкладень
• Повзунок навантаження фосфором з режимами повільного зростання та повільного зниження
• Показники рівня кисню та популяції риб у реальному часі

Акустична лінза: FDTD-симуляція хвиль

Акустичне фокусування — концентрація звукових хвиль у точці — лежить в основі медичного ультразвуку, сонарів та архітектурної акустики. Ця симуляція розв'язує 2D хвильове рівняння напряму методом Finite Difference Time Domain (FDTD), показуючи вигинання фронтів хвиль у реальному часі.

Leapfrog FDTD та межі Мура

2D хвильове рівняння тиску ∂²p/∂t² = c²(∂²p/∂x² + ∂²p/∂y²) дискретизується на сітці 240×144 за схемою leapfrog: тиск у момент n+1 залежить від моментів n і n-1 та просторового лапласіана в n. Умова стабільності CFL вимагає c·dt/dx ≤ 1/√2. Поглинальна межова умова Мура першого порядку усуває відбиття від країв сітки, узгоджуючи швидкість вихідної хвилі: p_край(t+dt) = p_{край+1}(t) + (c·dt−dx)/(c·dt+dx)·(p_{край+1}(t+dt) − p_край(t)).

Чотири типи лінз

• Опукла лінза: вища швидкість звуку в області лінзи (матеріал нижчої густини) сходить фронт хвилі
• Увігнута лінза: той самий матеріал, але увігнута геометрія розсіює промінь
• Лінза Френеля: дизайн зонної пластини — концентричні кільця з чергуванням високої/низької швидкості, що складаються синфазно у фокусній точці
• GRIN (параболічна): швидкість змінюється як c(r) = c₀ / (1 + (r/r₀)²) — лінза з градієнтним показником заломлення без дискретних меж

• Дивергентна теплова карта: поле тиску синій/червоний з адитивним світінням
• Накладання променів за законом Снелла: наближення геометричної оптики поверх хвильової симуляції
• EMA тиску в передбаченій фокусній точці підтверджує фокусну відстань за параксіальною формулою

Планетарні кільця: кеплерівські частинки та резонансні проміжки

Кільця Сатурна — одна з найфотогенічніших структур Сонячної системи, і їхня структура контролюється орбітальною механікою. Проміжок Кассіні, помітний розрив між кільцями A та B, спорожнений орбітальним резонансом 2:1 із місяцем Мімасом.

Кеплерівські частинки та резонансні поштовхи

До 10 000 тестових частинок розміщуються на кеплерівських орбітах з великою піввіссю a, розподіленою як n ∝ a^(−3/2) (плоска поверхнева густина). Кожна частинка отримує гравітаційний поштовх від місяця на кшталт Мімаса при кожному з'єднанні. На радіусі резонансу 2:1 (де частинка робить рівно 2 оберти на один оберт місяця) ці поштовхи завжди трапляються в тій самій орбітальній фазі — вони накопичуються когерентно й підвищують ексцентриситет, доки орбіта не почне перетинатись з іншими, після чого частинка швидко розсіюється. Після тисяч орбіт радіус резонансу показує розрив густини, що збігається з проміжком Кассіні.

Межа Роша та місяці-пастухи

Межа Роша R_Роша = a_планета · (2ρ_планета/ρ_супутник)^(1/3) позначає межу, за якою припливні сили перевищують власну гравітацію — місяці всередині неї подрібнюються на матеріал кілець. Симуляція також підтримує пари місяців-пастухів: два малих місяці, що обрамляють вузьке кільце, обмінюються моментом імпульсу з частинками кільця, утримуючи кільце чітко звуженим — механізм, що стоїть за структурою F-кільця Сатурна.

• Повзунок N: від 100 до 10 000 частинок
• Маса місяця та орбітальний радіус налаштовуються — спостерігайте, як формується й рухається проміжок
• Забарвлення за ексцентриситетом для візуалізації резонансного збудження
• Перемикач пари місяців-пастухів для звуження вузького кільця

Що далі

Хвиля 73 — спецвипуск WebGL/GLSL: усі п'ять симуляцій використовують власні шейдери, raymarching або GPGPU-техніки. Далі на черзі: raymarching 3D-фрактала Mandelbulb з оцінювачем відстані, SDF губки Менгера з розмірністю Хаусдорфа, вкладення параболоїда Флемма для кротовини Морріса-Торна, сферичний добуток суперформи Гіліса та GPGPU-симуляція на 250 000 частинок з ping-pong render targets.