Девлог #86 – Хвиля 65: Капілярність, A*, параметричні криві, вихровий пучок, колоїди та маршрутизація в інтернеті

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

💧 Капілярність — закон Жюрена і меніск

Поверхневий натяг і кут контакту

Поверхневий натяг γ (Н/м) виникає через асиметричні міжмолекулярні сили на межі рідина-повітря: молекули на поверхні мають менше сусідів, ніж молекули об'єму, і перебувають у стані з вищою енергією. Коли рідина контактує з твердим тілом, конкурують три поверхневі енергії: рідина-повітря (γ_LA), тверде тіло-рідина (γ_SL) та тверде тіло-повітря (γ_SA). Їхній баланс визначає кут контакту θ через рівняння Юнга:

Рівняння Юнга (баланс сил на лінії контакту):
  γ_SA = γ_SL + γ_LA · cos θ
  cos θ = (γ_SA − γ_SL) / γ_LA

Класи кута контакту:
  θ = 0°:    ідеальне змочування (рідина розтікається повністю)
  0 < θ < 90°: змочування (увігнутий меніск, капілярне підняття)
  θ = 90°:   нейтральний
  90° < θ < 180°: незмочування (випуклий меніск, капілярне пригнічення)
  θ = 180°:  ідеальне незмочування (супергідрофобний стан Кассі-Бакстера)

Вода на склі:  θ ≈ 20°  (змочування)
Ртуть на склі: θ ≈ 140° (незмочування)

Закон Жюрена: капілярне підняття

У трубці радіусом r меніск створює перепад тиску на вигнутій поверхні (рівняння Юнга-Лапласа). Для сферичного меніска в круглій трубці цей перепад тиску підтримує стовп рідини висотою h:

Тиск Юнга-Лапласа на вигнутій межі:
  ΔP = 2γ / R_c
  де R_c = r / cos θ   (радіус кривизни меніска в трубці)

Баланс тиску (капілярний тиск = гідростатичний тиск):
  2γ cos θ / r = ρ g h

Закон Жюрена:
  h = 2γ cos θ / (ρ g r)

Приклад — вода в склі (r = 0.5 мм):
  γ = 0.073 Н/м, θ = 20°, ρ = 1000 кг/м³, g = 9.81 м/с²
  h = 2 × 0.073 × cos(20°) / (1000 × 9.81 × 0.0005) ≈ 28 мм

Симуляція рендерить поруч трубки різних радіусів з анімованим підняттям рідини. Форма меніска малюється як дуга кола з кривизною 1/R_c. Живий графік h проти 1/r підтверджує лінійну залежність, передбачену законом Жюрена. Для ртуті кут контакту θ > 90° дає cosθ < 0, тож h від'ємне — ртуть опускається нижче рівня резервуара.

Температурна залежність поверхневого натягу

Поверхневий натяг спадає приблизно лінійно з температурою для більшості рідин. Для води γ(T) можна наблизити правилом Етвеша, і вона табульована від 0°C (0.0756 Н/м) до 100°C (0.0589 Н/м). Повзунок температури симуляції регулює γ, використовуючи реальні дані по воді, і перераховує висоту підняття в усіх трубках одночасно, демонструючи, що нагріта вода піднімається нижче — актуально для капілярних явищ при кипінні та в мікрофлюїдиці.

🗺️ Пошук шляху A* — евристики, оптимальність і відкритий набір

Алгоритм A*

A* (Гарт, Нільссон, Рафаель, 1968) шукає у зваженому графі найкоротший шлях від старту s до цілі g, оцінюючи кожен вузол n за оцінкою вартості f(n) = g(n) + h(n):

f(n) = g(n) + h(n)

  g(n) = точна вартість найкоротшого знайденого шляху від s до n
  h(n) = евристична оцінка вартості від n до цілі g

A* підтримує дві множини:
  Відкритий набір:   виявлені, але ще не розгорнуті вузли (мін-купа за значенням f)
  Закритий набір: уже розгорнуті вузли

Основний цикл:
  поки відкритий набір не порожній:
    n = вузол у відкритому наборі з найнижчим f(n)   // O(log N) з бінарною купою
    якщо n == ціль: повернути шлях
    додати n до закритого набору
    для кожного сусіда m вузла n:
      tentative_g = g(n) + edge_cost(n, m)
      якщо m у закритому наборі і tentative_g ≥ g(m): продовжити
      якщо tentative_g < g(m):
        g(m) = tentative_g
        f(m) = g(m) + h(m)
        parent(m) = n
        додати m до відкритого набору (або зменшити ключ, якщо вже присутній)

Допустимість евристики й оптимальність

A* гарантовано знаходить оптимальний шлях тоді й лише тоді, коли евристика h(n) є допустимою: вона ніколи не переоцінює справжню вартість досягнення цілі. Три реалізовані евристики:

Координати сітки (x, y) для клітини, ціль у (gx, gy):

Манхеттенська відстань (норма L1):
  h(n) = |x − gx| + |y − gy|
  Допустима для 4-зв'язних сіток (без діагонального руху)
  Недопустима для 8-зв'язних сіток (недооцінює діагональні скорочення)

Евклідова відстань (норма L2):
  h(n) = √((x − gx)² + (y − gy)²)
  Допустима для будь-якого руху; має тенденцію розгортати менше вузлів, ніж
  Манхеттенська на 8-зв'язних сітках, але вимагає порівнянь з плаваючою комою

Діагональна (Чебишева) відстань:
  h(n) = max(|x − gx|, |y − gy|)
  Допустима для 8-зв'язних сіток з одиничною діагональною вартістю
  Точна для рівномірного діагонального руху; розгортає найменше вузлів

Порівняння з Дейкстрою та жадібним найкращим-першим

Панель порівняння алгоритмів запускає всі три алгоритми одночасно на тому самому лабіринті й показує кількість розгорнутих вузлів поруч:

• Дейкстра (h = 0): розгортається концентричними кільцями; гарантовано оптимальний; найбільше розгорнутих вузлів.
• Жадібний найкращий-перший (f = лише h, ігнорує g): кидається до цілі; найменше розгорнутих вузлів; не оптимальний — може знайти довший шлях.
• A* з Манхеттенською/Евклідовою/Діагональною: оптимальний і зазвичай розгортає набагато менше вузлів, ніж Дейкстра; Діагональна евристика найближча до ідеалу на 8-зв'язних сітках.

Реалізація на бінарній купі тримає кожну ітерацію на рівні O(log N), де N — розмір відкритого набору, роблячи A* практичним на сітках до 500×500 (250 000 вузлів) з інтерактивною швидкістю. Лічильник розгортання вузлів оновлюється в реальному часі на кожному кадрі покрокової анімації.

✏️ Параметричні криві — Без'є, B-сплайни та NURBS

Криві Без'є та побудова де Кастельжо

Крива Без'є степеня n визначається n+1 контрольними точками P₀, …, P_n і параметризована t ∈ [0, 1]. Алгоритм де Кастельжо обчислює будь-яку точку на кривій повторною лінійною інтерполяцією:

Алгоритм де Кастельжо:
  P_i^(0) = Pᵢ             (ініціалізація контрольними точками)
  P_i^(r) = (1−t) · P_i^(r−1) + t · P_{i+1}^(r−1)   для r = 1, …, n

  Точка кривої при параметрі t: B(t) = P_0^(n)

Представлення в базисі Бернштейна (еквівалентне):
  B(t) = Σᵢ₌₀ⁿ  C(n,i) · tⁱ · (1−t)^(n−i) · Pᵢ
  де C(n,i) = n! / (i! (n−i)!)  (біноміальний коефіцієнт)

Степінь 3 (кубічна) Без'є — найпоширеніша:
  B(t) = (1−t)³P₀ + 3t(1−t)²P₁ + 3t²(1−t)P₂ + t³P₃

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

B-сплайни та рекурентність Кокса-де Бура

B-сплайни узагальнюють криві Без'є з локальним контролем: переміщення однієї контрольної точки впливає лише на криву в локальній околиці, а не на всю криву. Вони визначаються вектором вузлів T = {t₀, …, t_m} і базисними функціями N_i,k(t), обчисленими за Коксом-де Буром:

Рекурентність Кокса-де Бура:
  N_{i,0}(t) = 1  якщо tᵢ ≤ t < t_{i+1},  інакше 0
  N_{i,k}(t) = (t − tᵢ)/(t_{i+k} − tᵢ) · N_{i,k−1}(t)
             + (t_{i+k+1} − t)/(t_{i+k+1} − t_{i+1}) · N_{i+1,k−1}(t)

Крива B-сплайну степеня k з n+1 контрольними точками:
  C(t) = Σᵢ₌₀ⁿ  N_{i,k}(t) · Pᵢ

Властивості:
  Локальна підтримка: N_{i,k}(t) ≠ 0 лише на [tᵢ, t_{i+k+1}]
  Розбиття одиниці: Σᵢ N_{i,k}(t) = 1 для всіх t
  Неперервність C^(k−1) скрізь (C^(k−2) на повторних вузлах)

Редактор вектора вузлів дозволяє перетягувати окремі значення вузлів або вставляти/видаляти вузли. Повторні вузли локально знижують неперервність — вузол, повторений k разів, дає кут при цьому значенні параметра, корисно для загострень і різких деталей. Симуляція кольорово кодує кожну базисну функцію N_i,k(t) окремо й наносить усі їх на смугу під основним полотном.

NURBS і точне коло

Нерівномірні раціональні B-сплайни (NURBS) додають вагу w_i до кожної контрольної точки, забезпечуючи точне представлення конічних перерізів — включно з колами та дугами — які поліноміальні сплайни можуть лише наблизити:

Крива NURBS:
  C(t) = Σᵢ N_{i,k}(t) · wᵢ · Pᵢ
         ─────────────────────────
         Σᵢ N_{i,k}(t) · wᵢ

Точне чверть-коло з 3 контрольними точками (степінь 2):
  P₀ = (1, 0),  w₀ = 1
  P₁ = (1, 1),  w₁ = √2/2 ≈ 0.7071   (вага проміжної точки)
  P₂ = (0, 1),  w₂ = 1

  Вектор вузлів: [0, 0, 0, 1, 1, 1]

Повне коло: 9 контрольних точок (три чверть-дуги, з'єднані з C¹ неперервністю)
             Ваги чергуються між 1 і √2/2

Панель NURBS показує контрольний многокутник з повзунками ваги для кожної точки. Зменшення ваги нижче 1 відтягує криву від цієї контрольної точки; збільшення понад 1 притягує криву до неї. Точне коло — стандартний приклад NURBS за замовчуванням — зміна центральної ваги з √2/2 на 1 деформує його в параболу, демонструючи роль ваг.

🌀 Вихровий пучок — моди Лагерра-Гаусса та орбітальний кутовий момент

Моди пучка Лагерра-Гаусса

Моди Лагерра-Гаусса (LG) — точні розв'язки параксіального хвильового рівняння в циліндричних координатах (ρ, φ, z). Вони утворюють ортогональний базис для пучків з обертальною симетрією. Пучок LG з радіальним індексом p та азимутальним індексом ℓ має амплітуду електричного поля:

LG_{p,ℓ}(ρ, φ, z):

  u(ρ, φ, z) = C_{p,ℓ} · (ρ√2/w(z))^|ℓ| · L_p^|ℓ|(2ρ²/w²(z))
              · exp(−ρ²/w²(z))                     ← гаусова обвідна
              · exp(i ℓ φ)                          ← спіральна фаза
              · exp(−i k ρ²/(2R(z)))               ← кривизна хвильового фронту
              · exp(i (2p+|ℓ|+1) ζ(z))            ← фаза Гоуі

де:
  w(z)  = талія пучка в позиції z: w(z) = w₀√(1 + (z/z_R)²)
  R(z)  = радіус кривизни: R(z) = z(1 + (z_R/z)²)
  ζ(z)  = фаза Гоуі: arctan(z/z_R)
  z_R   = довжина Релея: πw₀²/λ
  L_p^|ℓ| = узагальнений поліном Лагерра

Спіральна фаза й орбітальний кутовий момент

Азимутальний фазовий множник exp(iℓφ) намотує фазу на 2πℓ радіанів за один повний оберт по φ. Цей спіральний фронт хвилі означає, що пучок несе орбітальний кутовий момент (OAM) ℓℏ на фотон — відмінний від спінового кутового моменту ±ℏ циркулярно поляризованого світла:

OAM на фотон: L_z = ℓ ħ

Профіль інтенсивності: I ∝ |u|² ∝ ρ^(2|ℓ|) · [L_p^|ℓ|(…)]² · exp(−2ρ²/w²)

Для p = 0 (без радіальних вузлів):
  Одне яскраве кільце («пончик») при ρ_max = w(z) · √(|ℓ|/2)
  Темне ядро (фазова сингулярність) при ρ = 0

Для p > 0:
  p концентричних темних кілець усередині зовнішнього кільця
  Загалом p+1 яскравих кілець

Топологічний заряд ℓ:
  |ℓ| = кількість фазових намоток 2π за азимутальний оберт
  Знак визначає хіральність (за годинниковою стрілкою проти проти годинникової рампи фази)

Симуляція рендерить 2D-профіль інтенсивності як теплову карту та накладає карту фази як кольорове колесо на будь-якій обраній z-площині. Панель інтерференції показує характерні спіральні (одна щілина) або вилкоподібні (дві щілини) візерунки, що використовуються в лабораторії для вимірювання топологічного заряду ℓ. Коли пучок LG інтерферує з плоскою хвилею, отриманий візерунок має |ℓ| додаткових чи відсутніх смуг у точці вилки.

🧪 Колоїди — теорія DLVO та колоїдна стабільність

Броунівський рух і рівняння Стокса-Ейнштейна

Колоїдна частинка (радіус r, 1 нм-1 μм), суспендована в рідині, зазнає броунівського руху, тому що теплові флуктуації безперервно її штовхають. Коефіцієнт дифузії D пов'язує середньоквадратичне зміщення з часом і задається рівнянням Стокса-Ейнштейна:

Рівняння Стокса-Ейнштейна:
  D = k_B T / (6π η r)

де:
  k_B = 1.38 × 10⁻²³ Дж/К  (стала Больцмана)
  T   = температура (К)
  η   = динамічна в'язкість розчинника (Па·с)
  r   = радіус частинки (м)

Співвідношення Ейнштейна (дифузія-рухливість):
  D = μ k_B T   де μ = рухливість = 1/(6π η r)

Середньоквадратичне зміщення:
  ⟨r²⟩ = 6 D t   (3D броунівський рух)
  ⟨r²⟩ = 4 D t   (2D)

Приклад — частинка 100 нм у воді при 25°C:
  D = (1.38×10⁻²³ × 298) / (6π × 8.9×10⁻⁴ × 50×10⁻⁹) ≈ 4.9 × 10⁻¹² м²/с

Теорія DLVO: стабільність і флокуляція

Теорія DLVO (Дерягін, Ландау, Фервей, Овербек) описує повний потенціал взаємодії між двома зарядженими колоїдними частинками як суму електростатичного відштовхування та притягання ван дер Ваальса:

Повний потенціал DLVO:
  V_total(h) = V_EDL(h) + V_vdW(h)

Електростатичне відштовхування подвійного шару (лінеаризоване Пуассона-Больцмана):
  V_EDL(h) = 64π ε ε₀ r (k_B T / ze)² tanh²(zeψ₀/4k_BT) · exp(−κh)

  κ = довжина екранування Дебая-Хюккеля:
  κ⁻¹ = √(ε ε₀ k_B T / (2 N_A e² I))
  I = іонна сила = ½ Σ cᵢ zᵢ²

Притягання ван дер Ваальса (Гамакер):
  V_vdW(h) = −A_H / (12π h²) · f(r, h)
  A_H = стала Гамакера (типове значення ~10⁻²⁰ Дж для оксид/вода/оксид)

Бар'єр стабільності:
  Якщо V_max > ~15 k_B T: стабільний колоїд (кінетично загальмована флокуляція)
  Якщо V_max < ~5 k_B T:  швидка коагуляція (швидка агрегація Смолуховського)
  Додавання солі стискає довжину Дебая → знижує бар'єр → флокуляція

Повзунок іонної сили стискає довжину Дебая κ⁻¹ приблизно з 100 нм (низька сіль) до ~1 нм (фізіологічна), зменшуючи енергетичний бар'єр V_max і викликаючи флокуляцію. Симуляція візуалізує це живою кривою V(h) з висотою бар'єра, позначеною в одиницях k_BT. Коли бар'єр падає нижче ~5k_BT, частинки починають агрегувати на полотні, а гістограма розподілу розміру агрегатів оновлюється в реальному часі.

Седиментація

Частинки більші за ~1 μм осідають під дією гравітації. Швидкість седиментації для сфери в режимі Стокса:

Швидкість седиментації Стокса:
  v_s = 2r²(ρ_p − ρ_f) g / (9η)

Число Пекле (адвекція проти дифузії):
  Pe = v_s · r / D = v_s · 6π η r² / (k_B T)

Pe < 1: домінує дифузія → стабільна суспензія
Pe > 1: домінує седиментація → гравітаційне осідання

Приклад — сфера кремнезему 1 мкм у воді:
  ρ_p = 2200 кг/м³, ρ_f = 1000 кг/м³
  v_s = 2 × (500×10⁻⁹)² × 1200 × 9.81 / (9 × 8.9×10⁻⁴) ≈ 0.73 мкм/с
  D = 4.9×10⁻¹³ м²/с  →  Pe ≈ 1.5 (граничний випадок)

Симуляція забарвлює частинки за їхнім поточним розміром агрегатного кластера (одинак = синій, димер = зелений, …, великий кластер = червоний) і накладає вектори швидкості осідання для частинок, де Pe > 1.

🌐 Маршрутизація в інтернеті — вектор відстані, стан каналу та BGP

Маршрутизація вектора відстані Беллмана-Форда

Протоколи вектора відстані (наприклад, RIP) змушують кожен маршрутизатор періодично транслювати всю свою таблицю маршрутизації прямим сусідам. Кожен маршрутизатор оновлює свою таблицю, використовуючи рекурентність Беллмана-Форда:

Рекурентність Беллмана-Форда:
  d_x(y) = min по всіх сусідах v { c(x,v) + d_v(y) }

  d_x(y) = оцінена вартість від маршрутизатора x до призначення y
  c(x,v) = вартість прямого каналу x → v
  d_v(y) = оголошена сусідом v вартість до y

Правило оновлення (коли x отримує оновлення від сусіда v):
  для кожного призначення y:
    якщо c(x,v) + d_v(y) < d_x(y):
      d_x(y) = c(x,v) + d_v(y)
      next_hop_x(y) = v    // оновити таблицю маршрутизації

Збіжність: O(діаметр) раундів оновлення для коректної топології

Проблема рахування до нескінченності

Коли канал виходить з ладу, протоколи вектора відстані можуть впасти в патологічну петлю, де два маршрутизатори по черзі збільшують оголошену вартість до призначення, відскакуючи один від одного, доки відстань не «порахує до нескінченності» (зазвичай обмежено 16 у RIP):

Топологія: A — B — C  (канал B-C виходить з ладу; лишається лише канал A-B)

До відмови: d_B(C) = 1, d_A(C) = 2
Після відмови:
  B помічає, що канал B-C неактивний; d_B(C) = ∞
  Але таблиця B досі показує d_A(C) = 2 (застаріла)
  B оновлює: d_B(C) = 1 + d_A(C) = 3  (неправильно! через A, але A маршрутизує через B)
  A оновлює: d_A(C) = 1 + d_B(C) = 4
  B оновлює: d_B(C) = 5 … і так далі, аж до ∞

Виправлення — розділений горизонт:
  Маршрутизатор A ніколи не оголошує B маршрут, вивчений від B
  → оголошення A для C (вивчене через B) пригнічується для B
  → B коректно вивчає d_B(C) = ∞ після відмови каналу

Виправлення — отруєний зворотний зв'язок:
  Агресивніше: A явно оголошує d_A(C) = ∞ назад до B
  Усуває рахування до нескінченності для дводольних петель негайно

Симуляція дозволяє викликати відмову каналу, клацнувши будь-яке ребро, а потім пройти покроково через раунди оновлення Беллмана-Форда. З вимкненим розділеним горизонтом петля рахування до нескінченності чітко видима в таблицях маршрутизаторів. Увімкнення розділеного горизонту показує ту саму відмову, вирішену за один раунд.

Стан каналу Дейкстри (OSPF) та перезбіжність

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

Затоплення стану каналу:
  Кожен маршрутизатор генерує LSA: {ID маршрутизатора, [сусід, вартість, ...]}
  LSA затоплюється до всіх маршрутизаторів у зоні (номери послідовності запобігають петлям)
  Усі маршрутизатори будують ідентичну базу даних стану каналу (LSDB)

Дейкстра на LSDB (для кожного маршрутизатора):
  Ідентична Дейкстрі із самим маршрутизатором як джерелом
  Результат: повне дерево найкоротших шляхів до всіх призначень
  Встановлюється в таблиці маршрутизації як записи пересилання

Перезбіжність після відмови каналу:
  Маршрутизатор/канал, що вийшов з ладу, генерує нове LSA з нескінченною вартістю
  Затоплює мережу → всі маршрутизатори оновлюють LSDB → повторно запускають Дейкстру
  Збіжність за O(E log V) на маршрутизатор проти O(діаметр) раундів для BF

BGP (протокол граничного шлюзу):
  Протокол вектора шляху: оголошує повний шлях AS, а не лише відстань
  Усуває петлі маршрутизації: відхиляє оновлення, що містять власний AS у шляху
  На основі політики: оператори контролюють, які маршрути приймаються/пріоритетні

Анімація перезбіжності показує поширення потоку LSA як хвилю по діаграмі мережі, потім підсвічує кожен маршрутизатор, що повторно запускає Дейкстру, коли потік досягає його. Таблиці маршрутизаторів оновлюються візуально на панелі праворуч. Вкладка BGP демонструє механізм вектора шляху зі спрощеною топологією з трьох AS, показуючи, як атрибут AS-шляху запобігає проблемі рахування до нескінченності, яка мучить протоколи вектора відстані.

Що далі

Хвиля 66 готується з шістьма подальшими симуляціями: епіцикли Фур'є, що анімують будь-яку замкнену криву як суму обертових кіл, магнітний маятник з фрактальними басейнами притягання, змочування поверхні та гістерезис кута контакту, кінематика зубчастих передач, задача рюкзака 0/1, розв'язана динамічним програмуванням, та кодування Гаффмана з межами ентропії. Бібліотека наближається до 600 симуляцій.