Девлог #84 – Хвиля 63: Течія Стокса, упаковка сфер, крива Гільберта, спіраль Улама, максимальний потік і крижаний ореол

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

🌊 Течія Стокса — фізика повзучого руху

Рівняння Стокса (Re → 0)

При дуже низькому числі Рейнольдса Re = ρUL/μ ≪ 1 інерційні сили нехтовно малі порівняно з в'язкими силами, і рівняння Нав'є-Стокса зводяться до лінійних рівнянь Стокса:

∇p = μ ∇²u       (імпульс: градієнт тиску врівноважує в'язку дифузію)
∇ · u = 0        (неперервність: нестисливий потік)

У 2D формі функції течії (u = ∂ψ/∂y, v = −∂ψ/∂x):
  ∇⁴ψ = 0         (бігармонічне рівняння)

Оскільки рівняння лінійні й не містять часових похідних (час з'являється лише через граничні умови), течія Стокса є квазістатичною: поле швидкості миттєво адаптується до граничних умов. Це має глибокий наслідок — оборотність.

Функція течії для потоку навколо сфери

Для рівномірного потоку U навколо сфери радіусом R точна функція течії Стокса в сферичних координатах (r, θ):

ψ(r, θ) = U · sin²θ · [r²/2  -  3R·r/4  +  R³/(4r)]

Компоненти швидкості:
  u_r = (1/r² sinθ) · ∂ψ/∂θ = U cosθ · [1 - 3R/(2r) + R³/(2r³)]
  u_θ = -(1/r sinθ) · ∂ψ/∂r = -U sinθ · [1 - 3R/(4r) - R³/(4r³)]

Сила опору (опір Стокса):
  F_drag = 6π μ R U

Лінії течії симетричні спереду й позаду сфери — властивість, унікальна для течії Стокса. Натомість при високому Re слід порушує симетрію, і опір різко зростає.

Демонстрація оборотності

Найяскравіша риса симуляції — демонстрація оборотності. Кольорові плями барвника впорскуються й переносяться вперед полем швидкості Стокса. Коли потік обертається (U → −U), частинки барвника точно повторюють свої шляхи назад і повторно сходяться до початкових позицій — система не має ефективного «змішування», тому що лінійність течії Стокса означає відсутність хаотичного розтягування й складання. Це знаменитий експеримент Тейлора з кукурудзяним сиропом 1966 року. Та сама демонстрація з течією Нав'є-Стокса при Re = 100 показала б плями, незворотно розмазані турбулентним змішуванням.

Теорема про гребінець

Одна горизонтальна лінія барвника поділяється на N сегментів, кожен з яких проштовхується крізь окремий паралельний циліндр. Після проходження крізь N циліндрів N окремих ниток барвника возз'єднуються в одну горизонтальну лінію — ідентичну оригіналу. Це наслідок лінійності й оборотності течії Стокса: шлях кожної частинки рідини крізь перешкоду в течії Стокса однозначно визначається її позицією входу й оборотний, тож позиції виходу гарантовано возз'єднуються в порядку.

🔵 Упаковка сфер — гіпотеза Кеплера і коефіцієнти упаковки

Коефіцієнти упаковки в 2D

У 2D відомі дві регулярні (ґраткові) упаковки. Квадратна упаковка розміщує центри кіл на сітці одиничного квадрата; гексагональна (трикутна) упаковка зсуває чергові ряди на половину діаметра:

Квадратна упаковка:
  центри в (i, j) для i,j ∈ ℤ; радіус r = 0.5
  площа кола / площа одиничної клітини = π r² / 1² = π/4 ≈ 0.7854

Гексагональна упаковка:
  центри в (i + j/2, j·√3/2); радіус r = 0.5
  площа одиничної клітини = √3/2 (ромб)
  коефіцієнт упаковки = π r² / (√3/2) = π/(2√3) ≈ 0.9069

Випадкове послідовне додавання (RSA):
  Розміщувати диски по одному у випадкових позиціях, якщо немає перекриття.
  Коефіцієнт упаковки заклинювання ≈ 0.547 (стала паркування Реньї)

Гексагональна упаковка доведено оптимальна в 2D — жодне розташування однакових кіл не може досягти вищої густини. Це довів Туе в 1910 році. 2D-симуляція показує живий коефіцієнт упаковки та характерний гексагональний порожній візерунок між колами.

3D-кристалічні ґратки та гіпотеза Кеплера

3D-задача набагато складніша. Чотири класичні ґраткові упаковки, реалізовані в симуляції:

Проста кубічна (SC):
  a = 2r;  коефіцієнт упаковки = π/6 ≈ 0.5236;  координаційне число Z = 6

Об'ємноцентрована кубічна (BCC):
  a = 4r/√3;  коефіцієнт упаковки = π√3/8 ≈ 0.6802;  Z = 8

Гранецентрована кубічна (FCC):
  a = r√8;  коефіцієнт упаковки = π/(3√2) ≈ 0.7405;  Z = 12

Гексагональна щільна упаковка (HCP):
  Той самий коефіцієнт упаковки, що й FCC = π/(3√2) ≈ 0.7405;  Z = 12
  Відрізняється від FCC лише послідовністю укладання (ABAB проти ABCABC)

Кеплер припустив у 1611 році, що π/(3√2) ≈ 74.05% — максимально можливий коефіцієнт упаковки для однакових сфер. Це формально довів Томас Гейлс у 1998 році, використовуючи комп'ютерно-асистований доказ на 250 сторінках плюс 3 ГБ коду — один з найдовших доказів у математиці. І FCC, і HCP досягають цього максимуму.

3D-проєкція та візуалізація контактів

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

🌀 Крива Гільберта — криві, що заповнюють простір, і локальність

Побудова рекурсивною підстановкою

Крива Гільберта генерується граматикою підстановки, застосованою до початкового патерну. У порядку 1 це ∪-подібна форма, що з'єднує 4 точки. У порядку n кожен з 4 підквадратів з порядку n−1 замінюється поверненою/відображеною копією кривої порядку n−1, з'єднаних мостовими сегментами. Крива порядку n відвідує 4ⁿ точок і має загальну довжину (4ⁿ − 1) × (розмір одиничної клітини), яка розходиться при n → ∞, тоді як крива залишається неперервною та обмеженою одиничним квадратом.

// відображення d-на-(x,y) (маніпуляція бітами, за Скіллінгом 2004):
function d2xy(n, d):
    rx, ry, s = 0, 0, 1
    x, y = 0, 0
    поки s < n:
        rx = (d >> 1) & 1
        ry = (d ^ rx) & 1
        rotate(s, x, y, rx, ry)
        x += s * rx
        y += s * ry
        d >>= 2
        s <<= 1
    повернути x, y

// Обернене (x,y)-на-d:
function xy2d(n, x, y):
    rx, ry, s = 0, 0, n >> 1
    d = 0
    поки s > 0:
        rx = (x & s) > 0 ? 1 : 0
        ry = (y & s) > 0 ? 1 : 0
        d += s * s * ((3 * rx) ^ ry)
        rotate(s, x, y, rx, ry)
        s >>= 1
    повернути d

Збереження локальності

Ключова властивість, яка робить криву Гільберта практично корисною (в управлінні кешем CPU, R-деревах баз даних, просторовому індексуванні ГІС), полягає в тому, що точки, близькі одна до одної на 1D-кривій, також близькі одна до одної в 2D-просторі, і навпаки. Це збереження локальності можна кількісно оцінити: для двох 2D-точок з евклідовою відстанню d їхні індекси Гільберта відрізняються щонайбільше на O(d). Натомість рядкове («растрове») впорядкування має найгірші пари на протилежних кінцях простору індексів, розділені лише одним рядком у 2D, а Z-порядок (Мортон) є проміжним.

Криві Пеано та Мортона

Крива Пеано (1890, перша коли-небудь відкрита крива, що заповнює простір) поділяє кожен квадрат на 9 підквадратів у зміяному візерунку, відвідуючи 9ⁿ точок у порядку n. Крива Z-порядку Мортона перемежовує біти координат x та y:

// Кодування Мортона (Z-порядок):
morton(x, y) = spread_bits(x) | (spread_bits(y) << 1)

spread_bits(v):  // вставити біт 0 між кожним бітом v
    v = (v | (v << 8)) & 0x00FF00FF
    v = (v | (v << 4)) & 0x0F0F0F0F
    v = (v | (v << 2)) & 0x33333333
    v = (v | (v << 1)) & 0x55555555
    повернути v

Підказка при наведенні показує як пряме відображення (відстань кривої d → позиція сітки (x,y)), так і обернене, дозволяючи одним поглядом побачити, як зберігається чи порушується локальність для будь-якої обраної точки.

🔢 Спіраль Улама — прості візерунки й поліноміальні хребти

Побудова спіралі Улама

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

n=1 → (0,0)
n=2 → (1,0)  n=3 → (1,1)  n=4 → (0,1)  n=5 → (-1,1)
n=6 → (-1,0) n=7 → (-1,-1) n=8 → (0,-1) n=9 → (1,-1) ...

Загальна намотка (рух вправо 1, вгору 1, вліво 2, вниз 2, вправо 3, ...):
  довжини сегментів: 1,1,2,2,3,3,4,4,...
  напрямки: Сх,Пн,Зх,Пд,Сх,Пн,Зх,Пд,...

Симуляція обчислює позиції спіралі для всіх N до ~90 000, використовуючи попередньо обчислену таблицю намотки — достатньо швидко, щоб відрендерити всю сітку з частотою 60 fps з інтерактивним масштабуванням і панорамуванням.

Чому діагоналі? Квадратичні поліноми

Діагоналі відповідають квадратичним поліномам від n. Уздовж однієї головної діагоналі: значення f(k) = 4k² + 2k + 1 для k = 0, 1, 2, … Уздовж іншої: f(k) = 4k² − 2k + 1. Будь-який квадратичний поліном, що дає багато простих чисел, засвітиться яскравою діагоналлю на спіралі Улама. Найзнаменитіший приклад — поліном Ейлера f(n) = n² + n + 41, який простий для n = 0, 1, 2, …, 39 (40 послідовних значень) — рекорд, що протримався сторіччями. На спіралі це виглядає як безпомилково яскравий діагональний хребет, перемикається в симуляції.

Спіраль Сакса

Спіраль Сакса (Роберт Сакс, 1994) відображає ціле число n на полярні координати, розміщуючи його під кутом θ = 2π√n і радіусом r = √n — тож повні квадрати лежать на додатній осі x, а спіраль робить рівно один оберт за кожен повний квадрат:

θ(n) = 2π √n
r(n) = √n
x(n) = r cosθ,  y(n) = r sinθ

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

📦 Максимальний потік — Едмондс-Карп і теорема про мінімальний розріз

Потокові мережі та задача максимального потоку

Потокова мережа — це орієнтований граф G = (V, E), де кожне ребро (u, v) має невід'ємну пропускну здатність c(u, v). Потік f(u, v) на кожному ребрі має задовольняти: (1) обмеження пропускної здатності f(u, v) ≤ c(u, v), та (2) збереження потоку: потік до будь-якого внутрішнього вузла дорівнює потоку з нього. Задача максимального потоку запитує максимальний загальний потік від джерела s до стоку t.

Теорема максимальний потік = мінімальний розріз (Форд-Фалкерсон, 1956):
  Максимальний s-t потік дорівнює мінімальній ємності s-t розрізу.

s-t розріз: розбиття (S, T) V з s ∈ S, t ∈ T
Ємність розрізу = Σ c(u,v) для ребер (u,v) з u ∈ S, v ∈ T

Алгоритм Едмондса-Карпа

Алгоритм Едмондса-Карпа (1972) — це Форд-Фалкерсон з конкретним вибором доповнювального шляху: завжди використовувати найкоротший шлях (найменше ребер) від s до t у залишковому графі, знайдений за допомогою BFS. Цей вибір гарантує поліноміальну складність за часом O(VE²), тоді як загальний Форд-Фалкерсон може займати O(E × max_flow) кроків при невдалому виборі шляху.

Едмондс-Карп:
1. Побудувати залишковий граф G_f: для кожного ребра (u,v) з ємністю c і потоком f,
       додати пряме ребро з залишковою ємністю c - f
       додати зворотне ребро з залишковою ємністю f
2. Знайти найкоротший доповнювальний шлях s→t через BFS у G_f
3. Вузьке місце Δ = мінімальна залишкова ємність вздовж шляху
4. Доповнити: збільшити потік на кожному прямому ребрі на Δ, зменшити на зворотному на Δ
5. Повторювати, доки не залишиться доповнювальних шляхів (BFS не знаходить шляху s→t)

Складність: O(VE²) — щонайбільше VE/2 доповнень, кожен BFS за O(E)

Ідентифікація мінімального розрізу та анімація

Після завершення алгоритму мінімальний розріз ідентифікується BFS у залишковому графі, починаючи з s: усі вершини, досяжні з s, формують множину S; усі інші формують T. Кожне ребро, що перетинає межу від S до T на повній ємності, є ребром розрізу — намальоване червоним. Симуляція показує кожен крок доповнення BFS по одному: обраний шлях підсвічується жовтим, значення потоку на кожному ребрі оновлюються чисельно, а залишкові ємності анімуються. Після збіжності анімація мінімального розрізу малює межу розбиття та позначає її ємністю (яка дорівнює максимальному потоку).

• Пресет 1 — класичний підручниковий приклад; максимальний потік = мінімальний розріз = 23.
• Пресет 2 — ромбоподібна мережа; демонструє розріз розміром 3, що відокремлює джерело від стоку.
• Пресет 3 — трубопровід з вузьким ребром; єдине ребро розрізу тривіально є мінімальним розрізом.
• Пресет 4 — формулювання двочасткового узгодження; максимальний потік дає максимальний розмір узгодження.
• Пресет 5 — велика випадкова мережа; демонструє продуктивність з 20 вузлами та 40 ребрами.

🌈 Крижаний ореол — оптика атмосферних кристалів льоду методом Монте-Карло

Фізика крижаних ореолів

Атмосферні крижані ореоли викликані заломленням і відбиттям сонячного світла крізь грані гексагональних призм льоду (кристалічна структура Ih). Гексагональна призма має дві гексагональні базальні грані та шість прямокутних граней призми. Ключовий кут заломлення залежить від кута призми α між двома гранями та показника заломлення n льоду для видимого світла (n ≈ 1.31 при 550 нм, дисперсне від ~1.307 для червоного до ~1.317 для фіолетового).

Закон Снелла і мінімальне відхилення

Для променя, що входить у грань 1 і виходить з грані 2 призми з кутом клина α, кут відхилення δ:

При вході:   n_air · sin(θ_i) = n_ice · sin(θ_r)
При виході:    n_ice · sin(θ_r') = n_air · sin(θ_e)

Геометрія:   θ_r + θ_r' = α   (кут між гранями призми)
Повне відхилення: δ = θ_i + θ_e - α

Мінімальне відхилення (δ_min) відбувається при θ_i = θ_e (симетричний промінь):
  sin((δ_min + α)/2) = n · sin(α/2)

Для ореолу 22°: α = 60° (чергові грані призми), n = 1.31
  sin((δ_min + 60°)/2) = 1.31 · sin(30°) = 0.655
  (δ_min + 60°)/2 = arcsin(0.655) ≈ 40.9°
  δ_min ≈ 21.8° ≈ 22°

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

Трасування променів методом Монте-Карло

Замість аналітичного інтегрування розподілу відхилення симуляція випускає N = 100 000 променів на кадр, кожен крізь випадково орієнтований (чи обмежений) кристал льоду:

для кожного променя:
    1. Вибрати орієнтацію кристала:
         випадкова → усі промені ореолу 22°
         пластинка (нахил σ ≈ 1°) → бічне сонце / парагелій
         колонка (нахил σ ≈ 1°) → навколозенітна дуга
    2. Обчислити грань входу: перетнути промінь з геометрією призми
    3. Застосувати закон Снелла на вході → напрямок заломлення всередині кристала
    4. Поширити до грані виходу; застосувати закон Снелла знову
    5. Обчислити (азимут, висоту) вихідного променя
    6. Накопичити в 2D кутовій гістограмі (проєкція неба)
    7. Застосувати дисперсію довжини хвилі: повторити для λ = 400, 550, 700 нм

Отримана гістограма нормалізується й рендериться як зображення неба з кутовою системою координат. Ореол 22° з'являється як дифузне яскраве кільце на кутовому радіусі 22° від симульованого сонця, яскравіше на внутрішньому краю (мінімальне відхилення). Червоний на внутрішній стороні кільця (менше n → більше δ_min), синій на зовнішній — протилежно до райдуги.

Бічні сонця та навколозенітна дуга

Бічні сонця (парагелії) викликані кристалами-пластинками, що осідають з вертикальною c-віссю через аеродинамічний опір. Промені входять крізь вертикальну грань призми та виходять крізь іншу вертикальну грань під кутом 60°, але тепер лише кристали з правильним азимутальним кутом дають внесок, концентруючи світло на ±22° від сонця вздовж горизонту:

Умова парагелію: світло входить/виходить крізь вертикальні грані призми
  Нахил кристала від горизонталі: σ ≈ 0°–2° (осідання пластинок)
  Кутова позиція: при сонячній висоті h бічне сонце зміщене на
    Δaz = arctan(tan(22°)/cos(h))
  При h = 0°: Δaz = 22°;  при h = 30°: Δaz ≈ 26°

Навколозенітна дуга (CZA) утворюється, коли світло входить крізь верхню базальну грань і виходить крізь вертикальну грань призми кристалів-пластинок. Геометрія дає вигнуту дугу на 46° над сонцем, з червоним на зовнішній (нижній) стороні та фіолетовим на внутрішній — утворюючи найінтенсивніше забарвлену дугу в атмосферній оптиці, часто насиченішу, ніж райдуга.

Що далі

Хвиля 64 у розробці з кандидатами, що охоплюють аеродинаміку коефіцієнта опору, геометрію геодезичного купола, стабільне узгодження (Гейл-Шеплі), картини інтерференції голограми, алгоритм виявлення спільнот Лувена та структуру зерна / кристалографію Вороного. Бібліотека на шляху перевищити 600 симуляцій до літа.