неділю, 3 грудня 2017 р.

Розподіл порядкових статистик

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

В цьому дописі я розгляну лише випадкові вибірки отримані з неперервних розподілів (тобто функція розподілу ймовірностей неперервна). Нехай $X_1, \dots, X_n$ - це випадкова вибірка розміру $n$ з неперервного розподілу з функцією розподілу ймовірностей $F(x)$. Впорядкувавши цю вибірку у зростному порядку отримуємо $Y_1, \dots, Y_n$.

Порядкова статистика $Y_i$ називається $i$-тою порядковою статистикою. Тут ми припускаємо, що $Y_1 < Y_2 < \dots < Y_n$, тобто збігів бути не може.

Функція розподілу ймовірностей

Якщо має місце подія $Y_i \le y$ тоді ми знаємо, що щонайменше $i$ з $X_j$ лежать ліворуч від $y$. Розглянемо подію, що $Y_j \le y$ як успіх, тоді має місце $n$ проб Бернуллі з імовірністю успіху $F(X \le y)$. Нас цікавить імовірність мати не менше ніж $i$ успіхів: \begin{equation} F_{Y_i}(y) = P(Y_i \le y) = \sum_{k=i}^n\binom{n}{k}F(y)^k(1-F(y))^{n-k} \label{eq:cdf} \end{equation}

Густина ймовірності

Для знаходження густини ймовірності нам знадобиться така формула: \begin{equation} F_{Y_i}(y) = F_{Y_{i-1}}(y) - \binom{n}{i-1}F(y)^{i-1}(1-F(y))^{n-i+1} \label{eq:cdf_diff} \end{equation} Я одразу наведу формулу густину ймовірності, а потім ми доведемо її індукційно. \begin{equation} f_{Y_i}(y) = \frac{n!}{(i-1)!(n-i)!}F(y)^{i-1}(1-F(y))^{n-i}f_X(y)\label{eq:pdf} \end{equation} Почнемо з бази індукції. Розглянемо $i = 1$. Отже, нас цікавить іморність того, що хоча б одна проба менша ніж $y$. Доповняльною подією буде те, що усі проби більші ніж $y$. Тому, $F_{Y_1}(y) = (1 - F(y))^n$. Для отримання густини ймовірності нам потрібно обчислити похідну: $$f_{Y_1}(y) = n(1 - F(y))^{n-1}f_X(y)$$ Припустимо, що ми довели \ref{eq:pdf} для $Y_{i-1}$: $$f_{Y_{i-1}}(y) = \frac{n!}{(i-2)!(n-i+1)!}F(y)^{i-2}(1-F(y))^{n-i+1}f_X(y)$$ Тепер розглянемо випадок для $Y_i$. Тут нам можна скористатись рівнянням \ref{eq:cdf_diff}: \begin{align*} f_{Y_i}(y) &= f_{Y_{i-1}}(y)\\ &\phantom{=}-\binom{n}{i-1}(i-1)F(y)^{i-2}f_X(y)(1-F(y))^{n-i+1}\\ &\phantom{=}-\binom{n}{i-1}F(y)^{i-1}(n-i+1)(1-F(y))^{n-i}f_X(y) \end{align*} Трошки алгебри і правий бік рівняння перетвориться в \ref{eq:pdf}.

Коментар

Розглянемо другу половину \ref{eq:pdf}: $$F(y)^{i-1}f_X(y)(1-F(y))^{n-i}$$ тут перший множник - це імовірність, що $i-1$ проб менші ніж $y$, другий - імовірність, що одна проба поблизу $y$ і третій - імовірність, що $n-i$ проб більші ніж $y$. Таким чином \ref{eq:pdf} - це такий мультиноміальний розподіл: \begin{equation} f_{Y_i}(y) = \frac{n!}{(i-1)!1!(n-i)!}F(y)^{i-1}f_X(y)(1-F(y))^{n-i}\label{eq:pdf_multi} \end{equation}

четвер, 2 листопада 2017 р.

Двоїстість у лінійному програмуванні

Згадаємо формулювання прямої і двоїстої програм у стандартній формі.
\begin{align} \text{максимізувати}\ &\vec{c}^T\vec{x}\\ \text{за умов}\ &A\vec{x}\le \vec{b}\\ &\vec{x}\ge 0 \end{align} \begin{align} \text{мінімізувати}\ &\vec{b}^T\vec{y}\\ \text{за умов}\ &A^T\vec{y} \ge \vec{c} \label{eq:dual_constr}\\ &\vec{y}\ge 0 \end{align}

Фізичне тлумачення

  • Застосувати поле $\vec{c}$ до кульки в середині обмеженого багатогранника $A\vec{x} \le \vec{b}$.
  • Зрештою, кулька зупиниться біля краю багатогранника (в точці $\vec{x}'$).
  • Грань $\vec{a}_i\vec{x}'\le b_i$ накладає певну силу $-y_i\vec{a}_i$ на кульку.
  • Кулька стабільна, тому $\vec{c}^T = \sum_i y'_i\vec{a}_i = \vec{y}'^TA$.
  • Двоїстість можна уявити як намагання мінімізувати "роботу" $\sum_i y_ib_i$, щоб перенести м'ячик у початок координат рухаючи багатогранник.
  • В оптимальній позиції тиснуть лише стіни суміжні з кулькою. (На рисунку лише $y_1, y_2$ не нульові.)
Звернемо увагу на те, що на кульку тиснуть лише грані до яких вона притуляється. Тобто для кожного $i$, або вона перебуває на лінії $\vec{a}_i^T \vec{x} = b_i$ або тиск дорівнює нулю: \begin{align} \vec{a}_i\vec{x} = b_i, y_i > 0&\ \text{торкається грані, тиск є} \label{eq:compslack1}\\ \vec{a}_i\vec{x} < b_i, y_i = 0&\ \text{не торкається грані, тиску немає}\label{eq:compslack2} \end{align} Тепер покажемо, що наш $\vec{x}'$ і є оптимальним розв'язком. Згідно з принципом слабкої двоїстості $\vec{c}^T\vec{x} \le \vec{b}^T\vec{y}$. Отже, якщо ми покажемо, що $\vec{c}^T\vec{x}' = \vec{b}^T\vec{y}'$, то це і є оптимальний розв'язок. Cтабільність кульки разом із третім законом Ньютона даюсть нам, що $A^T \vec{y}' = \vec{c}$. Завдяки рівнянням (\ref{eq:compslack1}, \ref{eq:compslack2}) ми можемо записати, що $\vec{a}_i\vec{x}'y'_i = b_iy'_i$, а отже $$\vec{c}^T\vec{x}' = \vec{x}'^T\vec{c} = \vec{x}'^TA^T\vec{y}' = (A\vec{x}')^T\vec{y}' = \sum_i\vec{a}_i\vec{x}'y'_i = \sum_i b_iy'_i = \vec{b}^T\vec{y}'$$ Отже, $\vec{x}'$ - оптимальний розв'язок, а $\vec{y}'$ - двоїстий до нього.

Економічне тлумачення

Розглянемо задачу оптимального виробництва:
  • $n$ продуктів, $m$ сирих матеріалів.
  • Продукт $j$ потребує $a_{ij}$ одиниць сирого матеріалу $i$.
  • Всього доступно $b_i$ одиниць матеріалу $i$.
  • Продукт $j$ дає $c_j$ прибутку на одиницю.
  • Майстерня хоче максимізувати прибуток за наявних матеріалів
\begin{matrix} \text{максимізувати}&\vec{c}^T\vec{x}&\\ \text{за умов}&\vec{a}_i^T\vec{x}\le b_i,&\ i = 1, \dots, m.\\ &x_j\ge 0,& j = 1,\dots,n. \end{matrix} Якщо задача розв'язна і розв'язок скінченний, то допустима область - це опуклий багатогранник. Завжди існує щонайменше одна вершина в якій цільова функція набуває оптимального значення. Для задачі оптимального виробництва це означає, що для $n$ продуктів і $m$ сирих матеріалів, існує оптимальний план з щонайбільше $m$ продуктами.

Тлумачення двоїстих змінних

Для кожної одиниці продукту $j, j \in \{1,n\}$, нам потріно $a_{ij}, i \in \{1, m\}$. Для того, щоб було вигідно випускати продукт, його ціна має бути не меншою ніж сумарна ціна всіх складових. Позначимо ціну продукту через $c_p(j)$ і ціну матеріалу через $c_m(i)$, тоді $$c_p(j) \ge \sum_i a_{ij} c_m(i)$$ Але ж з точністю до знаку це те саме, що і обмеження двоїстої програми (\ref{eq:dual_constr}). Поглянемо на це з іншого боку. За умови що $\vec{x}$ і $\vec{y}$ оптимальні розв'язки відповідно прямої і двоїстої задач згідно з принципом сильної двоїстості маємо: $$\vec{c}^T\vec{x}=\vec{b}^T\vec{y}$$
  • Ліворуч маємо максимальний виторг.
  • Праворуч - $\sum_i (\text{наявність матеріалу}\ i) \times (\text{виторг на одиницю матеріалу}\ i)$
Інакше кажучі: Значення $y_i$ в оптимумі - це двоїста ціна ресурсу $i$.

Якщо ж значення не оптимальні, тобто

$$\vec{c}^T\vec{x}<\vec{b}^T\vec{y}$$ Розв'язок не оптимальний, бо ресурси не використані повністю.

суботу, 7 жовтня 2017 р.

Найкоротший шлях у графі. Рішення через прямий-двоїстий алгоритм.

Розглянемо добре знайому задачу знаходження найкоротшого шляху. Маючи зв'язний граф $G = (V, E)$, цінову функцію $\omega: E \to R_+$ і дві вершини $s$ і $t$, знайти шлях найменшої вартості, що сполучає ці дві вершини в $G$. Цю задачу можна ефективно розв'язати за допомогою алгоритма Дейкстри. Але тут ми побудуємо прямий-двоїстий алгоритм, що ефективно обчислить найкоротший шлях.

Введемо такі позначення: $\mathcal S = \{ S \mid S \subset V, s \in S, t \notin S\}$ і, для кожної множини $S \subset V$, через $\delta(S)$ позначимо множину ребер, що мають одну вешину в $S$, а другу в $V\setminus S$.

Для цього ми роглянемо лінійну програму в формі. Взагалі-то, нам потрібно, щоб $x_e$ набували лише ${0,1}$, але цілочисельні лінійні програми NP-складні, тому ми ослаблюємо програму до нецілочельної лінійної програми. \begin{align} \min \sum_{e\in E}x_e\cdot\omega(e) \end{align} за умови, що \begin{align} \forall S \in \mathcal{S}, \sum_{e\in \delta(S)} x_e\ge 1\\ \forall e \in E, x_e \ge 0 \end{align} Наступний крок - це створення двоїстої лінійної програми. Для того, щоб це зробити спочатку згадаємо як виглядає лінійна програма в стандартній формі: \begin{align*} &\max \sum_{j=1}^n c_j x_j\\ \text{за умови, що}&\sum_{j=1}^n a_{ij}x_j \le b_i, i = 1,2,\dots, m\\ &x_j \ge 0, j = 1,2, \dots, n \end{align*} Часто зручніше записувати в матричному вигляді: \begin{align*} &\max c^Tx\\ \text{за умови, що}&\ Ax \le b, x \ge 0 \end{align*} Тоді двоїста програма така: \begin{align*} &\min b^T y\\ \text{за умови, що}&\ A^Ty \ge c, y \ge 0 \end{align*} Тепер ми можемо записати двоїсту програму:
\begin{align} \max \sum_{S\in\mathcal{S}}y_S \end{align} за умови, що \begin{align} \sum_{S:e\in \delta(S)} y_S\le c_e, \forall e \in E \label{eq:dual_constraint}\\ \forall S \in \mathcal{S}, y_S \ge 0 \end{align} Тобто, тепер ми максимізуємо суму нових зміних $y_S$. У (\ref{eq:dual_constraint}) ми для кожного ребра, сумуємо $y_S$ з усіх $S$, що перетинають це ребро.

Алгоритм

Роглянемо такий алгоритм:
  1. $F\gets \varnothing$
  2. $y \gets 0$
  3. Допоки в $F$ ще нема шляху, що поєднує $s$ і $t$
    • Нехай $C$ - компонента зв'язності графу $G' = (V, F)$, що містить $s$
    • Збільшити двоїсту змінну $y_C$ так, щоб вперше виконалось $\sum\limits_{S \in \mathcal{S}~:~e \in \delta(S)} y_S = w(e')$
    • Додати $e'$ до $F$
  4. Повернути шлях $P$ з $F$, що поєднує $s$ і $t$.

Коректність

На кожній ітерації ми збільшуємо іншу $y_S$, бо $C$ постійно змінюється. Отже виконання алгоритму завершиться не пізніше ніж через $|V|$ кроків.

Перебіг виконання

Спочатку сформулюємо лему, яка нам допоможе зрозуміти як виконується алгоритм.
Лема
На кожному кроці алгоритму множина $F$ - це дерево.
Важливо, що ми повертаємо не все дерево, а лише шлях, що поєднує дві вершини. Читач може переконатись, що порядок додавання ребер в дереву той самий, що й у алгоритму Дейкстри.

суботу, 23 вересня 2017 р.

Геометричний погляд на перцептрони

Двопогоровий нейрон

  • Спочатку обчислити зважену суму входових даних від інших нейронів (плюс зсув).
  • Тоді видати 1, якщо зважена сума більше нуля.
\begin{equation}\label{eq:z_def} z = b + \sum_i x_i w_i \end{equation} \begin{equation}\label{eq:y_def} y = \begin{cases} 1, \mbox{якщо}\ z \ge 0\\ 0, \mbox{інакше} \end{cases} \end{equation} Тут ми не зважатимемо на зсув, бо ми можемо розглядати зсув як вагу для ще одного входового елементу, що постійно дає 1.

Простір-ваг.

  • Цей простір має один вимір на кожну вагу.
  • Вектор (точка) в цьому просторі представляє якийсь набір ваг.
Припускаючи, що ми не зважаємо на зсуви, ми можемо розглядати кожен тренувальний випадок як гіперплощину крізь початок координат. Ваги мусять лежати з одного боку цієї гіперплощини, щоб відповідь була правильною.
Погляньмо, що ж це означає матиематично. припустимо, що у нас на вході є вектор з координатами $\vec{v} = (a, b)$. Щоб цей вектор повертав правильну відповідь, згідно з \ref{eq:z_def} і \ref{eq:y_def} нам потрібно, щоб $\vec{v}\cdot\vec{w} \ge 0$. Тут $\vec{w}$ - це вектор ваг з координатами $(x ,y)$. Тобто, нам треба, щоб $ax+by \ge 0$. Якщо переписати нерівність як $y = -\frac{a}{b} x$, то стає зрозуміло, що підхожі набори ваг лежать з боку цієї гіперплощини куди вказує входовий вектор.

Конус допустимих розв'язків

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

Процедура навчання

Вибираємо тренувальні випадки за алгоритмом згідно з яким зрештою ми переберемо усі тренувальні випадки.
  • Якщо вихід правильний, залишаємо ваги без змін.
  • Якщо результат - це помилковий 0, додаємо входовий вектор до ваг.
  • Якщо результат - це помилкова 1, віднімаємо входовий вектор від ваг.
Ця процедура гарантує знаходження набору ваг, який видає правильний резльтат для всіх тренувальних випадків, якщо такий набір існує.

Чому процедура працює

Розглянемо піднесену до квадрату відстань $d_a^2+d_b^2$ між будь-яким допустимим вектором ваг і поточним вектором ваг. Нам би хотілось щоб кожного разу коли перцептрон помиляється ми наближали поточний вектор ваг до всіх допустимих векторів ваг. Але це не так. Тому ми розглянемо щедро допустимі вектори ваг, що лежать в допустимомій області із відступом не менше ніж довжина входового вектора, який визначає гіперплощину-обмеження.

Кілька слів про доведення

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

Приклад програми, що реалізує цей алгоритм навчання.

Тут я подаю приклад програми, яка вчиться розрізняти лінійно роздільну множину. Насправді більшість коду - це виведення графіки в gif-файл.
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
ds = [1 1; 2 -2; -1 -1.5; -2 -1; -2 1; 3 -1];
class = [1 -1 -1 -1 1 1];

ds_neg = ds(class==-1,:);
ds_pos = ds(class==1,:);

w = [0 1 0.5]';
eta = 0.2;

count = size(ds, 1);

y = @(w, x) -(w(2)*x+w(1))/w(3);

h = figure;
epoch = 0;
% тут ми припускаємо, що множини лінійно роздільні, тому
% ми зупинимось лише коли не буде помилок
stillMistakes = true;
while stillMistakes
    for i = 1:count 
        x = [1 ds(i, :)];
        if class(i) * x * w < 0
            w = w + class(i) * eta * x';
        end

        hold off
        plot([-5, 5], [y(w, -5), y(w, 5)]);
        title(strcat('епоха', {' '}, num2str(epoch)));
        hold on
        scatter(ds_neg(:,1), ds_neg(:,2), 'o');
        scatter(ds_pos(:,1), ds_pos(:,2), '+');
        plot(x(2),x(3), 'o', 'MarkerSize',15);
        axis([-5, 5, -5, 5]);

        drawnow
        frame = getframe(1);
        im{epoch * count + i} = frame2im(frame);
    end
    epoch = epoch + 1;
    
    % як множину для перевірки ми використовуємо усю тренувальну множину
    stillMistakes = false;
    for i = 1:count
        x = [1 ds(i, :)];
        if class(i) * x * w < 0
            stillMistakes = true;
            break
        end
    end
end
close;

fn = 'linearlySeparable.gif';
for idx = 1:(epoch * count)
    [A,map] = rgb2ind(im{idx},256);
    if idx == 1
        imwrite(A,map,fn,'gif','LoopCount',Inf,'DelayTime',1);
    else
        imwrite(A,map,fn,'gif','WriteMode','append','DelayTime',1);
    end
end

суботу, 9 вересня 2017 р.

Нотатки на курс "Вступ до комплексного аналізу"

Аналітичність

Теорема
Нехай $f = u + iv$ визначена в області $D \subset C$. Тоді $f$ - аналітична в $D$ ТТТ, якщо $u(x, y)$ і $v(x, y)$ мають неперервні перші частинні похідні в $D$, які задовольняють рівнянням Коші-Рімана.

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

Означення
Нехай $D \subset C$ буде областю визначення і нехай $f : D \to C$ буде неперервною функцією. Первісна $f$ у $D$ - це аналітична функція $F : D \to C$ така, що $F' = f$ у $D$.
Теорема
Якщо $f$ - це неперервна функція в області $D$ і, якщо $f$ має первісну $F$ в $D$, тоді для будь-якої кривої $\gamma : [a, b] \to D$ маємо, що $$\int_{\gamma}f(z)dz = F(\gamma(b)) − F(\gamma(a)).$$

Коли $f$ має первісну?

Теорема (Коші для трикутників)
Нехай $D$ буде відкритою множиною в $C$ і нехай $f$ буде аналітичною в $D$. Нехай $T$ буде трикутником, що вміщається в $D$ (включно з границею) і нехай $\delta T$ буде його границею, орієнтованою позитивно. Тоді $$\int_{\delta T}f(z)dz = 0.$$
Теорема (Морери)
Якщо $f$ неперервна у однозв'язній області $D$ і, якщо $\int_{\gamma}f(z)dz = 0$ для кожної трикутної кривої в $D$, тоді $f$ має первісну в $D$.
Теорема (Ґурсата)
Нехай $D$ буде однозв'язною областю в $C$, і нехай $f$ буде аналітичною в $D$. Тоді $f$ має первісну в $D$. Більше того, первісна задається явно вибором точки $z_0 \in D$ і покладенням $$F(z) = \int_{z_0}^z f(w)dw,$$ де інтеграл береться по довільній кривій в $D$ від $z_0$ до $z$.
Теорема (Коші для однозв'язних областей)
Нехай $D$ буде однозв'язною областю в $C$ і нехай $f$ - аналітична в $D$. Нехай $\gamma : [a, b] \to D$ буде кусково гладкою, замкнутою кривою в $D$ (тобто $\gamma(b) = \gamma(a)$). Тоді $$\int_{\gamma}f(z)dz = 0.$$
Наслідок
Нехай $\gamma_1$ і $\gamma_2$ - це дві прості замкнуті криві (тобто жодна з них не перетинає саму себе), орінтовані проти годинникової стрілки, при чому $\gamma_2$ всередині $\gamma_1$. Якщо $f$ аналітична в області $D$, яка містить обидві криві і область між ними, тоді $$\int_{\gamma_1}f(z)dz =\int_{\gamma_2}f(z)dz.$$
Теорема (інтегральна формула Коші)
Нехай $D$ - це однозв'язна область, обмежена кусково гладкою кривою $\gamma$, і нехай $f$ аналітична на множині $U$, що містить в собі замкнення $D$ (тобто $D$ і $\gamma$). Тоді $$f(w) = \frac{1}{2\pi i}\int_{\gamma} \frac{f(z)}{z-w}dz$$ для всіх $w \in D.$
Теорема
Якщо $f$ аналітична у відкритій множині $U$, тоді $f'$ також аналітична в $U$.
Теорема (інтегральна формула Коші для похідних)
Нехай $D$ - це однозв'язна область, обмежена кусково гладкою кривою $\gamma$, і нехай $f$ аналітична на множині $U$, що містить в собі замкнення $D$ (тобто $D$ і $\gamma$). Тоді $$f(w)^{(k)} = \frac{k!}{2\pi i}\int_{\gamma}\frac{f(z)}{(z − w)^{k+1}}dz$$ для всіх $w \in D, k ≥ 0$.
Це дуже захоплююча теорема, бо для обчислення значення $f$ або її похідної в будь-якій точці з області обмеженої кривою $\gamma$ нам необхідно знати лише значення функції на самій кривій.
Теорема (оцінка Коші)
Припустимо, що $f$ аналітична у відкритій множині, що містить $B_r(z_0)$, і що $|f(z)| \le m$ виконується на $\delta B_r(z_0)$ для деякої сталої $m$. Тоді для всіх $k \ge 0$, $$|f^{(k)}(z_0)| \le \frac{k!m}{r^k}.$$
Теорема (Ліувілля)
Нехай $f$ аналітична на всій комплексній площині. Якщо $f$ обмежена, тоді $f$ мусить бути сталою.
З цієї теореми випливає, що раз $\sin z$ аналітична на всій комплексній площині і вона не костантна, то значить, що вона сягає $\infty$ в якомусь напрямку. І дійсно, такий напрямок може бути $ni$.
Теорема (Принцип максимума)
Нехай $f$ аналітична в $D$ і припустимо, що існує точка $z_0 \in D$ така, що $|f(z)| \le |f(z_0)|$ для всіх $z \in D$. Тоді $f$ константна в $D$.
Наслідок
Якщо $D \subset C$ - це обмежена область і якщо $f : D \to C$ неперервна і аналітична в $D$, тоді $|f|$ досягає максимума на $\delta D$.

Ряди

Теорема (Про розкладення в ряд Лорана)
Якщо $f:U\to\mathbb{C}$ аналітична і $\{r < |z-z_0|< R\} \subset U$, тоді $f$ можна розкласти в ряд Лорана: $$f(z) = \sum_{k=-\infty}^{\infty}a_k(z-z_0)^k,$$ який збігається в кожній точці кільця і збігається абсолютно і рівномірно в кожному підкільці $\{s\le |z-z_0|\le t\}$, де $ r < s < t < R $.
Теорема
Якщо $f$ аналітична в $\{r < |z-z_0| < R\}$, тоді $$f(z) = \sum_{k=-\infty}^{\infty}a_k(z-z_0)^k,$$ де $$a_k = \frac{1}{2\pi i}\int_{|z-z_0|=s}\frac{f(z)}{(z-z_0)^{k+1}}dz$$ для будь-якого $s$ між $r$ і $R$ і всіх $k\in \mathbb{Z}$.
Це не видається дуже корисним для знаходження значень $a_k$, але це може допомогти оцінити значення деяких з них.
Означення
Точка $z_0$ - ізольована сингулярність $f$, якщо $f$ - аналітична в проколотому диску $\{0 < |z − z_0| < r\}$ з центром $z_0$.
Якщо $f$ має ізольовану сингулярність в $z_0$, тоді $f$ має розклад Лорана. Члени ряду з від'ємними степенями називаються головною частиною, а з додатніми - правильною частиною ряду. Можливі три типи поведінки:
  1. Відсутність від'ємних сетепенів $z$: $f(z) = \frac{\cos z - 1}{z^2} = \frac{1}{z^2}\left(-\frac{z^2}{2!}+\frac{z^4}{4!}-+\cdots\right)$.
  2. Скінченна кількість від'ємних степенів $z$: $f(z) = \frac{\cos z}{z^4} = \frac{1}{z^4}\left(1-\frac{z^2}{2!}+\frac{z^4}{4!}-+\cdots\right)$.
  3. Нескінченна кількість від'ємних степенів $z$: $f(z) = \cos\left(\frac{1}{z}\right) = 1 - \frac{1}{2!} \frac{1}{z^2} + \frac{1}{4!}\frac{1}{z^4} - \frac{1}{6!}\frac{1}{z^6} +- \cdots$
Означення
Припустимо, що $z_0$ це ізольована сингулярність аналітичної функції $f$ з рядом Лорана $$\sum^{\infty}_{k=−\infty} a_k (z − z_0)^k, 0 < |z − z_0| < r.$$ Тоді сингулярність $z_0$ є
  • усувною, якщо $a_k = 0$ для всіх $k < 0$.
  • полюсом, якщо існує $N > 0$ таке, що $a_{−N} \neq 0$, але $a_k = 0$ для всіх $k < −N$. Індекс $N$ - порядок полюса.
  • істотною, якщо $a_k \neq 0$ для нескінченної кількості $k < 0$.
Теорема (Касораті-Веєрштраса)
Припустимо, що $z_0$ це істотна сингулярність $f$. Тоді для кожної точки $w_0 \in \mathbb{C}$ існує послідовність $\{z_n\}$ із $z_n \to z_0$ така, що $f(z_n) \to w_0$ коли $n \to ∞$.
Означення
Якщо $f$ має ізольовану сингулярність в $z_0$ з таким рядом Лорана $$f(z) = \sum^{\infty}_{k=−\infty} a_k (z − z_0)^k, 0 < |z − z_0| < r$$, тоді лишок $f$ в $z_0$ це $\operatorname{Res}(f, z0) = a_{-1}$.
Теорема (про лишки)
Нехай $D$ буде однозв'язною областю і нехай $f$ аналітична $D$, окрім як в ізольованих сингулярностях. Нехай $C$ - проста замкнута крива в $D$ (орієнтована проти годинникової стрілки), і нехай $z_1, \dots , z_n$ - ізольовані сингулярності $f$, які оточені $C$. Тоді $$\int_C f(z)dz = 2\pi i\sum_{k=1}^n \operatorname{Res}(f, z_k)$$.
  • Лишок в усувній сингулярності дорінює $0$.
  • Лишок в полюсі можна знайти за формулою $$\operatorname{Res}(f, z_0) = a_{−1} = \frac{1}{(n − 1)!} \lim_{z\to z_0}\frac{d^{n−1}}{dz^{n−1}}\left((z − z_0)^nf(z)\right).$$

неділю, 30 липня 2017 р.

Підсилення і комплексне підсилення

Формула експоненційного відгуку (ФЕВ)

Припустимо, що ми маємо рівняння $P(D)x = e^{rt}.$ Давайте спробуємо експоненційний розв'язок $x = A e^{rt}:$ $$e^{rt} = P(D)x = P(D)(Ae^{rt}) = AP(r)e^{rt}.$$ Отже, $A = 1/P(r)$ і розв'язком буде $$x_p(t) = \frac{e^{rt}}{P(r)}.$$ Це виконується якщо $P(r) \ne 0.$

Комплексна заміна.

ФЕВ чарівно поєднується з формулою Ейлера. Якщо коефіцієнти многочлена $P(r)$ дійсні, то $$\operatorname{Re}(P(D)f(t)) = P(D)\operatorname{Re}(f(t)).$$ З цього випливає, що розв'язок лінійного диф. рівняння $$P(D)x = \cos (\omega t)$$ зі сталими дійсними коефіцієнтами це дійсна частина розв'язку іншого рівняння $$P(D)z = e^{i\omega t}.$$ Це називається комплексною заміною.

Узагальнена ФЕВ

Випадок коли $P(r) = 0.$ Гляньмо ще раз на головне рівняння використане ФЕВ $$P(D)e^{rt} = P(r)e^{rt}.$$ Використаємо як змінну $r$ і продиференціюємо щодо $r$. $$ \begin{align*} P(D)te^{rt} &= P'(r)e^{rt} + P(r)re^{rt} &&\\ &= P'(r)e^{rt} && \mbox{бо} \ P(r) = 0 \end{align*} $$ Раз це виконується, то $\frac{te^{rt}}{P'(r)}$ це розв'язок $P(D)x = e^{rt}.$

Ми можемо продовжити процес: Якщо $P(r) = \cdots = P^{(k-1)} = 0$, але $P^{(k)}\ne 0,$ тоді

\begin{equation}\frac{t^ke^{rt}}{P^{(k)}(r)}\label{gen_erf}\end{equation} це розв'язок $$P(D)x = e^{rt}.$$

Комплексне підсилення

Поєднання комплексної заміни і ФЕВ надає критичну інформацію про синусоїдальний розв'язок рівнянь у вигляді $$P(D)x = Q(D)(\mbox{синусоїда}).$$ Цей метод дає синусоїдальний розв'язок у вигляді $$x_p(t) = \operatorname{Re}(Ge^{i\omega t})$$ де $\omega$ - це кутова швидкість входової синусоїди, а $G$ - це деяка комплексна стала.

Якщо виразити $G$ в полярній формі:

$$G=|G|e^{-i\phi },\quad \mbox{отже}\quad \phi =-\mathrm{arg}(G).$$ Тоді дійсна частина становитиме: $$|G|\cos(\omega t - \phi).$$ У випадку рівняння $$P(D)z = Q(D)e^{i\omega t},$$ знаходимо, що $$G(\omega)=\frac{Q(i\omega)}{P(i\omega)}.$$ Це комплексне число називається комплексним підсиленням. Воно містить в собі підсилення і фазове запізнення.

Трошки термінології

$\omega_n$ - природна частота. У випадку коли входова частота збігається з природною, то кажуть, що система в чистому резонансі, і тоді $\omega_n$ - це резонантна частота.

Приклад використання ФЕВ і перетворення Лапласа

Розглянемо рівняння: $$\ddot x + 9x = 9\cos(3t).$$ Для розв'язання ми можемо використати комплексне підсилення разом ФЕВ або перетворення Лапласа. Для використання перетворення Лапласа нам необхідно задати початкові умови. Нехай $x(0) = 2, \dot x(0) = 0$. Проведемо комлексну заміну, тобто ми шукатимемо дійсну частину розв'язку такого рівняння: $$\ddot x + 9x = 9e^{i3t}$$ Тут $P(D) = D^2 + 9$, отже маємо, що $P(i3) = 0.$ Тому нам необхідно використати узагальнену ФЕВ. Використавши \ref{gen_erf}, маємо: $$x=Re\left(\frac{9te^{i3t}}{6i}\right) = \frac{3}{2}t\sin(3t).$$ Для задоволення початкових умов, нам потрібен розв'язок $a\cos(3t)+b\sin(3t)$ відповідного однорідного рівняння $\ddot x + 9x = 0.$ Із нашими початковими умовами $a = 2, b = 0.$ І повний розв'язок такий: $\boxed{x = \frac{3}{2}t\sin(3t) + 2\cos(3t)}.$

Тепер використаємо перетворення Лапласа. Спочатку застосуємо перетворення Лапласа до обох боків рівняння. Таким чином ми перейдемо від диференціальної до алгебричної задачі. \begin{align*} \mathcal{L}(9\cos(3t)) &= \frac{9s}{s^2 + 9}\\ \mathcal{L}(x) &= X\\ \mathcal{L}(\ddot x) &= s^2X - sx(0) - \dot x(0)= s^2X - 2s\\ \end{align*} Тепер ми маємо алгебричне рівняння від $X$: \begin{align*} s^2X - 2s + 9X &= \frac{9s}{s^2+9}\\ X(s^2+9) &= \frac{9s}{s^2+9} +2s\\ X &= \frac{9s}{(s^2+9)^2} + \frac{2s}{s^2+9}\\ \end{align*} Зараз нам треба знайти зворотнє перетворення Лапласа для $X$. Із другим доданком все просто, він з точністю до сталої - табличний. $$\mathcal L^{-1}(\frac{2s}{s^2+9}) = 2\cos(3t)$$ З першим доданком трошки цікавіше. $$\mathcal L^{-1}\left(\frac{9s}{(s^2+9)^2}\right) = L^{-1}\left(\frac{3}{2}\left(\frac{3}{s^2+9}\right)'\right) = \frac{3}{2}t\sin(3t)$$ Отже, $\boxed{x = \frac{3}{2}t\sin(3t) + 2\cos(3t)}.$
Так, обидва методи дали ту саму відповідь, що ми й очікували.

Ремарка на перетворення Лапласа

Звернімо увагу, що зворотне перетворення Лапласа не однозначне, бо інтервал інтегування це $0, +\infty$. Тому частина, що менше або рівна нулю може бути будь-якою. Через це
Якщо задача починається зараз і продовжується в майбутньому і вам не потрібно нічого знати про минуле, це задача на перетворення Лапласа. Якщо вам треба знати про минуле, то це задача на перетворення Фур'є.

вівторок, 27 червня 2017 р.

Векторна тотожність $\vec{u}\times \vec{\omega}=\vec{\nabla}(\frac{1}{2}\vec{u}\cdot \vec{u})-\vec{u}\cdot \vec{\nabla}\vec{u}.$

У цьому дописі я розгляну векторну тотожність $$\vec{u}\times \vec{\omega}=\vec{\nabla}(\frac{1}{2}\vec{u}\cdot \vec{u})-\vec{u}\cdot \vec{\nabla}\vec{u}.$$ Почнемо з того, що $\vec\omega = \frac{1}{2}\vec\nabla\times\vec u,$ бо абсолютна величина кутової швидкості дорівнює половині абсолютної величини ротора. Отже, ми можемо переписати $$\vec{u}\times (\vec\nabla\times\vec{u})=\vec{\nabla}(\frac{1}{2}\vec{u}\cdot \vec{u})-\vec{u}\cdot \vec{\nabla}\vec{u}.$$ Або, скориставшись антикомутативністю векторного добутку: \begin{align} \vec\nabla\times\vec{u}\times\vec{u}=\vec{u}\cdot \vec{\nabla}\vec{u}-\vec{\nabla}(\frac{1}{2}\vec{u}\cdot \vec{u})\label{prepared}. \end{align}

Допоміжні поняття

Тут нам знадобиться поняття символу перестановки: $$ e_{ijk} = \begin{cases} +1 & \text{якщо } i,j,k - \text{це парна перестановка } 1,2,3 \\ -1 & \text{якщо } i,j,k - \text{це непарна перестановка } 1,2,3 \\ 0 & \text{інакше} \end{cases} $$ Розглянемо матрицю $A$ з елементами $a_{ir}$. За допомогою нашого символу ми можемо записати її визначник як \begin{align} \det A &= \frac{1}{3!}e_{ijk}e_{rst}a_{ir}a_{js}a_{kt}\nonumber,\\ \det A &= e_{rst}a_{1r}a_{2s}a_{3t}\label{detordered}. \end{align}

Тепер зауважимо, що \begin{align} e_{ijm}e_{rsm} = e_{1ij}e_{1rs}-e_{i2j}e_{r2s}+e_{ij3}e_{rs3}=e_{ij}e_{rs}\label{deltacancel}. \end{align} І в свою чергу \begin{align} e_{ij}e_{rs} = \delta_{ir}\delta_{js} - \delta_{is}\delta_{jr}, \end{align} де $\delta_{ij} -$ це символ Кронекера.

Доведення

Повернемось до \ref{prepared} і запишемо її у координатному вигляді. Для цього нам потрібно розібратись із $\vec{\nabla}\times\vec u.$ Використовуючи \ref{detordered} ми можемо записати: \begin{align} \vec{\nabla}\times\vec u &= e^{ijk}\vec{e}_i\frac{\partial}{\partial x_j}u_k.\label{curl} \end{align} \begin{align} \vec\nabla\times\vec{u}\times\vec{u} &= e_{lmn}\vec{e}_l(e_{mjk}\frac{\partial}{\partial x_j}u_k)u_n\nonumber\\ &= -\vec{e}_le_{lnm}e_{jkm}\frac{\partial}{\partial x_j}u_ku_n\nonumber\\ &= -\vec{e}_le_{ln}e_{jk}\frac{\partial}{\partial x_j}u_ku_n&&\text{скористалися }\ref{deltacancel}\nonumber\\ &= -\vec{e}_l\left(\delta_{lj}\delta_{nk}\frac{\partial}{\partial x_j}u_ku_n-\delta_{lk}\delta_{nj}\frac{\partial}{\partial x_j}u_ku_n\right)\nonumber\\ &= -\vec{e}_l\left(\frac{\partial}{\partial x_l}u_ku_k-\frac{\partial}{\partial x_j}u_lu_j\right)\label{proofongoing} \end{align} Тут ми трошки пригальмуємо і роззирнемось. Щоб рухатись далі згадаємо, що таке градієнт вектора: $$\vec\nabla\vec{u} = \frac{\partial u_i}{\partial x_j} \vec{e}_i\otimes\vec{e}_j.$$ Тепер нам варто розглянути як в індексному записі ми можемо представити величини, що ми зустріли: \begin{align*} \vec\nabla\vec{u} &= \frac{\partial u_i}{\partial x_j},&& ij\text{-й елемент}\\ \vec{\nabla}(\frac{1}{2}\vec{u}\cdot \vec{u}) &= \frac{1}{2}\vec\nabla\vec{u}_i\vec{u}_i + \frac{1}{2}\vec{u}_i\vec\nabla\vec{u}_i = \frac{\partial u_i}{\partial x_k}u_i, && k\text{-й елемент}\\ \vec{u}\cdot\vec\nabla\vec{u} &= (\vec{u}\cdot\vec\nabla)\vec{u} = u_i\frac{\partial u_k}{\partial x_i} && k\text{-й елемент} \end{align*} Використаємо нові знання про індексні записи у \ref{proofongoing}: \begin{align} \vec\nabla\times\vec{u}\times\vec{u} &= -\vec{e}_l\left(\frac{\partial u_k}{\partial x_l}u_k-u_j\frac{\partial u_l}{\partial x_j}\right)\nonumber\\ &= -\vec{\nabla}(\frac{1}{2}\vec{u}\cdot \vec{u}) + \vec{u}\cdot\vec\nabla\vec{u}\nonumber \end{align} Це і є те, що ми збирались довести. $\square$

четвер, 22 червня 2017 р.

Баланс мас у механіці суцільних середовищ

Баланс -- це питання швидкості зміни густини маси. У цьому дописі ми розглянемо рівняння балансу мас у механіці суцільних середовищ. Тут ми припускаємо, що якщо ми відкотимо деформацію, то ми отримаємо ту саму масу.

Розглянемо те як змінюється густина по мірі деформації тіла. Важливим припущенням для нас є те, що у густина -- це неперервна функція, і вона залишається неперервною незалежно від сили деформації, тобто ми не розглядаємо можливого утворення дірок.
Рисунок 1. Перетворення околу точки $\mathbf X,$ який ми позначатимемо $\mathcal N (\mathbf X)$ у окіл точки $\mathbf x$ -- $\mathcal n (\mathbf x).$
Для цього зосередимось на околі точки $\mathbf X$ у початковій конфігурації тіла - $\Omega_0.$ З плином часу цей окіл деформується в окіл точки $\mathbf x = \Phi(\mathbf X, t).$

Зауважимо, що $\rho_0(\mathbf X) = \lim_{\operatorname{vol}({\mathcal N(\mathbf X))}\to 0}\frac{\operatorname{m}(\mathcal N(\mathbf X))}{\operatorname{vol}(\mathcal N(\mathbf X))}.$

Тепер розглянемо $$ \begin{align*} \rho(\mathbf x, t) &= \lim_{\operatorname{vol}({\mathcal n(\mathbf x))}\to 0}\frac{\operatorname{m}(\mathcal n(\mathbf x))}{\operatorname{vol}(\mathcal n(\mathbf x))}\\ &=\lim_{\operatorname{vol}({\mathcal N(\mathbf x))}\to 0}\frac{\operatorname{m}(\mathcal n(\mathbf x))}{\operatorname{vol}(\mathcal N(\mathbf X))}\frac{\operatorname{vol}(\mathcal N(\mathbf X))}{\operatorname{vol}(\mathcal n(\mathbf x))}\\ &=\rho_0(\mathbf X) \det (\mathbf F^{-1}). \end{align*} $$ Тут ми використали градієнт деформації $F = \frac{\partial \mathbf x}{\partial \mathbf X},$ за допомогою якого ми можемо перейти від початкового об'єму до об'єму в час $ t$. А саме, ми знаємо, що $\operatorname{vol}(\mathcal n(\mathbf X)) = \det (\mathbf F) \operatorname{vol}(\mathcal N(\mathbf X))$. $$ \begin{equation} {\rho(\mathbf x, t)\Big\vert\,}_{\mathbf x = \Phi(\mathbf X, t)} = \rho_0(\mathbf X)\det \left(\mathbf F^{-1}(\mathbf X, t)\right) = \rho_0(\mathbf X) \mathcal J^{-1}(\mathbf X, t). \label{eq:rho} \end{equation} $$ Відображаючи умову вказану на початку допису, а саме, того що густина початкової конфігурації не залежить від часу, маємо: $$ \begin{equation} \frac{\partial \rho_0(\mathbf X)}{\partial t} = 0. \end{equation} $$ Як щодо $\rho (\mathbf x, t)$? $$\rho_0(\mathbf X) = {\rho(\mathbf x, t)\Big\vert\,}_{\mathbf x = \Phi(\mathbf X, t)} \mathcal J (\mathbf X, t).$$ $$ \begin{align} {\frac{\partial \rho_0(\mathbf X)}{\partial t}\Bigg\vert}_{\mathbf X - стала} &= \frac{\mathrm D \rho(\mathbf x, t)}{\mathrm D t} \mathcal J (\mathbf X, t) + \rho(\mathbf x, t) {\frac{\partial \mathcal J (\mathbf X, t)}{\partial t}\Bigg\vert}_{\mathbf X - стала}\nonumber\\ &=\left(\frac{\partial \rho(\mathbf x, t)}{\partial t} + \nabla \rho(\mathbf x, t) \cdot \mathbf{v}(\mathbf x, t)\right)\mathcal J (\mathbf X, t) + \rho(\mathbf x, t)\dot{\mathcal J}(\mathbf X, t),\label{eq:rho0dert} \end{align} $$ Тут $\frac{\mathrm D \rho(\mathbf x, t)}{\mathrm D t}$ позначає матеріальну похідну. Тепер розглянемо $\dot{\mathcal J}(\mathbf X, t)$ : $$ \begin{align} \dot{\mathcal J}(\mathbf X, t) &= \frac{\partial \det\left(\mathbf F(\mathbf X, t)\right)}{\partial t}\nonumber\\ &= \frac{\partial \det(\mathbf F(\mathbf X, t))}{\partial \mathbf F(\mathbf X, t)} \mathbf{:} \frac{\partial \mathbf F(\mathbf X, t)}{\partial t}\nonumber\\ &= \frac{\partial \det(\mathbf F(\mathbf X, t))}{\partial F(\mathbf X, t)_{iI}} \frac{\partial F(\mathbf X, t)_{iI}}{\partial t}\label{eq:dotIpre}. \end{align} $$ Тут значок $\mathbf{:}$ позначає згортку тензорів. Градієнт деформації має індекси, що складаються з великої і малої літери бо це двоточковий тензор, тобто він слугує для зв'язку між початковою і поточною конфігураціями. З цього моменту, усюди де це само собою зрозуміло я не писатиму параметри функції. $$\frac{\partial \det(\mathbf F)}{\partial \mathbf F_{iI}} = \frac{\partial \frac{1}{6}\varepsilon_{abc}\varepsilon_{ABC}F_{aA}F_{bB}F_{cC}}{\partial F_{iI}} = \operatorname{Cof}(\mathbf F)_{iI}$$ Де $\varepsilon_{ijk}$ -- це символ Леві-Чивіти, а $\operatorname{Cof}(\mathbf F)_{iI}$ -- це кофактор.

Згадаємо формулу для оберненої матриці: $A^{-1}_{ij}=\frac{\operatorname{Cof}(A^T_{ij})}{\det{A}},$ або в матричному вигляді: $A^{-1} = \frac{\operatorname{Cof}A^T}{\det A}$.

Отже, $\ref{eq:dotIpre}$ можна записати як: $$\dot{\mathcal J} = \mathcal J \mathbf F^{-T} \mathbf{:} \dot{\mathbf F} = \mathcal{J} \dot{F}_{iI} F^{-1}_{Ii} = \mathcal{J} \operatorname{trace}(\dot{\mathbf F} \mathbf F^{-1}) = \mathcal J \mathbb{I}\mathbf{:}(\dot{\mathbf F} \mathbf F^{-1}).$$ Для того, щоб рухатись далі нам потрібно розглянути градієнти матеріальної і просторової швидкостей.

Градієнт матеріальної швидкості $\mathbf V$: $\frac{\partial}{\partial \mathbf X}\frac{\partial \Phi}{\partial t} = \frac{\partial}{\partial t}\frac{\partial \Phi}{\partial \mathbf X} = \frac{\partial}{\partial t} \mathbf F = \dot{\mathbf F}.$

Градієнт просторової швидкості $\mathbf v$: ${\frac{\partial}{\partial \mathbf x}\frac{\partial \Phi}{\partial t}\Bigg\vert}_{\mathbf x = \Phi(\mathbf X, t)} = \frac{\partial}{\partial t}\frac{\partial \Phi}{\partial \mathbf X}\frac{\partial \mathbf X}{\partial \mathbf x} = \dot{\mathbf F}\mathbf F^{-1}.$

Введемо такі позначення: $$ \begin{align} \operatorname{GRAD}\mathbf V &= \frac{\partial \mathbf V}{\partial \mathbf X} = \dot{\mathbf F},\nonumber\\ \mathbf\nabla\mathbf v &= \frac{\partial \mathbf v}{\partial \mathbf x} = \dot{\mathbf F}{\mathbf F^{-1}}\label{eq:spatialgrad}. \end{align} $$ Отже, використовуючи \ref{eq:spatialgrad} маємо, що $\dot{\mathcal J} = \mathcal J \operatorname{trace}(\mathbf\nabla \mathbf v).$ Підставляючи в \ref{eq:rho0dert} отримуємо: $$ \frac{\partial \rho_0}{\partial t} = \left(\frac{\partial \rho}{\partial t} + \nabla \rho\cdot \mathbf v + \rho \operatorname{trace}(\mathbf\nabla\mathbf v)\right)\mathcal J = 0. $$ Тут ми знаємо, що $\mathcal J \ne 0,$ властивість непроникності речовини, отже ми можемо поділити на $\mathcal J$ і отримати: $$ \frac{\partial \rho}{\partial t} + \nabla \rho\cdot \mathbf v + \rho \operatorname{trace}(\mathbf\nabla\mathbf v) = 0. $$ Щоб спростити далі зауважимо, що $\mathbf\nabla\mathbf v = \frac{\partial \mathbf v_{i}}{\partial \mathbf x_j}\mathbf e_i \otimes \mathbf e_j,$ відповідно слід цієї матриці дорівнює $\nabla \cdot \mathbf v.$ І, нарешті, рівняння збереження маси в поточній конфігурації таке: $$ \begin{equation} \boxed{ \frac{\partial \rho}{\partial t} + \nabla \rho\cdot \mathbf v + \rho \nabla \cdot \mathbf v = 0. }\label{eq:massconscur} \end{equation} $$ У випадку нестисних плинів $\nabla \cdot \mathbf v = 0,$ тому $$ \begin{equation} \boxed{ \frac{\partial \rho}{\partial t} + \nabla \rho\cdot \mathbf v = 0. }\label{eq:massconscurincompress} \end{equation} $$ Зауважте, що ми можемо поєднати два останні доданки у \ref{eq:massconscur}: $$ \boxed{ \frac{\partial \rho}{\partial t} + \nabla \cdot( \rho \mathbf v) = 0. } $$

Див. також

вівторок, 30 травня 2017 р.

Оператор Лапласа: від декартових до сферичних координат

Оператор Лапласа дуже важливий, він повсюдно зустрічається у фізиці. Його форма в декартових координатах проста і симетрична: $$ \begin{equation} \nabla^2 = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2} \label{eq:cartesian}\end{equation} $$

Вираження у сферичних координатах:

Перехід до декартових від сферчних: $$ \begin{align} x &=r \sin{\theta}\cos{\phi},\label{eq:x}\\ y &=r \sin{\theta}\sin{\phi},\label{eq:y}\\ z &=r \cos{\theta}.\label{eq:z} \end{align} $$ І навпаки: $$ \begin{align} r &= \sqrt{x^2+y^2+z^2}, \label{eq:r}\\ \cos\phi &= \frac{x}{\sqrt{x^2+y^2}}, \label{eq:cosphi}\\ \sin\phi &= \frac{y}{\sqrt{x^2+y^2}}, \label{eq:sinphi}\\ \cos\theta &= \frac{z}{\sqrt{x^2+y^2+z^2}}, \label{eq:costheta}\\ \sin\theta &= \frac{\sqrt{x^2+y^2}}{\sqrt{x^2+y^2+z^2}}. \label{eq:sintheta} \end{align} $$ Використовуючи ланцюгове правило маємо: $$ \begin{align} \frac{\partial}{\partial x} = \frac{\partial}{\partial r}\frac{\partial r}{\partial x} + \frac{\partial}{\partial \phi}\frac{\partial \phi}{\partial x} + \frac{\partial}{\partial \theta}\frac{\partial \theta}{\partial x},\\ \frac{\partial}{\partial y} = \frac{\partial}{\partial r}\frac{\partial r}{\partial y} + \frac{\partial}{\partial \phi}\frac{\partial \phi}{\partial y} + \frac{\partial}{\partial \theta}\frac{\partial \theta}{\partial y},\\ \frac{\partial}{\partial z} = \frac{\partial}{\partial r}\frac{\partial r}{\partial z} + \frac{\partial}{\partial \phi}\frac{\partial \phi}{\partial z} + \frac{\partial}{\partial \theta}\frac{\partial \theta}{\partial z}. \end{align} $$ Знайти $\frac{\partial r}{\partial x}, \frac{\partial r}{\partial y}, \frac{\partial r}{\partial z}$ не дуже складно. Використовуємо $\ref{eq:x}, \ref{eq:y}, \ref{eq:z}:$ $$ \begin{align} \frac{\partial r}{\partial x}&=\frac{x}{\sqrt{x^2+y^2+z^2}}=\sin\theta\cos\phi,\\ \frac{\partial r}{\partial y}&=\frac{y}{\sqrt{x^2+y^2+z^2}}=\sin\theta\sin\phi,\\ \frac{\partial r}{\partial z}&=\frac{z}{\sqrt{x^2+y^2+z^2}}=\cos\theta.\\ \end{align} $$ Трошки більше проблем з $\frac{\partial \phi}{\partial x}, \frac{\partial \phi}{\partial y}, \frac{\partial \phi}{\partial z}.$ Нам буде потрібен той факт, що $\sqrt{x^2+y^2} = r\sin\theta.$ Застосовуючи неявне диференціювання, почергово вважаючи $y$ і $x$ сталими:
$$ \begin{align*} -\sin\phi\ d\phi &= \left(\frac{1}{\sqrt{x^2+y^2}}-\frac{x^2}{(x^2+y^2)^{3/2}}\right)dx = \frac{y^2}{(x^2+y^2)^{3/2}}dx = \frac{\sin^2\phi}{r\sin\theta}dx,\\ \cos\phi\ d\phi &= \left(\frac{1}{\sqrt{x^2+y^2}}-\frac{y^2}{(x^2+y^2)^{3/2}}\right)dy = \frac{x^2}{(x^2+y^2)^{3/2}}dy = \frac{\cos^2\phi}{r\sin\theta}dy, \end{align*} $$ отже $$ \begin{align} \frac{\partial \phi}{\partial x}&=-\frac{\sin\phi}{r\sin\theta},\\ \frac{\partial \phi}{\partial y}&=\frac{\cos\phi}{r\sin\theta},\\ \frac{\partial \phi}{\partial z}&=0. \end{align} $$ Останнє рівняння має місце бо $\phi$ не залежить від $z.$

Використовуючи почергово $\ref{eq:costheta}, \ref{eq:costheta}$ і $\ref{eq:sintheta}$, подібним чином обчислюємо $\frac{\partial \theta}{\partial x}, \frac{\partial \theta}{\partial y}, \frac{\partial \theta}{\partial z}:$ $$ \begin{align*} -\sin\theta\ d\theta &= -\frac{zx}{(x^2+y^2+z^2)^{3/2}} dx = -\frac{\cos\ \theta\sin\ \theta\cos\ \phi}{r}dx,\\ -\sin\theta\ d\theta &= -\frac{zy}{(x^2+y^2+z^2)^{3/2}} dy = -\frac{\cos\ \theta\sin\ \theta\sin\ \phi}{r}dy,\\ \cos\theta\ d\theta &= -\frac{z\sqrt{x^2+y^2}}{(x^2+y^2+z^2)^{3/2}} dz = \frac{\cos\ \theta\sin\ \theta}{r}dz, \end{align*} $$ отже $$ \begin{align} \frac{\partial \theta}{\partial x}&=\frac{\cos\ \theta\cos\ \phi}{r},\\ \frac{\partial \theta}{\partial y}&=\frac{\cos\ \theta\sin\ \phi}{r},\\ \frac{\partial \theta}{\partial z}&=-\frac{\sin\ \theta}{r}. \end{align} $$

Тепер нам потрібно обчислити $\frac{\partial^2}{\partial x^2}, \frac{\partial^2}{\partial y^2}$ і $\frac{\partial^2}{\partial z^2}.$ Почнемо з найпростішого: $$ \frac{\partial^2}{\partial z^2} = \left(\cos\ \theta\frac{\partial}{\partial r} - \frac{\sin\ \theta}{r}\frac{\partial}{\partial \theta}\right)\left(\cos\ \theta\frac{\partial}{\partial r} - \frac{\sin\ \theta}{r}\frac{\partial}{\partial \theta}\right). $$ Далі гірше: $$ \begin{align*} \frac{\partial^2}{\partial x^2} &= \left(\sin\ \theta\cos\ \phi\frac{\partial}{\partial r} - \frac{\sin\ \phi}{r\sin\ \theta}\frac{\partial}{\partial \phi} + \frac{\cos\ \theta\cos\ \phi}{r}\frac{\partial}{\partial \theta}\right)\\ &\quad \times\left(\sin\ \theta\cos\ \phi\frac{\partial}{\partial r} - \frac{\sin\ \phi}{r\sin\ \theta}\frac{\partial}{\partial \phi} + \frac{\cos\ \theta\cos\ \phi}{r}\frac{\partial}{\partial \theta}\right),\\ \frac{\partial^2}{\partial y^2} &= \left(\sin\ \theta\sin\ \phi\frac{\partial}{\partial r} + \frac{\cos\ \phi}{r\sin\ \theta}\frac{\partial}{\partial \phi} + \frac{\cos\ \theta\sin\ \phi}{r}\frac{\partial}{\partial \theta}\right)\\ &\quad \times\left(\sin\ \theta\sin\ \phi\frac{\partial}{\partial r} + \frac{\cos\ \phi}{r\sin\ \theta}\frac{\partial}{\partial \phi} + \frac{\cos\ \theta\sin\ \phi}{r}\frac{\partial}{\partial \theta}\right). \end{align*} $$ Якщо ми будемо просто розкривати дужки, то збожеволіємо. Тому ми будемо розглядати доданки зі спільним множником. Почнемо з $\frac{\partial^2}{\partial r^2}$: $$[\cos^2\ \theta + \sin^2\ \theta(\cos^2\ \phi + \sin^2\ \phi)] = 1.$$ Себто, перший доданок результату -- це $\boxed{\frac{\partial^2}{\partial r^2}}.$ Гляньмо на те, що відбувається з доданками де кожен множник містить $\frac{\partial}{\partial \phi}.$ Розпишемо ті з них, які в рівняннях із $\frac{\partial}{\partial x^2}$ і $\frac{\partial}{\partial y^2}$ відповідно: $$ \begin{align*} \left(-\frac{\sin\ \phi}{r\sin\ \theta}\frac{\partial}{\partial \phi}\right)\cdot\left(-\frac{\sin\ \phi}{r\sin\ \theta}\frac{\partial}{\partial \phi}\right) &=\frac{\sin\ \phi\cos\ \phi}{r^2\sin^2\ \theta}\frac{\partial}{\partial \phi}+\frac{\sin^2\ \phi}{r^2\sin^2\ \theta}\frac{\partial^2}{\partial \phi^2},\\ \left(-\frac{\cos\ \phi}{r\sin\ \theta}\frac{\partial}{\partial \phi}\right)\cdot\left(-\frac{\cos\ \phi}{r\sin\ \theta}\frac{\partial}{\partial \phi}\right) &=-\frac{\cos\ \phi\sin\ \phi}{r^2\sin^2\ \theta}\frac{\partial}{\partial \phi}+\frac{\cos^2\ \phi}{r^2\sin^2\ \theta}\frac{\partial^2}{\partial \phi^2}. \end{align*} $$ У сумі виходить $\boxed{\frac{1}{r^2\sin^2\ \theta}\frac{\partial^2}{\partial \phi^2}}.$ Переходимо до множення додантків із $\frac{\partial}{\partial \theta}:$ $$ \begin{align*} &\left(\frac{\sin\ \theta\cos\ \theta}{r^2} - \frac{\cos\ \theta\sin\ \theta\cos^2\ \phi}{r^2} - \frac{\cos\ \theta\sin\ \theta\sin^2\ \phi}{r^2}\right)\frac{\partial}{\partial \theta}\\ &+\left(\frac{\sin^2\ \theta}{r^2} + \frac{\cos^2\ \theta\cos^2\ \phi}{r^2} + \frac{\cos^2\ \theta\sin^2\ \phi}{r^2}\right)\frac{\partial^2}{\partial \theta^2}. \end{align*} $$ Отже, маємо $\boxed{\frac{1}{r^2}\frac{\partial^2}{\partial \theta^2}}.$

Попереду ще багатенько, хай йому грець! Добре, наступне коло, беремося за всі доданки з $\frac{\partial}{\partial r}$ і $\frac{\partial}{\partial \phi}:$ $$ \begin{align*} &\left(\frac{\sin\ \theta\cos\ \phi\sin\ \phi}{r^2\sin\ \phi} - \frac{\sin\ \theta\sin\ \phi\cos\ \phi}{r^2\sin\ \phi}\right)\frac{\partial}{\partial \phi}\\ &+\left(\frac{\sin^2\ \phi\sin\ \theta}{r\sin\ \theta} + \frac{\cos^2\ \phi\sin\ \theta}{r\sin\ \theta}\right)\frac{\partial}{\partial r}\\ &+\bigg(-\frac{\sin\ \theta\cos\ \phi\sin\ \phi}{r\sin\ \theta} - \frac{\sin\ \phi\sin\ \theta\cos\ \phi}{r\sin\ \theta}\\ &\quad\quad + \frac{\sin\ \theta\sin\ \phi\cos\ \phi}{r\sin\ \theta} + \frac{\cos\ \phi\sin\ \theta\sin\ \phi}{r\sin\ \theta}\bigg)\frac{\partial}{\partial r}\frac{\partial}{\partial \phi}\\ =&\boxed{\frac{1}{r}\frac{\partial}{\partial r}}. \end{align*} $$ Тепер черга за доданками з $\frac{\partial}{\partial r}$ і $\frac{\partial}{\partial \theta}:$ $$ \begin{align*} &\left(\frac{\cos\ \theta\sin\ \theta}{r^2}-\frac{\sin\ \theta\cos^2\ \phi\cos\ \theta}{r^2}-\frac{\sin\ \theta\sin^2\ \phi\cos\ \theta}{r^2}\right)\frac{\partial}{\partial \theta}\\ &+\left(\frac{\sin^2\ \theta}{r} + \frac{\cos^2\ \theta\cos^2\ \phi}{r} + \frac{\cos^2\ \theta\sin^2\ \phi}{r} \right)\frac{\partial}{\partial r}\\ &+\bigg(-\frac{\cos\ \theta\sin\ \theta}{r}-\frac{\sin\ \theta\cos\ \theta}{r}\\ &\quad\quad + \frac{\sin\ \theta\cos\ \theta\cos^2\ \phi}{r} + \frac{\sin\ \theta\cos\ \theta\sin^2\ \phi}{r}\\ &\quad\quad + \frac{\sin\ \theta\cos\ \theta\sin^2\ \phi}{r} + \frac{\sin\ \theta\cos\ \theta\sin^2\ \phi}{r}\bigg)\frac{\partial}{\partial r}\frac{\partial}{\partial \theta}\\ =&\boxed{\frac{1}{r}\frac{\partial}{\partial r}}. \end{align*} $$ І остання двійка множників -- це $\frac{\partial}{\partial \phi}$ і $\frac{\partial}{\partial \theta}:$ $$ \begin{align*} &\left(\frac{\sin^2\ \phi\cos\ \theta}{r^2\sin\ \theta}+\frac{\cos^2\ \phi\cos\ \theta}{r^2\sin\ \theta}\right)\frac{\partial}{\partial \theta}\\ &+\left(\frac{\cos^2\ \theta\cos\ \phi\sin\ \phi}{r\sin^2\ \theta} - \frac{\cos^2\ \theta\sin\ \phi\cos\ \phi}{r\sin^2\ \theta}\right)\frac{\partial}{\partial \phi}\\ &+\bigg(-\frac{\sin\ \phi\cos\ \theta\cos\ \phi}{r^2\sin\ \theta}-\frac{\cos\ \theta\cos\ \phi\sin\ \phi}{r^2\sin\ \theta}\\ &\quad\quad + \frac{\cos\ \phi\cos\ \theta\sin\ \phi}{r^2\sin\ \theta} + \frac{\sin\ \phi\cos\ \theta\cos\ \phi}{r^2\sin\ \theta}\bigg)\frac{\partial}{\partial \phi}\frac{\partial}{\partial \theta}\\ =&\boxed{\frac{\cos\ \theta}{r^2\sin\ \theta}\frac{\partial}{\partial \theta}}. \end{align*} $$ Тепер збираємо всі доданки обведені прямокутниками, щоб отримати формулу оператора Лапласа у сферичних координатах: $$\nabla^2 = \frac{\partial^2}{\partial r^2} + \frac{2}{r}\frac{\partial}{\partial r} + \frac{1}{r^2}\bigg[\frac{\partial^2}{\partial \theta^2} +\frac{\cos\ \theta}{\sin\ \theta}\frac{\partial}{\partial \theta} + \frac{1}{\sin^2\ \theta}\frac{\partial^2}{\partial \phi^2}\bigg].$$ Ми можемо переписати її в трошки охайнішій формі: $$\boxed{\nabla^2 = r^2\frac{\partial}{\partial r}\frac{1}{r^2}\frac{\partial}{\partial r} + \frac{1}{r^2\sin\ \theta}\frac{\partial}{\partial\theta}\sin\ \theta\frac{\partial}{\partial\theta} + \frac{1}{r^2\sin^2\ \theta}\frac{\partial^2}{\partial\phi^2}}.$$ Зауважте, що домноження на $r^2$ повністю розділяє кутову частину від радіальної. Ось чому уся робота була варта того.

пʼятницю, 26 травня 2017 р.

Лапласіан радіальної функції

Нехай $f:\mathbb{R}^n\to\mathbb{R} - $ це радіальна функція, тобто $f(x)=f(r)$ де $r:=\left\|x\right\|_2$. Тоді $$\Delta f=\frac{1}{r^{n-1}}\frac{d}{dr}(r^{n-1}f').$$

Доведення:

З одного боку, $$\frac{1}{r^{n-1}}\frac{d}{dr}\left(r^{n-1}f'\right) = \frac{1}{r^{n-1}} \left( (n-1)r^{n-2} f' + r^{n-1} f'' \right) = \frac{n-1}{r}f' + f''$$ З іншого боку, градієнт $f - $ це $$ \begin{align*} \displaystyle \nabla f& = \left(\frac{\partial}{\partial x_1}f,\dots,\frac{\partial}{\partial x_n}f\right)\\ & = \Bigg(\frac{\partial r}{\partial x_1}\frac{\partial f}{\partial r},\dots,\frac{\partial r}{\partial x_n}\frac{\partial f}{\partial r}\Bigg)\\ & = \Bigg(\frac{x_1}{r}f',\dots,\frac{x_n}{r}f'\Bigg)\\ & = \frac{\vec{r}}{r}f',\\ \end{align*} $$ де $\frac{\partial r}{\partial x_i} = \left(\sqrt{\sum_{j=1}^n x_j^2}\right)'_{x_i} = \frac{2x_i}{2\sqrt{\sum_{j=1}^n x_j^2}} = \frac{x_i}{r}.$

Тепер ми можемо обчислити Лапласіан $f.$

$$ \begin{align*} \nabla\cdot\nabla f& = \nabla\cdot\left(\frac{x_1}{r}f',\dots,\frac{x_n}{r}f'\right) = \sum_{i=1}^n \left(\frac{x_i}{r}f'\right)'_{x_i}\\ & = \sum_{i=1}^n \left(\frac{1\cdot r - x_i\cdot x_i/r}{r^2}f' + \frac{x_i}{r}f''\frac{\partial r}{\partial x_i}\right)\\ & = \left( \frac{n}{r} - \frac{\vec{r}\cdot\vec{r}}{r^3}\right) f' + \frac{\vec{r}\cdot\vec{r}}{r^2}f'' = \frac{n-1}{r}f' + f'' \end{align*} $$ Отже, $\boxed{\displaystyle \Delta f = \frac{1}{r^{n-1}}\frac{d}{dr}\left(r^{n-1}f'\right).}$

неділю, 5 лютого 2017 р.

По дорозі до Нав'є — Стокса

Загальна транспортна теорема

Нехай $\mathbf F$ буде гладким векторним або скалярним полем в області $\mathcal R(t)$ чиї межі це $\mathcal S(t),$ і нехай $\mathbf W$ буде полем швидкостей часозалежного руху $\mathcal S(t).$ Тоді $$\frac{d}{dt}\int_{\mathcal R(t)}\mathbf F(x, t)\ dV = \int_{\mathcal R(t)}\frac{\partial \mathbf F}{\partial t}\ dV + \int_{\mathcal S(t)}\mathbf F\mathbf W \cdot \mathbf n\ dA.$$

Транспортна теорема Рейнольдса

Нехай $\mathbf\Phi$ буде гладким векторним або скалярним полем, і припустимо, що $\mathcal R(t)$ це пакунок плину з поверхнею $\mathcal S(t),$ що мандрує зі швидкістю потоку $\mathbf U.$ Тоді $$\frac{D}{Dt}\int_{\mathcal R(t)}\mathbf\Phi\ dV =\int_{\mathcal R(t)}\frac{\partial\mathbf\Phi}{\partial t}\ dV +\int_{\mathcal S(t)}\mathbf\Phi \mathbf U \cdot \mathbf n\ dA.$$

Виведення з Загальної транспортної теореми.

Якщо наша область рухається зі швидкістю потоку, то $\mathbf W = \mathbf U$ і загальна похідна перетворюється на матеріальну.1

Рівняння неперервності для контрольного об'єму

Почнемо з розглядання фіксованої маси плину $m,$ що міститься в довільній області $\mathcal R(t).$ $$m=\int_{\mathcal R(t)}\rho\ dV.$$ $\mathcal R(t), m$ можуть змінюватись з часом, але маса повинна залишатись незмінною. Наприклад, можна розглянути повітряну кульку з теплим повітрям у холодному середовищі. $$\frac{dm}{dt}=\frac{d}{dt}\int_{\mathcal R(t)}\rho\ dV.$$ Застосувавши загальну транспортну теорему маємо: $$\int_{\mathcal R(t)}\frac{\partial \rho}{\partial t}\ dV + \int_{\mathcal S(t)}\rho \mathbf W\cdot \mathbf n\ dA = 0.$$ У плинових системах зазвичай зручно вважати, що поле швидкостей $\mathbf W$ було полем швидкостей плину, що відповідає розгляданню на локальному рівні нашої області як пакунок плину. $$\int_{\mathcal R(t)}\frac{\partial \rho}{\partial t}\ dV + \int_{\mathcal S(t)}\rho \mathbf U\cdot \mathbf n\ dA = 0.$$ Використавши формула Остроградського-Гауса ми позбуваємось двох типів інтегралів залишаючи лише інтеграл по об'єму: $$\int_{\mathcal R(t)}\frac{\partial \rho}{\partial t} + \nabla \cdot \rho \mathbf U\ dV = 0.$$ З того, що область $\mathcal R(t)$ було обрана довільно ми можемо вважати, що це діє і для наскільки завгодно малої області (в рамках нашої гіпотези неперервності), тобто інтеграл мусить бути нулем усюди в $\mathcal R(t)$. Якщо це не так ми можемо розділити $\mathcal R(t)$ на дві підобласті додатну і від'ємну і показити цим порушення рівності нулю. Отже, $$\begin{equation}\frac{\partial \rho}{\partial t} + \nabla \cdot \rho \mathbf U=0\label{eq:cont_eq_df}\end{equation}.$$ Це диференціальна форма рівняння неперервності.
Якщо розписати оператор дивергенції, то маємо $$\begin{equation} \rho_t + (\rho u)_t + (\rho v)_t + (\rho w)_t = 0\label{eq:cont_eq_div_decomp}\end{equation}.$$ Тепер можна використати добутку для похідних $$\frac{D\rho}{Dt} + \rho\nabla\cdot \mathbf U = 0\label.$$ У випадку нестисного потоку, $\rho - $ стале, маємо $$\begin{equation}u_x + v_y + w_z = 0\label{eq:incompres_uvw}\end{equation}.$$

Аналіз у випадку контрольного об'єму

Оскільки $\ref{eq:cont_eq_df}$ виконується для будь-якого об'єму, то ми можемо записати його і для пакунку плину $$\int_{\mathcal R(t)} \frac{\partial \rho}{\partial t} + \nabla \cdot \rho \mathbf U\ dV = 0.$$ Тепер ми використаємо формулу Остроградського-Гауса в іншому напрямку $$\int_{\mathcal R(t)} \frac{\partial \rho}{\partial t}\ dV + \int_{\mathcal S(t)} \rho \mathbf U \cdot \mathbf n = 0.$$ Тепер Загальну транспортну теорему $$\frac{d}{dt} \int_{\mathcal R(t)}\rho\ dV + \int_{\mathcal S(t)}\rho(\mathbf U - \mathbf W)\cdot \mathbf n \ dA = 0.$$ Якщо $\mathcal S_e(t)$ - зона входів і виходів через які плин може втікати і витікати, то $$\frac{d}{dt} \int_{\mathcal R(t)}\rho\ dV + \int_{\mathcal S_e(t)}\rho(\mathbf U - \mathbf W)\cdot \mathbf n \ dA = 0.$$ Це рівняння називається рівняння неперервності для контрольного об'єму і виражає очевидний факт, що
$\{$ швидкість збільшення маси
контрольного об'єму
$\}$ = $\{$ сумарний потік
крізь поверхню
$\}$

Спрощення

Якщо густина незмінна з часом, а контрольний об'єм нерухомий, тобто $\mathbf W = 0,$ тоді можна записати $$\begin{equation}\int_{\mathcal S_e(t)}\rho\mathbf U \cdot \mathbf n\ dA = 0\label{eq:cv_simple}\end{equation}.$$ Рівняння $\ref{eq:cv_simple}$ вимірюється в одиницях $\rho U A,$ тобто $\frac{M}{L^3}\frac{L}{T}L^2 = \frac{M}{T}.$ Ми називаємо це масовою ви́тратою $\dot m = \rho U A.$ У випадку часозалежної поверхні маємо $$\dot m = \int_{\mathcal S_e(t)} \rho(\mathbf U - \mathbf W)\cdot n\ dA.$$

Рівняння імпульсу потоку плину

Тут нам знадобиться закон Ньютона, який ми запишемо в альтернативній формі $$F/V = \frac{d}{dt}(\rho u).$$ На відміну від випадку з рівнянням неперевності наша область одразу буде пакунком плину. Якщо ми використовуємо Ейлерову точку зору, то швидкість зміни вектора імпульсу $\frac{D}{Dt}\int_{\mathcal R(t)}\rho \mathbf U\ dV.$ Маємо основні два типи сил, які треба врахувати:
  • масова сила - $\int_{\mathcal R(t)}\mathbf F_B\ dV,$
  • поверхнева сила - $\int_{\mathcal S(t)}\mathbf F_S\ dA.$
$$\begin{equation}\frac{D}{Dt}\int_{\mathcal R(t)}\rho \mathbf U\ dV = \int_{\mathcal R(t)}\mathbf F_B\ dV + \int_{\mathcal S(t)}\mathbf F_S\ dA\label{eq:momentum}\end{equation}.$$ Не дуже зручно мати в одному рівнянні два типи інтегралів: по поверхні і по об'єму. Також хотілось би використовувати змінні, які звичні інженерних обчислень, тому потрібно виразити об'єм і поверхню через ці змінні. Зауважимо, що масові сили досить прості, тобто зазвичай це лише виштовхувальна сила через гравітаційне прискорення, $\mathbf F_B = \rho \mathbf g.$ Натомість робота з поверхневими силами вимагає певних серйозних зусиль, тому спочатку ми спростимо рівняння внесенням диференціювання під знак інтеграла. Для цього ми використаємо транспортну теорему Рейнольдса. $$\frac{D}{Dt}\int_{\mathcal R(t)}\rho u\ dV =\int_{\mathcal R(t)}\frac{\partial \rho u}{\partial t}\ dV +\int_{\mathcal S(t)}\rho u \mathbf U \cdot \mathbf n\ dA.$$ Знов використаємо формулу Остроградського-Гауса, щоб отримати інтеграл по об'єму праворуч, під знаком інтегралу маємо $$\frac{\partial \rho u}{\partial t} + \nabla \cdot (\rho u \mathbf U) = \frac{\partial \rho u}{\partial t} + \mathbf U \cdot \nabla (\rho u) + \rho u \nabla \cdot \mathbf U.$$ Тому що ми маємо справу з нестисним потоком, ми можемо використати $\ref{eq:incompres_uvw}$ $$\begin{equation}\frac{D}{Dt}\int_{\mathcal R(t)}\rho u\ dV =\int_{\mathcal R(t)}\rho\frac{\partial u}{\partial t} + \rho\mathbf U \cdot \nabla u\ dV =\int_{\mathcal R(t)}\rho\frac{D u}{D t}\ dV \label{eq:for_momentum_under_int}\end{equation},$$ де $u$ - це $x-$компонент швидкості. Перезапишемо $\ref{eq:momentum}$ використавши $\ref{eq:for_momentum_under_int}$: $$\begin{equation}\int_{\mathcal R(t)}\rho \frac{D\mathbf U}{Dt}\ dV = \int_{\mathcal R(t)}\mathbf F_B\ dV + \int_{\mathcal S(t)}\mathbf F_S\ dA\label{eq:momentum_under_int}\end{equation}.$$

Ну що ж, нам залишились два кроки: перевести рівняння $\ref{eq:momentum_under_int}$ до диференціальної форми і розібратись із поверхневими силами. Зауважимо, що $\mathbf F_S$ - це вектор і це говорить, що має існувати матриця $\Lambda,$ така що $\mathbf F_S = \mathbf{\Lambda} \cdot \mathbf n.$ З того, що $\mathbf n$ - це чисто геометрична величина, нормаль до поверхні, то фізична сутність поверхневих сил мусить зберігатись в матриці $\mathbf{\Lambda}.$

Тепер ми готові записати базову форму диференціального рівняння імпульсу: $$\int_{\mathcal R(t)}\rho \frac{D\mathbf U}{Dt}\ dV = \int_{\mathcal R(t)}\mathbf F_B\ dV + \int_{\mathcal S(t)}\mathbf{\Lambda}\cdot\mathbf n\ dA,$$ Тому що $\mathcal R(t)$ - це довільний елемент плину і, отже, ми можемо обрати його як завгодно малим, то підінтегральний вираз також дорівнює нулю: $$\begin{equation}\rho \frac{D\mathbf U}{Dt} - \mathbf F_B - \nabla\cdot\mathbf{\Lambda} = 0\label{eq:momentum_dif_form}\end{equation}.$$ Це дуже загальний баланс імпульсу, що виконується для всіх точок будь-якого потоку плину (в межах гіпотези неперервності).

Трактування поверхневих сил

Поверхневе напруження складніша величина. Причиною цьому є те, що ми не можемо говорити про напруження в точці не визначивши поверхню на яку це напруження діє. $$\Lambda = \begin{bmatrix} -p + \tau_{xx} & \tau_{xy} & \tau_{xz} \\ \tau_{yx} & -p + \tau_{yy} & \tau_{yz} \\ \tau_{zx} & \tau_{zy} & -p + \tau_{zz} \\ \end{bmatrix} = -p\mathbf I + \mathbf\tau.$$ $\mathbf\tau$ відомий як в'язкісний тензор напружень.

Рівняння Нав'є-Стокса

$$ \begin{align*} \nabla\cdot\Lambda &= \big(-\frac{\partial p}{\partial x} + \frac{\partial\tau_{xx}}{\partial x} + \frac{\partial\tau_{yx}}{\partial y} + \frac{\partial\tau_{zx}}{\partial z}\big)\mathbf e_1 +\\ &= \big(-\frac{\partial p}{\partial y} + \frac{\partial\tau_{xy}}{\partial x} + \frac{\partial\tau_{yy}}{\partial y} + \frac{\partial\tau_{zy}}{\partial z}\big)\mathbf e_2 +\\ &= \big(-\frac{\partial p}{\partial z} + \frac{\partial\tau_{xz}}{\partial x} + \frac{\partial\tau_{yz}}{\partial y} + \frac{\partial\tau_{zz}}{\partial z}\big)\mathbf e_3. \end{align*} $$ Тепер ми можемо переписати $\ref{eq:momentum_dif_form}$ у вигляді $$\begin{equation}\rho\frac{Du}{Dt} = -\frac{\partial p}{\partial x} + \frac{\partial\tau_{xx}}{\partial x} + \frac{\partial\tau_{yx}}{\partial y} + \frac{\partial\tau_{zx}}{\partial z} + F_{B,x}\label{eq:momentum_df_tau}\end{equation}$$ Згадаємо багатовимірний закон в'язкості Ньютона, який надає такі формули для компонентів напруження в'язкості: $$\tau_{xx} = 2\mu u_x, \tau_yx = \mu(u_y+v_x), \tau_zx = \mu=\mu(u_z+w_x).$$ Підставляючи це у другий, третій і четвертий доданки $\ref{eq:momentum_df_tau}$ маємо $$2\mu u_{xx}+\mu(u_y+v_x)_y+\mu(u_z+w_x)_z.$$ Ми можемо переформатувати це так: $$\mu (u_{xx} + u_{yy} + u_{zz}) + \mu (u_x + v_y + w_z)_x,$$ але ми маємо справу з нестисним потоком, тому згідно з $\ref{eq:incompres_uvw}$ дивергенція вектора швидкості дорівнює нулю і залишається лише перший доданок. $$\rho\frac{Du}{Dt}=-p_x+\mu(u_{xx}+u_{yy}+u_{zz})+F_{B,x}.$$ Зазвичай це рівняння ділять на $\rho$ і представляти його як $$u_t + uu_x+vu_y+wu_z = -\frac{1}{\rho}p_x+\nu\Delta u+\frac{1}{\rho}F_{B,x}.$$ Тут, $\nu$ - це кінематична в'язкість.

Рівняння для всіх трьох координат називаються рівняннями Нав'є-Стокса, і вони надають опис в кожній точці для по суті будь-якого нестисного потоку.

Аналіз рівнянь Нав'є-Стокса

$$\overbrace{\underbrace{u_t}_{\substack{локальне \\ прискорення}} + \underbrace{uu_x+vu_y+wu_z}_{\substack{конвекційне\\ прискорення}}}^{повне\ прискорення} = \overbrace{-\frac{1}{\rho}p_x}^{сили\ тиску}+ \overbrace{\nu\Delta u}^{сили\ в'язкості}+ \overbrace{\frac{1}{\rho}F_{B,x}}^{масові\ сили}.$$ Доданки ліворуч в контексті рівнянь Н.-С. часто називають інерційними.

Сили в'язкості

$$\nu\Delta u = \nu(u_{xx}+u_{yy}+u_{zz}).$$ Загалом, доданки у вигляді другої похідної у диференційному рівнянні зазвичай пов'язують із дифузією, і в в обох фізичному і математичному контекстах представляють розмазування, згладжування або змішування.

Ех, зібратись би колись і закінчити це

Примітки