Девлог #87 – Хвиля 66: Епіцикли Фур'є, магнітний маятник, змочування поверхні, зубчаста передача, рюкзак і Гаффман

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

🌀 Епіцикли Фур'є — малювання за допомогою DFT

Від намальованого шляху до комплексного DFT

Коли ви малюєте замкнену криву на полотні, симулятор вибирає з неї N рівновіддалених значень параметра, розглядаючи кожну точку вибірки як комплексне число z_n = x_n + i·y_n. Обчислюється комплексне дискретне перетворення Фур'є:

X[k] = Σ_{n=0}^{N-1}  z_n · exp(−2πi·k·n / N)    для k = 0, 1, …, N−1

Кожен бін k представляє один епіцикл:
  центральна частота : k  обертів за повний період малювання
  амплітуда        : |X[k]| / N
  початкова фаза    : arg(X[k])

Реконструкція (часткова сума M епіциклів):
  z(t) = (1/N) · Σ_{k=0}^{M-1}  X[k] · exp(2πi·k·t / N)

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

Чим це відрізняється від симулятора рядів Фур'є

Існуючий симулятор fourier-series платформи застосовує 1D DFT до періодичних хвильових форм (прямокутної, пилкоподібної, трикутної) і візуалізує амплітуди гармонік на стовпчастій діаграмі. Симулятор fourier-epicycles використовує 2D комплексний DFT на довільному намальованому шляху — і x(t), і y(t) кодуються разом як єдиний комплексний сигнал. Це означає, що кожен частотний бін дає внесок у вигляді 2D кругового руху (епіцикла), а не просто вертикального коливання. Набір обертових кіл — це візуальний доказ того, що будь-яку неперервну замкнену криву можна розкласти на кругові рухи.

Вбудовані пресети та їхні спектри

П'ять пресетів обрані для ілюстрації різних спектральних профілів:

• Коло — єдиний ненульовий бін при k = 1; одного епіцикла достатньо для ідеальної реконструкції.
• Зірка (5-променева) — домінантні гармоніки при k = 1, 4, 6, 9, 11, … (кратні 5 ± 1); дає характерні п'ятикутні спалахи.
• Серце — асиметрична форма вимагає багатьох гармонік; повільна збіжність розкриває загострення внизу.
• Вісімка — симетрична; лише непарні гармоніки присутні в x-компоненті, даючи чистий спектр.
• Лемніската (∞) — параметризована як x = cos t / (1 + sin²t), y = sin t cos t / (1 + sin²t); домінує k = 2 через подвійну петлю.

Рендеринг і продуктивність

Обчислення повного DFT при N = 512 коштує O(N²) множень. Симулятор використовує алгоритм швидкого перетворення Фур'є Кулі-Тюкі з основою 2 (доповнення до наступного степеня двійки), щоб зменшити це до O(N log N):

// Ітеративне на місці FFT Кулі-Тюкі
function fft(x):   // x: комплексний масив, довжина має бути степенем 2
    N = x.length
    перестановка бітового реверсу x
    для s = 1 до log2(N):
        m = 2^s
        w_m = exp(-2πi / m)
        для k = 0 до N-1 крок m:
            w = 1
            для j = 0 до m/2 - 1:
                t = w * x[k + j + m/2]
                u = x[k + j]
                x[k + j]       = u + t
                x[k + j + m/2] = u - t
                w = w * w_m

Анімація епіциклів рендериться з частотою 60 fps, перераховуючи лише позиції кінчиків обертових плечей кожен кадр (без FFT на кадр — DFT обчислюється один раз при відпусканні миші). Слід шляху малюється на позаекранне полотно й компонується, зберігаючи основний цикл легким.

🧲 Магнітний маятник — хаос і фрактальні басейни

Рівняння руху

Тягарець маятника моделюється як точкова маса m, підвішена на нитці довжиною L, обмежена рухом у горизонтальній площині безпосередньо над трьома магнітами. Висота z тягарця над площиною магнітів утримується фіксованою на h (висота спокою). Рівняння руху в декартових координатах (x, y):

Сили на тягарець:
  Відновлювальна гравітаційна:  F_g = −(mg/L) · (x, y)   (наближення малого кута)
  Демпфування:            F_d = −b · (ẋ, ẏ)
  Магнітна (магніт i у позиції r_i = (x_i, y_i)):
    d_i = sqrt((x − x_i)² + (y − y_i)² + h²)
    F_m_i = μ · (r_i − r) / d_i³

Повне ДР:
  ẍ = −(g/L)·x  −  b·ẋ  +  Σ_i  μ·(x_i − x) / d_i³
  ÿ = −(g/L)·y  −  b·ẏ  +  Σ_i  μ·(y_i − y) / d_i³

Параметри за замовчуванням: L = 0.5 м, h = 0.05 м, b = 0.2 (демпфування), μ = 1.0 (магнітна сила). Три магніти розташовані під кутами 0°, 120°, 240° на колі радіусом 0.3 м.

Інтегрування RK4

Кожна траєкторія інтегрується за допомогою класичного методу Рунге-Кутти 4-го порядку з фіксованим кроком часу Δt = 0.01 с. Вектор стану s = (x, y, ẋ, ẏ):

k1 = f(t,       s)
k2 = f(t + h/2, s + h·k1/2)
k3 = f(t + h/2, s + h·k2/2)
k4 = f(t + h,   s + h·k3)
s_{n+1} = s_n + (h/6)·(k1 + 2k2 + 2k3 + k4)

Симуляція зупиняється, коли швидкість тягарця падає нижче 0.001 м/с протягом 50 послідовних кроків, записуючи, який магніт найближчий, як «захоплений» магніт. Для карти басейну обчислюється сітка початкових позицій (до 400×400 = 160 000 траєкторій) за допомогою Web Workers для порційного рендерингу, тож інтерфейс залишається чуйним.

Фрактальні межі басейнів

Отримана карта басейну — канонічний приклад фракталу. На кожній видимій межі між двома кольоровими регіонами є третій регіон, вкладений нескінченно глибоко — межі мають розмірність Хаусдорфа, строго більшу за 1. Це пряме наслідок чутливості до початкових умов (додатні показники Ляпунова) поблизу сепаратрис. Наближення до межі розкриває самоподібну тонку структуру на кожному масштабі, роблячи це однією з найбільш візуально доступних демонстрацій хаосу на платформі.

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

💧 Змочування поверхні — від рівняння Юнга до ефекту лотоса

Рівняння Юнга (ідеальна плоска поверхня)

Коли крапля рідини лежить на плоскій твердій поверхні в парогазовому середовищі, три поверхневі напруги діють на лінії контакту: тверде тіло-пара (γ_SV), тверде тіло-рідина (γ_SL) та рідина-пара (γ_LV). Механічна рівновага вздовж поверхні вимагає:

Рівняння Юнга:
  γ_SV = γ_SL + γ_LV · cos(θ_Y)

  → cos(θ_Y) = (γ_SV − γ_SL) / γ_LV

Режими кута контакту:
  θ_Y = 0°         : повне змочування (рідина розтікається нескінченно)
  0° < θ_Y < 90°   : гідрофільна (вода змочує поверхню)
  90° < θ_Y < 150° : гідрофобна (вода збирається в краплі)
  θ_Y > 150°        : супергідрофобна (ефект лотоса)

Анімований поперечний переріз краплі неперервно змінюється разом із повзунком кута контакту. Три стрілки поверхонь оновлюють напрямок і мітку, щоб показати баланс сил.

Модель Венцеля (шорстка поверхня)

Реальні поверхні ніколи не є ідеально плоскими. Венцель (1936) ввів коефіцієнт шорсткості r (відношення реальної площі поверхні до проєктованої площі, завжди r ≥ 1) в енергетичний баланс:

Рівняння Венцеля:
  cos(θ_W) = r · cos(θ_Y)

  Наслідок:
    Якщо θ_Y < 90°  (гідрофільна): шорсткість робить її більш змочуваною  → θ_W < θ_Y
    Якщо θ_Y > 90°  (гідрофобна): шорсткість робить її менш змочуваною  → θ_W > θ_Y

Приклади коефіцієнта шорсткості:
  Полірований метал:  r ≈ 1.05
  Піскоструминна обробка:    r ≈ 1.5 – 2.0
  Лист лотоса:      r ≈ 2.5 (подвійна текстура мікро- та наномасштабу)

Візуалізація шорсткості поверхні показує збільшений поперечний переріз підкладки з регульованою амплітудою нерівностей і просторовою частотою. Кут контакту краплі оновлюється в реальному часі, коли змінюються θ_Y і r.

Модель Кассі-Бакстера та ефект лотоса

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

Рівняння Кассі-Бакстера:
  cos(θ_CB) = f₁ · cos(θ₁) + f₂ · cos(θ₂)

Для композиту тверде тіло-повітря (f₂ = частка повітря, θ₂ = 180°):
  cos(θ_CB) = f_s · cos(θ_Y) − (1 − f_s)

  де f_s = дробова площа твердого контакту

Лист лотоса: f_s ≈ 0.05 (лише 5% твердого контакту)
  cos(θ_CB) ≈ 0.05 · cos(105°) − 0.95 ≈ −0.963
  θ_CB ≈ 164°  (супергідрофобна)

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

⚙️ Зубчаста передача — передавальні числа, крутний момент і складені ступені

Основи передавального числа

Дві зачеплені зовнішні циліндричні шестерні з кількістю зубців N₁ і N₂ обертаються з кутовими швидкостями ω₁ і ω₂. Оскільки їхні початкові кола котяться без ковзання, тангенціальна швидкість у точці зачеплення рівна:

Радіуси початкових кіл: r₁ = N₁·m/2,  r₂ = N₂·m/2  (m = модуль = розмір зубця)

Умова зачеплення без ковзання: r₁·ω₁ = r₂·ω₂

  → передавальне число: ω₂/ω₁ = N₁/N₂  (величина)

Напрямок: зовнішнє зачеплення змінює обертання → ω₂/ω₁ = −N₁/N₂
           внутрішнє зачеплення (кільцева шестерня) зберігає обертання

Передавальне число крутного моменту (ідеальне, без тертя):
  T₂/T₁ = N₂/N₁   (крутний момент множиться на той самий коефіцієнт, на який ділиться швидкість)

Холості шестерні

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

Ведуча (N₁) → Холоста (N_i) → Ведена (N₂):
  ω_i / ω₁ = −N₁/N_i       (перше зачеплення, перевернуте)
  ω₂ / ω_i = −N_i/N₂       (друге зачеплення, знову перевернуте)
  ω₂ / ω₁  = N₁/N₂         (кількість зубців холостої скасовується)

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

Складені зубчасті передачі

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

3-ступенева складена передача:
  Ступінь 1: ведуча N₁ зачіплюється з шестернею N₂ (на валу A)
  Ступінь 2: шестерня N₃ (на валу A) зачіплюється з шестернею N₄ (на валу B)
  Ступінь 3: шестерня N₅ (на валу B) зачіплюється з виходом N₆

  Загальне передавальне число: ω₆/ω₁ = (N₁·N₃·N₅) / (N₂·N₄·N₆)

Приклад: N₁=10, N₂=40, N₃=12, N₄=60, N₅=15, N₆=75
  Передавальне число = (10·12·15)/(40·60·75) = 1800/180000 = 1/100
  Вхідна швидкість 1000 об/хв → вихід 10 об/хв, крутний момент ×100

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

🎒 Рюкзак — таблиця ДП 0/1 відкрита

Визначення задачі та складність грубої сили

Дано n предметів, кожен з вагою w_i та цінністю v_i, і рюкзак місткістю W, знайти підмножину предметів з максимальною сумарною цінністю за умови, що загальна вага ≤ W. Кожен предмет можна включити не більше одного разу (обмеження «0/1»). Підхід грубої сили перебирає всі 2ⁿ підмножини — доцільно лише для n ≤ 20 приблизно.

Динамічне програмування знизу вгору

Рекурентність ДП заповнює таблицю (n+1) × (W+1), де K[i][w] = максимальна цінність, досяжна з використанням перших i предметів у рюкзаку місткістю w:

Базовий випадок:
  K[0][w] = 0   для всіх w   (без предметів → нульова цінність)

Рекурентність (для i = 1…n, w = 0…W):
  якщо w_i > w:
      K[i][w] = K[i-1][w]            // предмет i занадто важкий: пропустити його
  інакше:
      K[i][w] = max(
          K[i-1][w],                 // не брати предмет i
          v_i + K[i-1][w - w_i]     // взяти предмет i
      )

Відповідь: K[n][W]
Складність: O(n·W) за часом, O(n·W) за пам'яттю (можна звести до O(W) з ковзним масивом)

Симуляція рендерить повну таблицю, кольорово закодовану за величиною цінності. Коли анімація заповнює кожну клітину, вона підсвічує дві порівнювані клітини (K[i-1][w] і K[i-1][w−w_i]) коротким спалахом, роблячи структуру рекурентності одразу видимою.

Зворотне відстеження для пошуку обраних предметів

Коли таблиця завершена, фактично обрані предмети відновлюються зворотним відстеженням від K[n][W]:

i = n, w = W
поки i > 0 і w > 0:
    якщо K[i][w] != K[i-1][w]:
        включити предмет i           // клітина змінилася → предмет i був узятий
        w = w − w_i
    i = i − 1

Обрані предмети підсвічуються золотою рамкою у списку предметів. Панель порівняння поруч запускає жадібний алгоритм (сортувати за відношенням цінність/вага, жадібно заповнювати) і позначає будь-яку розбіжність — для рюкзака 0/1 жадібний алгоритм часто пропускає оптимум.

Чому ДП не є справді поліноміальним

Складність O(nW) іноді називають «псевдополіноміальною», тому що W може бути експоненційно великим у кількості бітів, що використовуються для його кодування. Якщо W = 2^k, таблиця має n·2^k клітин — експоненційно від розміру входу. Справжні поліноміальні за часом алгоритми для рюкзака 0/1 невідомі (це NP-важка задача). Симуляція робить це конкретним за допомогою відображення «довжини бітів»: посуньте W від 100 до 100 000 і спостерігайте, як розмір таблиці (і час рендерингу) масштабується лінійно з W, але експоненційно з log W.

🌳 Кодування Гаффмана — оптимальне стиснення без префіксів

Побудова дерева з мін-купою

Дано алфавіт символів з відомими частотами (чи ймовірностями), алгоритм Гаффмана будує оптимальний бінарний префіксно-вільний код:

Вхід: символи s_1…s_k з частотами f_1…f_k

Алгоритм:
1. Вставити всі символи в чергу з мінімальним пріоритетом (мін-купу), ключ за частотою.
2. Поки розмір купи > 1:
     а. Витягнути два вузли з мінімальною частотою: лівий (f_L), правий (f_R)
     б. Створити внутрішній вузол з частотою f_L + f_R
     в. Вставити внутрішній вузол назад у купу
3. Вузол, що залишився — корінь дерева Гаффмана.

Призначення кодів:
  Обійти дерево від кореня: рух ліворуч додає '0', праворуч — '1'.
  Код кожного листа = бітовий рядок на шляху від кореня до цього листа.

Симуляція анімує кожен крок вилучення та об'єднання. Зростаюче дерево малюється з мітками частот на кожному вузлі, а два вилучені вузли світяться, коли об'єднуються.

Ентропія та середня довжина коду

Теорема Шеннона про кодування джерела стверджує, що жоден однозначно декодований код не може досягти середньої довжини коду, коротшої за ентропію H джерела:

Ентропія Шеннона:
  H = −Σ_{i=1}^{k}  p_i · log₂(p_i)    (бітів на символ)

Середня довжина коду Гаффмана:
  L̄ = Σ_{i=1}^{k}  p_i · len(code_i)

Гарантія:  H ≤ L̄ < H + 1

Для оптимального блочного коду довжиною n символів: L̄ → H при n → ∞

Живий графік на правій панелі наносить H (горизонтальна опорна лінія) та L̄ (стовпчик) поруч. Коли ви редагуєте вхідний текст, обидва значення оновлюються, а розрив між ними позначається в мілібітах. Для природного англійського тексту H ≈ 4.1 біта/символ, а Гаффман досягає L̄ ≈ 4.2 біта/символ — у межах приблизно 2% теоретичного мінімуму, використовуючи однозначні коди.

Доказ властивості префіксної свободи

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

Що далі

Хвиля 67 уже побудована і виводить бібліотеку на знакову віху: позначку 600 симуляцій. Шість доповнень охоплюють теорію хаосу (діаграма біфуркації логістичного відображення зі сталою Фейгенбаума), оптику (муарові візерунки й наліз Найквіста), обчислювальну геометрію (тріангуляція Делоне), моделювання транспорту (клітинний автомат трафіку Нагеля-Шреккенберга), алгоритми (просторове розбиття quadtree) та орбітальну механіку (приливні сили з циклами весняних/квадратурних припливів). Читайте про це все в Девлозі #88.