arma-thesis

git clone https://git.igankevich.com/arma-thesis.git
Log | Files | Refs | LICENSE
commit ba4f0f8b14352847b7541258949d3f78d622f959
parent 96b7528afc53b9c649f0d585fa3a9f8a1cb53616
Author: Ivan Gankevich <igankevich@ya.ru>
Date:   Mon, 13 Nov 2017 15:26:52 +0300

Reorder sections and parts.

Diffstat:

arma-thesis-ru.org | 1904+++++++++++++++++++++++++++++++++++++++----------------------------------------


arma-thesis.org | 1735+++++++++++++++++++++++++++++++++++++++----------------------------------------


2 files changed, 1781 insertions(+), 1858 deletions(-)

diff --git a/arma-thesis-ru.org b/arma-thesis-ru.org
@@ -633,827 +633,933 @@ Motion Programme (LAMP), программе для моделирования к
 при условии наличия соответствующих формул корректировки коэффициентов, что
 является целью дальнейших исследований.
 
-** Известные формулы определения поля давлений
-**** Теория волн малых амплитуд.
-В\nbsp{}cite:stab2012,degtyarev1998modelling,degtyarev1997analysis дается
-решение обратной задачи гидродинамики для случая идеальной несжимаемой жидкости
-в рамках теории волн малых амплитуд (в предположении, что длина волны много
-больше ее высоты: \(\lambda \gg h\)). В этом случае обратная задача линейна и
-сводится к уравнению Лапласа со смешанным граничным условием, а уравнение
-движения используется только для нахождения давлений по известным значениям
-производных потенциала скорости. Предположение о малости амплитуд волн означает
-слабое изменение локального волнового числа во времени и пространстве по
-сравнению с подъемом (аппликатой) взволнованной поверхности. Это позволяет
-вычислить производную подъема поверхности по \(z\) как \(\zeta_z=k\zeta\), где
-\(k\)\nbsp{}--- волновое число. В двухмерном случае решение записывается явной
-формулой
-\begin{align}
-    \left.\frac{\partial\phi}{\partial x}\right|_{x,t}= &
-        -\frac{1}{\sqrt{1+\alpha^{2}}}e^{-I(x)}
-            \int\limits_{0}^x\frac{\partial\dot{\zeta}/\partial
-                z+\alpha\dot{\alpha}}{\sqrt{1+\alpha^{2}}}e^{I(x)}dx,\label{eq-old-sol-2d}\\
-    I(x)= & \int\limits_{0}^x\frac{\partial\alpha/\partial z}{1+\alpha^{2}}dx,\nonumber
-\end{align}
-где \(\alpha\)\nbsp{}--- уклоны волн. В трехмерном случае решение записывается в
-виде эллиптического дифференциального уравнения в частных производных (ДУЧП):
-\begin{align*}
-    & \frac{\partial^2 \phi}{\partial x^2} \left( 1 + \alpha_x^2 \right) +
-    \frac{\partial^2 \phi}{\partial y^2} \left( 1 + \alpha_y^2 \right) +
-    2\alpha_x\alpha_y \frac{\partial^2 \phi}{\partial x \partial y} + \\
-    & \left(
-        \frac{\partial \alpha_x}{\partial z} +
-        \alpha_x \frac{\partial \alpha_x}{\partial x} +
-        \alpha_y \frac{\partial \alpha_x}{\partial y}
-    \right) \frac{\partial \phi}{\partial x} + \\
-    & \left(
-        \frac{\partial \alpha_y}{\partial z} +
-        \alpha_x \frac{\partial \alpha_y}{\partial x} +
-        \alpha_y \frac{\partial \alpha_y}{\partial y}
-    \right) \frac{\partial \phi}{\partial y} + \\
-    & \frac{\partial \dot{\zeta}}{\partial z} +
-    \alpha_x \dot{\alpha_x} + \alpha_y \dot{\alpha_y} = 0.
-\end{align*}
-Уравнение предполагается решать численно путем сведения к разностному.
-
-Как будет показано в разделе\nbsp{}[[#sec-compare-formulae]]
-формула\nbsp{}eqref:eq-old-sol-2d расходится при попытке вычислить поле
-скоростей для волн больших амплитуд, а значит не может быть использована
-совместно с моделью морского волнения, генерирующей волны произвольных амплитуд.
-
-**** Линеаризация граничного условия.
-:PROPERTIES:
-:CUSTOM_ID: linearisation
-:END:
-
-Модель ЛХ позволяет вывести явную формулу для поля скоростей путем линеаризации
-кинематического граничного условия. Формула для потенциала скорости запишется
-как
-\begin{equation*}
-\phi(x,y,z,t) = \sum_n \frac{c_n g}{\omega_n}
-     e^{\sqrt{u_n^2+v_n^2} z}
-     \sin(u_n x + v_n y - \omega_n t + \epsilon_n).
-\end{equation*}
-Формула дифференцируется для получения производных потенциала, которые
-подставляются в динамическое граничное условие для определения поля давлений.
-
-** Определение поля давлений под дискретно заданной взволнованной поверхностью
-Аналитические решения граничных задач для классических уравнений часто
-используются для исследования различных свойств уравнений, и для таких
-исследований запись формулы общего решения неудобна ввиду своей сложности и
-наличия интегралов от неизвестных функций. Одним из методов нахождения
-аналитических решений ДУЧП является метод Фурье. Основой метода служит
-преобразование Фурье, применение которого к некоторым ДУЧП позволяет свести их к
-алгебраическому, а его решение записывается как обратное преобразование Фурье от
-некоторой функции (которая может содержать преобразования Фурье от других
-функций). Поскольку эти преобразования не всегда можно записать аналитически, то
-вместо этого ищутся частные решения задачи и анализируется их поведение в
-различных областях. В то же время, численный расчет дискретных преобразований
-Фурье возможен для любой дискретно заданной функции, используя алгоритмы БПФ.
-Эти алгоритмы используют симметрию комплексных экспонент для понижения
-асимптотической сложности с \(\mathcal{O}(n^2)\) до \(\mathcal{O}(n\log_{2}n)\).
-Таким образом, даже если общее решение содержит преобразования Фурье от
-неизвестных функций, они все равно могут быть взяты численно, а использование
-алгоритмов БПФ делает этот подход эффективным.
-
-Альтернативным подходом к решению ДУЧП является их сведение к разностным
-уравнениям, решаемым с помощью построения различных численных схем. При этом
-решение получается приближенным, а асимптотическая сложность соответствующих
-алгоритмов сопоставима со сложностью алгоритма БПФ. Например, для стационарного
-эллиптического уравнения в частных производных строится неявная численная схема;
-эта схема расчитывается итерационным методом, на каждом шаге которого ищется
-решение трехдиагональной или пятидиагональной СЛАУ методом прогонки.
-Асимптотическая сложность алгоритма составляет \(\mathcal{O}({n}{m})\), где
-\(n\)\nbsp{}--- количество точек на сетке взволнованной поверхности, а
-\(m\)\nbsp{}--- число итераций.
-
-Несмотря на широкое распространение, итеративные алгоритмы неэффективно
-отображаются на архитектуру параллельных машин ввиду неизбежной синхронизации
-процессов после каждой итерации; в частности, отображение на сопроцессоры может
-включать в себя копирование данных на сопроцессор и обратно на каждой итерации,
-что отрицательно сказывается на их производительности. В то же время, наличие
-большого количества преобразований Фурье в решении является скорее
-преимуществом, чем недостатком. Во-первых, решения, полученные с помощью метода
-Фурье, явные, а значит хорошо масштабируются на большое количество параллельно
-работающих вычислительных ядер с использованием простейших приемов параллельного
-программирования. Во-вторых, для алгоритмов БПФ существуют готовые
-оптимизированные реализации для различных архитектур процессоров и сопроцессоров
-(GPU, MIC). Эти преимущества обусловливают использование метода Фурье для
-получения явного решения задачи определения давлений под взволнованной морской
-поверхностью.
-
-*** Двухмерное поле скоростей
-:PROPERTIES:
-:CUSTOM_ID: sec-pressure-2d
-:END:
+** Моделирование нелинейности морских волн
+Модель АРСС позволяет учесть асимметричность распределения волновых аппликат,
+т.е.\nbsp{}генерировать морские волны, закон распределения аппликат которых
+имеет ненулевой эксцесс и асимметрию. Такой закон распределения характерен для
+реальных морских волн\nbsp{}cite:longuet1963nonlinear и задается либо
+полиномиальной аппроксимацией натурных данных, либо аналитически.
+Асимметричность волн моделируется с помощью нелинейного безынерционного
+преобразования (НБП) случайного процесса, однако, любое нелинейное
+преобразование случайного процесса приводит к преобразованию его АКФ. Для того
+чтобы подавить этот эффект, необходимо предварительно преобразовать АКФ, как
+показано в\nbsp{}cite:boukhanovsky1997thesis.
 
-**** Формула для жидкости бесконечной глубины.
-Задача Робена для уравнения Лапласа в двух измерениях записывается как
-\begin{align}
-    \label{eq-problem-2d}
-    & \phi_{xx}+\phi_{zz}=0,\\
-    & \zeta_t + \zeta_x\phi_x = \frac{\zeta_x}{\sqrt{1 + \zeta_x^2}} \phi_x - \phi_z, & \text{на }z=\zeta(x,t).\nonumber
-\end{align}
-Для ее решения воспользуемся методом Фурье. Возьмем преобразование Фурье от
-обоих частей уравнений Лапласа и получим
-\begin{equation*}
-    -4 \pi^2 \left( u^2 + v^2 \right)
-    \FourierY{\phi(x,z)}{u,v} = 0,
-\end{equation*}
-откуда имеем \(v = \pm i u\). Здесь и далее будет использоваться следующая
-симметричная форма преобразования Фурье:
-\begin{equation*}
-    \FourierY{f(x,y)}{u,v} =
-    \iint\limits_{-\infty}^{\phantom{--}\infty}
-    f(x,y)
-    e^{-2\pi i (x u + y v)}
-    dx dy.
-\end{equation*}
-Решение уравнения будем искать в виде обратного преобразования Фурье
-\(\phi(x,z)=\InverseFourierY{E(u,v)}{x,z}\). Подставляя[fn::Выражение \(v={-i}{u}\)
-не подходит в данной задаче, поскольку потенциал скорости должен стремиться к
-нулю с увеличением глубины до бесконечности.] \(v={i}{u}\) в формулу, решение
-перепишется как
+**** Преобразование взволнованной поверхности.
+Формула \(z=f(y)\) преобразования взволнованной поверхности к необходимому
+одномерному закону распределения \(F(z)\) получается путем решения нелинейного
+трансцендентного уравнения \(F(z) = \Phi(y)\), где \(\Phi(y)\)\nbsp{}--- функция
+одномерного нормального закона распределения. Поскольку функция распределения
+аппликат морских волн часто задается некоторой аппроксимацией, основанной на
+натурных данных, то это уравнение целесообразно решать численно в каждой точке
+\(y_k|_{k=0}^N\) сетки сгенерированной поверхности относительно \(z_k\). Тогда
+уравнение запишется в виде
 \begin{equation}
-    \label{eq-guessed-sol-2d}
-    \phi(x,z) = \InverseFourierY{e^{2\pi u z}E(u)}{x}.
+    \label{eq-distribution-transformation}
+    F(z_k)
+    =
+    \frac{1}{\sqrt{2\pi}}
+    \int\limits_0^{y_k} \exp\left[ -\frac{t^2}{2} \right] dt
+    .
 \end{equation}
-Для того чтобы подстановка \(z=\zeta(x,t)\) не помешала использованию
-преобразований Фурье в решении, перепишем\nbsp{}eqref:eq-guessed-sol-2d в виде
-свертки:
+Поскольку функции распределения монотонны, для решения этого уравнения
+используется простейший численный метод половинного деления (метод бисекции).
+
+**** Предварительное преобразование АКФ.
+Для преобразования АКФ \(\gamma_z\) процесса ее необходимо разложить в ряд по
+полиномам Эрмита (ряд Грама---Шарлье)
 \begin{equation*}
-    \phi(x,z)
+    \gamma_z \left( \vec u \right)
     =
-    \Fun{z}
-    \ast
-    \InverseFourierY{E(u)}{x},
+    \sum\limits_{m=0}^{\infty}
+    C_m^2 \frac{\gamma_y^m \left( \vec u \right)}{m!},
 \end{equation*}
-где \(\Fun{z}\)\nbsp{}--- некоторая функция, вид которой будет определен в
-разделе\nbsp{}[[#sec-compute-delta]] и для которой выполняется соотношение
-\(\FourierY{\Fun{z}}{u}=e^{2\pi{u}{z}}\). Подставляя выражение для \(\phi\) в
-граничное условие, получим
+где
 \begin{equation*}
-    \zeta_t
-    =
-    \left( i f(x) - 1 \right)
-    \left[
-        \Fun{z}
-        \ast
-        \InverseFourierY{2\pi u E(u)}{x}
-    \right],
+    C_m = \frac{1}{\sqrt{2\pi}}
+  \int\limits_{0}^\infty
+    f(y) H_m(y) \exp\left[ -\frac{y^2}{2} \right],
 \end{equation*}
-где \(f(x) = {\zeta_x}/{\sqrt{1 + \zeta_x^2}} - \zeta_x\). Применяя преобразование
-Фурье к обеим частям, получаем выражение для коэффициентов \(E\):
+\(H_m\)\nbsp{}--- полином Эрмита, а \(f(y)\)\nbsp{}--- решение
+уравнения\nbsp{}eqref:eq-distribution-transformation. Воспользовавшись
+полиномиальной аппроксимацией \(f(y) \approx \sum\limits_i d_i y^i\) и
+аналитическими выражениями для полиномов Эрмита, формулу определения
+коэффициентов можно упростить, используя следующее равенство:
 \begin{equation*}
-    E(u) =
-    \frac{1}{2\pi u}
-    \frac{
-    \FourierY{\zeta_t / \left(i f(x) - 1\right)}{u}
-    }{
-    \FourierY{\Fun{z}}{u}
-    }
+    \frac{1}{\sqrt{2\pi}}
+    \int\limits_\infty^\infty
+    y^k \exp\left[ -\frac{y^2}{2} \right]
+    =
+    \begin{cases}
+        (k-1)!! & \text{для четных }k,\\
+        0       & \text{для нечетных }k.
+    \end{cases}
 \end{equation*}
-Выполняя подстановку \(z=\zeta(x,t)\) и подставляя полученное выражение
-в\nbsp{}eqref:eq-guessed-sol-2d, получаем окончательное выражение для
-\(\phi(x,z)\):
+Оптимальное количество коэффициентов \(C_m\) определяется путем вычисления их
+последовательно и критерий прекращения счета определяется совпадением дисперсий
+обоих полей с требуемой точностью \(\epsilon\):
 \begin{equation}
-    \label{eq-solution-2d}
-    \boxed{
-        \phi(x,z)
-        =
-        \InverseFourierY{
-            \frac{e^{2\pi u z}}{2\pi u}
-            \frac{
-            \FourierY{ \zeta_t / \left(i f(x) - 1\right) }{u}
-            }{
-            \FourierY{ \Fun{\zeta(x,t)} }{u}
-            }
-        }{x}.
-    }
+    \label{eq-nit-error}
+    \left| \Var{z} - \sum\limits_{k=0}^m
+    \frac{C_k^2}{k!} \right| \leq \epsilon.
 \end{equation}
 
-Множитель \(e^{2\pi u z}/(2\pi u)\) делает график функции, от которой берется
-обратное преобразования Фурье, несимметричным относительно оси \(OY\). Это
-затрудняет применение БПФ, поскольку оно требует периодичную функцию,
-принимающую нулевые значения на концах промежутка. В то же время, использование
-численного интегрирования вместо БПФ не позволит получить преимущество над
-решением всей системы уравнений с помощью разностных схем. Эту проблему можно
-обойти, используя формулу\nbsp{}eqref:eq-solution-2d-full для жидкости конечной
-глубины с заведомо большим значением глубины водоема \(h\). Вывод формулы дан в
-следующем разделе.
+В\nbsp{}cite:boukhanovsky1997thesis автор предлагает использовать полиномиальную
+аппроксимацию для \(f(y)\) также для преобразования поверхности, однако на
+практике в реализации взволнованной поверхности часто находятся точки,
+выпадающие за промежуток на котором построена аппроксимация, что приводит к
+резкому уменьшению ее точности. В этих точках
+уравнение\nbsp{}eqref:eq-distribution-transformation эффективнее решать методом
+бисекции. Использование полиномиальной аппроксимации в формулах для
+коэффициентов ряда Грама---Шарлье не приводит к аналогичным ошибкам.
 
-**** Формула для жидкости конечной глубины.
-На дне водоема вертикальная составляющая скорости перемещения жидкости должна
-равняться нулю, т.е.\nbsp{}\(\phi_z=0\) на \(z=-h\), где \(h\)\nbsp{}--- глубина
-водоема. В этом случае пренебречь равенством \(v = -i u\), полученным из
-уравнения Лапласа, нельзя, и решение ищется в виде
-\begin{equation}
-    \phi(x,z)
-    =
-    \InverseFourierY{
-        \left( C_1 e^{2\pi u z} + C_2 e^{-2\pi u z} \right)
-        E(u)
-    }{x}.
-    \label{eq-guessed-sol-2d-full}
-\end{equation}
-Подставляя \(\phi\) в условие на дне водоема, получим
-\begin{equation*}
-    C_1 e^{-2\pi u h} - C_2 e^{2\pi u h} = 0,
-\end{equation*}
-откуда имеем \(C_1=\frac{1}{2}C{e}^{2\pi{u}{h}}\) и
-\(C_2=-\frac{1}{2}C{e}^{-2\pi{u}{h}}\). Константа \(C\) здесь произвольна,
-поскольку при подстановке станет частью неизвестных коэффициентов \(E(u)\).
-Подставляя полученные выражения для \(C_1\) и \(C_2\)
-в\nbsp{}eqref:eq-guessed-sol-2d-full, получаем выражение
-\begin{equation*}
-    \phi(x,z) = \InverseFourierY{ \Sinh{2\pi u (z+h)} E(u) }{x}.
-\end{equation*}
-Подставляя \(\phi\) в граничное условие на свободной поверхности, получаем
-\begin{equation*}
-    \zeta_t = f(x) \InverseFourierY{ 2\pi i u \Sinh{2\pi u (z+h)} E(u) }{x}
-            - \InverseFourierY{ 2\pi u \SinhX{2\pi u (z+h)} E(u) }{x}.
-\end{equation*}
-Здесь \(\sinh\) и \(\cosh\) дают схожие результаты вблизи свободной поверхности, и,
-поскольку эта область является наиболее интересной с точки зрения практического
-применения, положим \(\Sinh{2\pi{u}(z+h)}\approx\SinhX{2\pi{u}(z+h)}\). Выполняя
-аналогичные предыдущему разделу операции, получаем окончательное выражение для
-\(\phi(x,z)\):
+**** Разложение в ряд Грама---Шарлье.
+В\nbsp{}cite:huang1980experimental было экспериментально показано, что
+распределение аппликат морской поверхности отличается от нормального ненулевым
+эксцессом и асимметрией. В\nbsp{}cite:rozhkov1996theory показано, что такое
+распределение раскладывается в ряд Грама---Шарлье (РГШ):
+\begin{align}
+    \label{eq-skew-normal-1}
+    & F(z; \mu=0, \sigma=1, \gamma_1, \gamma_2) \approx
+    \Phi(z; \mu, \sigma) \frac{2 + \gamma_2}{2}
+    - \frac{2}{3} \phi(z; \mu, \sigma)
+    \left(\gamma_2 z^3+\gamma_1
+    \left(2 z^2+1\right)\right) \nonumber
+    \\
+    & f(z; \gamma_1, \gamma_2) \approx
+    \frac{1}{\sigma\sqrt{2 \pi}}e^{-\frac{(z-\mu)^2}{2\sigma^2}}
+    \left[
+        1+
+        \frac{1}{6} \gamma_1 z \left(z^2-3\right)
+        + \frac{1}{24} \gamma_2 \left(z^4-6z^2+3\right)
+    \right],
+\end{align}
+где \(\Phi(z)\)\nbsp{}--- функция распределения (ФР) нормального закона,
+\(\phi\)\nbsp{}--- ФПР нормального закона, \(\gamma_1\)\nbsp{}--- асимметрия,
+\(\gamma_2\)\nbsp{}--- эксцесс, \(f\)\nbsp{}--- ФПР, \(F\)\nbsp{}--- функция
+распределения (ФР). Согласно\nbsp{}cite:rozhkov1990stochastic для аппликат
+морских волн значение асимметрии выбирается на интервале
+\(0,1\leq\gamma_1\leq{0,52}]\), а значение эксцесса на интервале
+\(0,1\leq\gamma_2\leq{0,7}\). Семейство плотностей распределения при различных
+параметрах показано на рис.\nbsp{}[[fig-skew-normal-1]].
+
+#+name: fig-skew-normal-1
+#+begin_src R :file build/skew-normal-1-ru.pdf
+source(file.path("R", "common.R"))
+x <- seq(-3, 3, length.out=100)
+params <- data.frame(
+  skewness = c(0.00, 0.52, 0.00, 0.52),
+  kurtosis = c(0.00, 0.00, 0.70, 0.70),
+  linetypes = c("solid", "dashed", "dotdash", "dotted")
+)
+arma.skew_normal_1_plot(x, params)
+legend(
+  "topleft",
+  mapply(
+    function (s, k) {
+      as.expression(bquote(list(
+        gamma[1] == .(arma.fmt(s, 2)),
+        gamma[2] == .(arma.fmt(k, 2))
+      )))
+    },
+    params$skewness,
+    params$kurtosis
+  ),
+  lty = paste(params$linetypes)
+)
+#+end_src
+
+#+caption[График плотности распределения на основе РГШ]:
+#+caption: График плотности распределения закона, основанного на РГШ,
+#+caption: при различных значениях асимметрии \(\gamma_1\) и эксцесса
+#+caption: \(\gamma_2\).
+#+name: fig-skew-normal-1
+#+RESULTS: fig-skew-normal-1
+[[file:build/skew-normal-1-ru.pdf]]
+
+**** Асимметричное нормальное распределение.
+Альтернативной аппроксимацией распределения волновых аппликат служит формула
+асимметричного нормального распределения (АНР):
+\begin{align}
+    \label{eq-skew-normal-2}
+    F(z; \alpha) & = \frac{1}{2}
+   \mathrm{erfc}\left[-\frac{z}{\sqrt{2}}\right]-2 T(z,\alpha ), \nonumber \\
+    f(z; \alpha) & = \frac{e^{-\frac{z^2}{2}}}{\sqrt{2 \pi }}
+   \mathrm{erfc}\left[-\frac{\alpha z}{\sqrt{2}}\right],
+\end{align}
+где \(T\)\nbsp{}--- функция Оуэна\nbsp{}cite:owen1956tables. Эта формула не
+позволяет задать значения асимметрии и эксцесса по отдельности\nbsp{}--- оба
+значения регулируются параметром \(\alpha\). Преимущество данной формулы лишь в
+относительной простоте вычисления: эта функция встроена в некоторые программы и
+библиотеки математических функций. График функции для разных значений \(\alpha\)
+представлен на рис.\nbsp{}[[fig-skew-normal-2]].
+
+#+name: fig-skew-normal-2
+#+begin_src R :file build/skew-normal-2-ru.pdf
+source(file.path("R", "common.R"))
+x <- seq(-3, 3, length.out=100)
+alpha <- c(0.00, 0.87, 2.25, 4.90)
+params <- data.frame(
+  alpha = alpha,
+  skewness = arma.bits.skewness_2(alpha),
+  kurtosis = arma.bits.kurtosis_2(alpha),
+  linetypes = c("solid", "dashed", "dotdash", "dotted")
+)
+arma.skew_normal_2_plot(x, params)
+legend(
+  "topleft",
+  mapply(
+    function (a, s, k) {
+      as.expression(bquote(list(
+        alpha == .(arma.fmt(a, 2)),
+        gamma[1] == .(arma.fmt(s, 2)),
+        gamma[2] == .(arma.fmt(k, 2))
+      )))
+    },
+    params$alpha,
+    params$skewness,
+    params$kurtosis
+  ),
+  lty = paste(params$linetypes)
+)
+#+end_src
+
+#+caption[График плотности АНР]:
+#+caption: График плотности АНР при различных значениях коэффициента
+#+caption: асимметрии \(\alpha\).
+#+name: fig-skew-normal-2
+#+RESULTS: fig-skew-normal-2
+[[file:build/skew-normal-2.pdf]]
+
+** Форма АКФ для разных волновых профилей
+:PROPERTIES:
+:CUSTOM_ID: sec-wave-acfs
+:END:
+
+**** Аналитический метод.
+Прямой способ нахождения АКФ, соответствующей заданному профилю морской волны,
+заключается в применении теоремы Винера---Хинчина. Согласно этой теореме
+автокорреляционная функция \(K\) функции \(\zeta\) равна преобразованию Фурье от
+квадрата модуля этой функции:
 \begin{equation}
-\boxed{
-    \phi(x,z,t)
-    =
-  \InverseFourierY{
-        \frac{\Sinh{2\pi u (z+h)}}{2\pi u}
-        \frac{
-            \FourierY{ \zeta_t / \left(i f(x) - 1\right) }{u}
-        }{
-            \FourierY{ \FunSecond{\zeta(x,t)} }{u}
-        }
-    }{x},
-}
-    \label{eq-solution-2d-full}
+  K(t) = \Fourier{\left| \zeta(t) \right|^2}.
+  \label{eq-wiener-khinchin}
 \end{equation}
-где \(\FunSecond{z}\)\nbsp{}--- некоторая функция, вид которой будет определен
-в\nbsp{}[[#sec-compute-delta]] и для которой выполняется соотношение
-\(\FourierY{\FunSecond{z}}{u}=\Sinh{2\pi{u}{z}}\).
+Если заменить \(\zeta\) на формулу для волнового профиля, то это выражение даст
+аналитическую формулу для соответствующей АКФ.
 
-**** Сведение к формулам линейной теории волн.
-Справедливость полученных формул проверим, подставив в качестве \(\zeta(x,t)\)
-известные аналитические выражения для плоских волн. Символьные вычисления
-преобразований Фурье в этом разделе производились с помощью пакета
-Mathematica\nbsp{}cite:mathematica10. В линейной теории широко используется
-предположение о малости амплитуд волн, что позволяет упростить исходную систему
-уравнений\nbsp{}eqref:eq-problem-2d до
-\begin{align*}
-    & \phi_{xx}+\phi_{zz}=0,\\
-    & \zeta_t = -\phi_z & \text{на }z=\zeta(x,t),
-\end{align*}
-решение которой запишется как
-\begin{equation*}
-    \phi(x,z,t)
-    =
-    -\InverseFourierY{
-        \frac{e^{2\pi u z}}{2\pi u}
-        \FourierY{\zeta_t}{u}
-    }{x}
-    .
-\end{equation*}
-Профиль прогрессивной волны описывается формулой
-\(\zeta(x,t)=A\cos(2\pi(kx-t))\). Подстановка этого выражения
-в\nbsp{}eqref:eq-solution-2d дает равенство
-\(\phi(x,z,t)=-\frac{A}{k}\sin(2\pi(kx-t))\Sinh{2\pi{k}{z}}\). Чтобы свести его
-к формуле линейной теории волн, представим гиперболический синус в
-экспоненциальной форме и отбросим член, содержащий \(e^{-2\pi{k}{z}}\), как
-противоречащий условию \(\phi\underset{z\rightarrow-\infty}{\longrightarrow}0\).
-После взятия действительной части выражения получится известная формула линейной
-теории \(\phi(x,z,t)=\frac{A}{k}e^{2\pi{k}{z}}\sin(2\pi(kx-t))\). Аналогично,
-предположение о малости амплитуд волн позволяет упростить
-формулу\nbsp{}eqref:eq-solution-2d-full до
-\begin{equation*}
-    \phi(x,z,t)
-    =
-    -\InverseFourierY{
-        \frac{\Sinh{2\pi u (z+h)}}{2\pi u \Sinh{2\pi u h}}
-        \FourierY{\zeta_t}{u}
-    }{x}.
-\end{equation*}
-Подстановка формулы для прогрессивной плоской волны вместо \(\zeta(x,t)\) дает
-равенство
+Для трехмерного волнового профиля (два пространственных и одно временное
+измерение) аналитическая формула представляет собой многочлен высокой степени, и
+ее лучше всего вычислять с помощью программы для символьных вычислений. Затем,
+для практического применения она может быть аппроксимирована суперпозицией
+экспоненциально затухающих косинусов (именно так выглядит АКФ стационарного
+процесса АРСС\nbsp{}cite:box1976time).
+
+**** Эмпирический метод.
+Впрочем, для трехмерного случая существует более простой эмпирический метод
+нахождения формы АКФ, не требующий использования сложного программного
+обеспечения. Известно, что АКФ, представляющая собой суперпозицию
+экспоненциально затухающих косинусов, является решением уравнения Стокса для
+гравитационных волн\nbsp{}cite:boccotti1983wind. Значит, если в моделируемом
+морском волнении важна только форма волны, а не точные ее характеристики, то
+заданный волновой профиль можно просто домножить на затухающую экспоненту, чтобы
+получить подходящую АКФ. Эта АКФ не отражает параметры волн, такие как высота и
+период, зато это открывает возможность моделировать волны определенных форм,
+определяя профиль волны дискретно заданной функцией, домножая его на экспоненту
+и используя результирующую функцию в качестве АКФ. Таким образом, эмпирический
+метод неточен, но более простой по сравнению с применением теоремы
+Винера---Хинчина; он, в основном, полезен для тестирования модели АРСС.
+
+**** АКФ стоячей волны.
+Профиль трехмерной плоской стоячей волны задается как
 \begin{equation}
-    \label{eq-solution-2d-linear}
-    \phi(x,z,t)=\frac{A}{k}
-    \frac{\Sinh{2 \pi k (z+h)}}{ \Sinh{2 \pi k h} }
-    \sin(2 \pi (k x-t)),
+  \zeta(t, x, y) = A \sin (k_x x + k_y y) \sin (\sigma t).
+  \label{eq-standing-wave}
 \end{equation}
-что соответствует формуле линейной теории для конечной глубины.
-
-Различные записи решения уравнения Лапласа, в которых затухающая экспонента
-может встречаться как со знаком "+", так и со знаком "-", могут стать причиной
-разницы между формулами линейно теории и формулами, выведенными в данной работе,
-где вместо \(\sinh\) используется \(\cosh\). Выражение
-\(\frac{\Sinh{2\pi{k}(z+h)}}{\Sinh{2\pi{k}{h}}}\approx\frac{\sinh(2\pi{k}(z+h))}{\sinh(2\pi{k}{h})}\)
-превращается в строгое равенство на поверхности, и разница между правой левой
-частью увеличивается при приближении к дну водоема (для достаточно большой
-глубины ошибка вблизи поверхности жидкости незначительна). Поэтому для
-достаточно большой глубины можно использовать любую из функций (\(\cosh\) или
-\(\sinh\)) для вычисления потенциала скорости вблизи взволнованной поверхности.
+Найдем АКФ с помощью аналитического метода. Домножив формулу на затухающую
+экспоненту (поскольку преобразование Фурье определено для функции \(f\), для
+которой справедливо \(f\underset{x\rightarrow\pm\infty}{\longrightarrow}0\)),
+получим
+\begin{equation}
+  \zeta(t, x, y) =
+  A
+  \exp\left[-\alpha (|t|+|x|+|y|) \right]
+  \sin (k_x x + k_y y) \sin (\sigma t).
+  \label{eq-decaying-standing-wave}
+\end{equation}
+Затем, применяя трехмерное преобразование Фурье к обоим частям уравнения с
+помощью программы для символьных вычислений, получим многочлен высокой степени,
+который аппроксимируем выражением
+\begin{equation}
+  K(t,x,y) =
+  \gamma
+  \exp\left[-\alpha (|t|+|x|+|y|) \right]
+  \cos \beta t
+  \cos \left[ \beta x + \beta y \right].
+  \label{eq-standing-wave-acf}
+\end{equation}
+Таким образом, после применения теоремы Винера---Хинчина получаем исходную
+формулу, но с косинусами вместо синусов. Это различие важно, поскольку значение
+АКФ в точке \((0,0,0)\) равно дисперсии процесса АРСС, которое при использовании
+синусов было бы неверным.
+
+Если попытаться получить ту же самую формулу с помощью эмпирического метода, то
+выражение\nbsp{}eqref:eq-decaying-standing-wave необходимо адаптировать для
+соответствия\nbsp{}eqref:eq-standing-wave-acf. Это можно осуществить либо,
+изменяя фазу синуса, либо заменой синуса на косинус, чтобы сдвинуть максимум
+функции в начало координат.
+
+**** АКФ прогрессивной волны.
+Профиль трехмерной плоской прогрессивной волны задается как
+\begin{equation}
+  \zeta(t, x, y) = A \cos (\sigma t + k_x x + k_y y).
+  \label{eq-propagating-wave}
+\end{equation}
+Для аналитического метода повторение шагов из предыдущих двух параграфов дает
+\begin{equation}
+  K(t,x,y) =
+  \gamma
+  \exp\left[-\alpha (|t|+|x|+|y|) \right]
+  \cos\left[\beta (t+x+y) \right].
+  \label{eq-propagating-wave-acf}
+\end{equation}
+Для эмпирического метода профиль волны можно просто домножить на затухающую
+экспоненту, не изменяя положение максимума АКФ (как это требуется для стоячей
+волны).
+
+**** Сравнение изученных методов.
+Итого, аналитический метод нахождения АКФ морских волн сводится к следующим
+шагам.
+- Обеспечить затухание выражения для профиля волны на \(\pm\infty\), домножив
+  его на затухающую экспоненту.
+- Взять преобразование Фурье от квадрата модуля получившегося профиля,
+  воспользовавшись программой для символьных вычислений.
+- Аппроксимировать получившийся многочлен подходящим выражением для АКФ.
+
+Два примера этого раздела показывают, что затухающие профили стоячих и
+прогрессивных волн схожи по форме с соответствующими АКФ с тем лишь различием,
+что максимум АКФ должен быть перенесен в начало координат, чтобы сохранить
+дисперсию моделируемого процесса. Применение эмпирического метода нахождения АКФ
+сводится к следующим шагам.
+- Обеспечить затухание выражения для профиля волны на \(\pm\infty\), домножив его
+  на затухающую экспоненту.
+- Перенести максимум получившейся функции в начало координат, используя свойства
+  тригонометрических функций для сдвига фазы.
+
+** Верификация модели АРСС
+:PROPERTIES:
+:CUSTOM_ID: sec-verification
+:END:
+
+Для модели АР в
+работах\nbsp{}cite:degtyarev2011modelling,degtyarev2013synoptic,boukhanovsky1997thesis
+экспериментальным путем были верифицированы
+- распределения различных характеристик волн (высоты волн, длины волн, длины
+  гребней, период волн, уклон волн, показатель трехмерности),
+- дисперсионное соотношение,
+- сохранение интегральных характеристик для случая смешанного волнения.
+В данной работе верифицируются как модель АР, так и СС путем сравнения
+распределений различных характеристик волн.
+
+*** Верификация интегральных характеристик взволнованной поверхности
+В\nbsp{}cite:rozhkov1990stochastic авторы показывают, что некоторые
+характеристики морских волн (перечисленные в табл.\nbsp{}[[tab-weibull-shape]])
+имеют распределение Вейбулла, а подъем взволнованной поверхности\nbsp{}---
+нормальное распределение. Для верификации генерируемых моделями АР и СС
+реализаций используются спрямленные диаграммы (графики, в которых по оси \(OX\)
+откладываются квантили функции распределения, вычисленные аналитически, а по оси
+\(OY\)\nbsp{}--- вычисленные экспериментально). Если экспериментально полученное
+распределение соответствует аналитическому, то график представляет собой прямую
+линию. Концы графика могут отклоняться от прямой линии, поскольку не могут быть
+надежно получены из реализации конечной длины.
+
+#+name: tab-weibull-shape
+#+caption[Значение коэффициента формы распределения Вейбулла]:
+#+caption: Значение коэффициента формы распределения Вейбулла для различных
+#+caption: характеристик волн.
+#+attr_latex: :booktabs t
+| Характеристика          | Коэффициент формы \(k\) |
+|-------------------------+-------------------------|
+|                         | <l>                     |
+| Высота волны            | 2                       |
+| Длина волны             | 2,3                     |
+| Длина гребня волны      | 2,3                     |
+| Период волны            | 3                       |
+| Уклон волны             | 2,5                     |
+| Показатель трехмерности | 2,5                     |
+
+Верификация производится для стоячих и прогрессивных волн. Соответствующие АКФ и
+спрямленные диаграммы распределений характеристик волн представлены на
+рис.\nbsp{}[[acf-slices]],\nbsp{}[[standing-wave-distributions]],\nbsp{}[[propagating-wave-distributions]].
+
+#+name: propagating-wave-distributions
+#+begin_src R :file build/propagating-wave-qqplots-ru.pdf
+source(file.path("R", "common.R"))
+par(pty="s", mfrow=c(2, 2))
+arma.qqplot_grid(
+  file.path("build", "propagating_wave"),
+  c("elevation", "heights_y", "lengths_y", "periods"),
+  c("подъем", "высота по Y", "длина по Y", "период"),
+  xlab="x",
+  ylab="y"
+)
+#+end_src
+
+#+caption[Спрямленные диаграммы для прогрессивных волн]:
+#+caption: Спрямленные диаграммы для прогрессивных волн.
+#+name: propagating-wave-distributions
+#+RESULTS: propagating-wave-distributions
+[[file:build/propagating-wave-qqplots.pdf]]
+
+#+name: standing-wave-distributions
+#+begin_src R :file build/standing-wave-qqplots-ru.pdf
+source(file.path("R", "common.R"))
+par(pty="s", mfrow=c(2, 2))
+arma.qqplot_grid(
+  file.path("build", "standing_wave"),
+  c("elevation", "heights_y", "lengths_y", "periods"),
+  c("подъем", "высота по Y", "длина по Y", "период"),
+  xlab="x",
+  ylab="y"
+)
+#+end_src
+
+#+caption[Спрямленные диаграммы для стоячих волн]:
+#+caption: Спрямленные диаграммы для стоячих волн.
+#+name: standing-wave-distributions
+#+RESULTS: standing-wave-distributions
+[[file:build/standing-wave-qqplots-ru.pdf]]
+
+#+name: acf-slices
+#+header: :width 6 :height 9
+#+begin_src R :file build/acf-slices-ru.pdf
+source(file.path("R", "common.R"))
+propagating_acf <- read.csv(file.path("build", "propagating_wave", "acf.csv"))
+standing_acf <- read.csv(file.path("build", "standing_wave", "acf.csv"))
+par(mfrow=c(5, 2), mar=c(0,0,0,0))
+for (i in seq(0, 4)) {
+  arma.wavy_plot(standing_acf, i, zlim=c(-5,5))
+  arma.wavy_plot(propagating_acf, i, zlim=c(-5,5))
+}
+#+end_src
+
+#+caption[Временные срезы АКФ для стоячих и прогрессивных волн]:
+#+caption: Временные срезы АКФ для стоячих (слева) и прогрессивных (справа) волн.
+#+name: acf-slices
+#+RESULTS: acf-slices
+[[file:build/acf-slices-ru.pdf]]
+
+Хвосты распределений на рис.\nbsp{}[[propagating-wave-distributions]] отклоняются от
+оригинального распределения для характеристик отдельных волн, поскольку каждую
+волну необходимо извлечь из полученной взволнованной поверхности, чтобы измерить
+ее длину, период и высоту. Алгоритм, который бы гарантировал безошибочное
+извлечение всех волн, не известен, поскольку волны могут и часто накладываются
+друг на друга. Правый хвост распределения Вейбулла отклоняется больше, поскольку
+он представляет редко возникающие волны.
+
+Степень соответствия для стоячих волн (рис.\nbsp{}[[standing-wave-distributions]])
+ниже для высот и длин, примерно одинакова для подъема поверхности и выше для
+периодов волн. Более низкая степень соответствия длин и высот может быть
+результатом того, что распределения были получены эмпирически для морских волн,
+которые, в основном, являются прогрессивными, и аналогичные распределения для
+стоячих волн могут отличаться. Более высокая степень соответствия периодов волн
+является следствием того, что периоды стоячих волн извлекаются более точно,
+поскольку волн не перемещаются вне моделируемой области взволнованной
+поверхности. Одинаковая степень соответствия для подъема поверхности получается
+из-за того, что это характеристика поверхности (и соответствующего процесса АР
+или СС), и она не зависит от типа волн.
+
+*** Верификация нелинейного безынерционного преобразования
+Для того чтобы измерить влияние НБП на форму результирующей взволнованной
+поверхности, было сгенерировано три реализации:
+- реализация с Гауссовым распределением (без НБП),
+- реализация с распределением на основе ряда Грама---Шарлье (РГШ), и
+- реализация с асимметричным нормальным распределением (АНР).
+Начальные состояния ГПСЧ были заданы одинаковыми для всех запусков программы,
+чтобы модель АРСС выдавала одни и те же значения для каждой реализации. Было
+проведено два эксперимента: для стоячих и прогрессивных волн с АКФ, заданными
+формулами из раздела\nbsp{}[[#sec-wave-acfs]].
+
+В то время как эксперимент показал, что применение НБП с распределением РГШ
+увеличивает крутизну волн, то же самое нельзя сказать об асимметричном
+нормальном распределении (рис.\nbsp{}[[fig-nit]]). Использование этого распределения
+приводит к взволнованной поверхности, в которой аппликаты всегда больше или
+равны нулю. Таким образом, асимметричное нормальное распределение не подходит
+для НБП. НБП увеличивает высоту и крутизну как прогрессивных, так и стоячих
+волн. Увеличение асимметрии или эксцесса РГШ приводит в увеличению как высоты,
+так и крутизны волн. Ошибка аппроксимации АКФ (ур.\nbsp{}eqref:eq-nit-error)
+принимает значения от 0,20 для РГШ до 0,70 для АНР (табл.\nbsp{}[[tab-nit-error]]).
+
+#+name: fig-nit
+#+header: :width 7 :height 7
+#+begin_src R :file build/nit.pdf
+source(file.path("R", "nonlinear.R"))
+par(mfrow=c(2, 1), mar=c(4,4,4,0.5), family='serif')
+args <- list(
+  graphs=c('Гауссово', 'РГШ', 'АНР'),
+  linetypes=c('solid', 'dashed', 'dotted')
+)
+args$title <- 'Прогрессивные волны'
+arma.plot_nonlinear(file.path("build", "nit-propagating"), args)
+args$title <- 'Стоячие волны'
+arma.plot_nonlinear(file.path("build", "nit-standing"), args)
+#+end_src
+
+#+name: fig-nit
+#+caption[Срезы взволнованной поверхности с различными распределениями]:
+#+caption: Срезы взволнованной поверхности с различными распределениями
+#+caption: волновых аппликат (Гауссово, РГШ и асимметричное нормальное).
+#+RESULTS: fig-nit
+[[file:build/nit.pdf]]
+
+#+name: tab-nit-error
+#+caption[Ошибки аппроксимации АКФ]:
+#+caption: Ошибки аппроксимации АКФ (разность дисперсий) для различных
+#+caption: распределений волновых аппликат. \(N\)\nbsp{}--- количество
+#+caption: коэффициентов аппроксимации АКФ.
+#+attr_latex: :booktabs t
+| Тип волн      | Распред. | \(\gamma_1\) | \(\gamma_2\) | \(\alpha\) | Ошибка | \(N\) | Высота волн |
+|---------------+----------+--------------+--------------+------------+--------+-------+-------------|
+|               |          |          <r> |          <r> |        <r> |    <r> |       |         <r> |
+| прогрессивные | Гауссово |              |              |            |        |       |        2,41 |
+| прогрессивные | РГШ      |         2,25 |          0,4 |            |   0,20 |     2 |        2,75 |
+| прогрессивные | АНР      |              |              |          1 |   0,70 |     3 |        1,37 |
+| стоячие       | Гауссово |              |              |            |        |       |        1,73 |
+| стоячие       | РГШ      |         2,25 |          0,4 |            |   0,26 |     2 |        1,96 |
+| стоячие       | АНР      |              |              |          1 |   0,70 |     3 |        0,94 |
+
+Таким образом, единственный тестовый сценарий, который показал приемлемые
+результаты\nbsp{}--- это реализации с распределением на основе РГШ для
+прогрессивных и стоячих волн. АНР искажает взволнованную поверхность для обоих
+типов волн. Реализации с распределением на основе РГШ характеризуются большой
+ошибкой аппроксимации АКФ, что приводит к увеличению высоты волн. Причина
+большой ошибки заключается в неточности аппроксимации РГШ, которая не сходится
+для всевозможных функций\nbsp{}cite:wallace1958asymptotic. Несмотря на большую
+ошибку, изменение высоты волн невелико (табл.\nbsp{}[[tab-nit-error]]).
+
+*** Нефизическая природа модели
+В силу своей нефизической природы модель АРСС не включает в себя понятие морской
+волны; вместо этого она моделирует взволнованную поверхность как единое целое.
+Движения отдельных волн и их форма часто получаются грубыми, а точное количество
+генерируемых волн неизвестно. Несмотря на это, интегральные характеристики
+взволнованной поверхности соответствуют реальным морским волнам.
 
-Сведение формул\nbsp{}eqref:eq-solution-2d и\nbsp{}eqref:eq-solution-2d-full к
-формулам линейной теории волн показывает, что формула\nbsp{}eqref:eq-solution-2d
-для жидкости бесконечной глубины не подходит для вычисления потенциала скорости
-с использованием метода Фурье, т.к.\nbsp{}не обладает необходимой для
-преобразования Фурье симметрией. Однако, для такого случая можно использовать
-формулу для конечной глубины, полагая \(h\) равным некоторому характерному
-значению глубины. Для стоячих волн сведение к формулам линейной теории
-происходит с аналогичными предположениями.
+Теоретически, профили самих морских волн могут быть использованы в качестве АКФ,
+если предварительно обеспечить их экспоненциальное затухание. Это может
+позволить генерировать волны произвольных профилей и является одной из тем
+дальнейших исследований.
 
-*** Трехмерное поле скоростей
-В трех измерениях исходная система уравнений\nbsp{}eqref:eq-problem
-переписывается как
+* Поле давлений под дискретно заданной взволнованной поверхностью
+** Известные формулы определения поля давлений
+**** Теория волн малых амплитуд.
+В\nbsp{}cite:stab2012,degtyarev1998modelling,degtyarev1997analysis дается
+решение обратной задачи гидродинамики для случая идеальной несжимаемой жидкости
+в рамках теории волн малых амплитуд (в предположении, что длина волны много
+больше ее высоты: \(\lambda \gg h\)). В этом случае обратная задача линейна и
+сводится к уравнению Лапласа со смешанным граничным условием, а уравнение
+движения используется только для нахождения давлений по известным значениям
+производных потенциала скорости. Предположение о малости амплитуд волн означает
+слабое изменение локального волнового числа во времени и пространстве по
+сравнению с подъемом (аппликатой) взволнованной поверхности. Это позволяет
+вычислить производную подъема поверхности по \(z\) как \(\zeta_z=k\zeta\), где
+\(k\)\nbsp{}--- волновое число. В двухмерном случае решение записывается явной
+формулой
 \begin{align}
-    \label{eq-problem-3d}
-    & \phi_{xx} + \phi_{yy} + \phi_{zz} = 0,\\
-    & \zeta_t + \zeta_x\phi_x + \zeta_y\phi_y
-    =
-    \frac{\zeta_x}{\SqrtZeta{1 + \zeta_x^2 + \zeta_y^2}} \phi_x
-    +\frac{\zeta_y}{\SqrtZeta{1 + \zeta_x^2 + \zeta_y^2}} \phi_y
-    - \phi_z, & \text{на }z=\zeta(x,y,t).\nonumber
+    \left.\frac{\partial\phi}{\partial x}\right|_{x,t}= &
+        -\frac{1}{\sqrt{1+\alpha^{2}}}e^{-I(x)}
+            \int\limits_{0}^x\frac{\partial\dot{\zeta}/\partial
+                z+\alpha\dot{\alpha}}{\sqrt{1+\alpha^{2}}}e^{I(x)}dx,\label{eq-old-sol-2d}\\
+    I(x)= & \int\limits_{0}^x\frac{\partial\alpha/\partial z}{1+\alpha^{2}}dx,\nonumber
 \end{align}
-Для ее решения также воспользуемся методом Фурье. Возьмем преобразование Фурье
-от обоих частей уравнений Лапласа и получим
-\begin{equation*}
-    -4 \pi^2 \left( u^2 + v^2 + w^2 \right)
-    \FourierY{\phi(x,y,z)}{u,v,w} = 0,
-\end{equation*}
-откуда имеем \(w=\pm{i}\sqrt{u^2+v^2}\). Решение уравнения будем искать в виде
-обратного преобразования Фурье
-\(\phi(x,y,z)=\InverseFourierY{E(u,v,w)}{x,y,z}\). Подставляя
-\(w=i\sqrt{u^2+v^2}=i\Kveclen\) в исходную формулу, получаем
-\begin{equation*}
-    \phi(x,y,z) = \InverseFourierY{
-        \left(
-            C_1 e^{2\pi \Kveclen z}
-            -C_2 e^{-2\pi \Kveclen z}
-        \right)
-        E(u,v)
-    }{x,y}.
-\end{equation*}
-Подставляя \(\phi\) в условие на дне водоема аналогично двухмерному случаю,
-получаем
-\begin{equation}
-    \label{eq-guessed-sol-3d}
-    \phi(x,y,z) = \InverseFourierY{
-        \Sinh{2\pi \Kveclen (z+h)} E(u,v)
-    }{x,y}.
-\end{equation}
-Подставляя выражение для \(\phi\) в граничное условие, получим
-\begin{equation*}
-    \arraycolsep=1.4pt
-    \begin{array}{rl}
-        \zeta_t = & i f_1(x,y) \InverseFourierY{2 \pi u \Sinh{2\pi \Kveclen (z+h)}E(u,v)}{x,y} \\
-        + & i f_2(x,y) \InverseFourierY{2 \pi v \Sinh{2\pi \Kveclen (z+h)}E(u,v)}{x,y} \\
-        - & \InverseFourierY{2 \pi \sqrt{u^2+v^2} \Sinh{2\pi \Kveclen (z+h)}E(u,v)}{x,y}
-    \end{array}
-\end{equation*}
-где \(f_1(x,y)={\zeta_x}/{\SqrtZeta{1+\zeta_x^2+\zeta_y^2}}-\zeta_x\) и
-\(f_2(x,y)={\zeta_y}/{\SqrtZeta{1+\zeta_x^2+\zeta_y^2}}-\zeta_y\).
+где \(\alpha\)\nbsp{}--- уклоны волн. В трехмерном случае решение записывается в
+виде эллиптического дифференциального уравнения в частных производных (ДУЧП):
+\begin{align*}
+    & \frac{\partial^2 \phi}{\partial x^2} \left( 1 + \alpha_x^2 \right) +
+    \frac{\partial^2 \phi}{\partial y^2} \left( 1 + \alpha_y^2 \right) +
+    2\alpha_x\alpha_y \frac{\partial^2 \phi}{\partial x \partial y} + \\
+    & \left(
+        \frac{\partial \alpha_x}{\partial z} +
+        \alpha_x \frac{\partial \alpha_x}{\partial x} +
+        \alpha_y \frac{\partial \alpha_x}{\partial y}
+    \right) \frac{\partial \phi}{\partial x} + \\
+    & \left(
+        \frac{\partial \alpha_y}{\partial z} +
+        \alpha_x \frac{\partial \alpha_y}{\partial x} +
+        \alpha_y \frac{\partial \alpha_y}{\partial y}
+    \right) \frac{\partial \phi}{\partial y} + \\
+    & \frac{\partial \dot{\zeta}}{\partial z} +
+    \alpha_x \dot{\alpha_x} + \alpha_y \dot{\alpha_y} = 0.
+\end{align*}
+Уравнение предполагается решать численно путем сведения к разностному.
 
-Также как и в разделе\nbsp{}[[#sec-pressure-2d]] мы предполагаем, что
-\(\Sinh{2\pi{u}(z+h)}\approx\SinhX{2\pi{u}(z+h)}\) вблизи свободной поверхности,
-однако в трехмерном случае этого недостаточно для решения задачи. Для того чтобы
-получить явную формулу для коэффициентов \(E\), мы должны предположить, что
-преобразования Фурье в равенстве имеют радиально симметричные ядра,
-т.е.\nbsp{}заменить \(u\) и \(v\) на \(\Kveclen\). Есть два момента,
-поддерживающих это предположение. Во-первых, в численной реализации
-интегрирование ведется по положительным волновым числам, так что знак \(u\) и
-\(v\) не влияет на решение. Во-вторых, скорость роста \(\cosh\) в ядре интеграла
-значительно выше, чем скорость роста \(u\) или \(\Kveclen\), так что замена
-слабо влияет на величину решения. Несмотря на эти два момента, использование
-более математически строго подхода было бы предпочтительнее.
+Как будет показано в разделе\nbsp{}[[#sec-compare-formulae]]
+формула\nbsp{}eqref:eq-old-sol-2d расходится при попытке вычислить поле
+скоростей для волн больших амплитуд, а значит не может быть использована
+совместно с моделью морского волнения, генерирующей волны произвольных амплитуд.
 
-Выполняя замену, применяя преобразование Фурье к обеим частям равенства, и,
-подставляя результат в\nbsp{}eqref:eq-guessed-sol-3d, получаем выражение для
-\(\phi\):
+**** Линеаризация граничного условия.
+:PROPERTIES:
+:CUSTOM_ID: linearisation
+:END:
+
+Модель ЛХ позволяет вывести явную формулу для поля скоростей путем линеаризации
+кинематического граничного условия. Формула для потенциала скорости запишется
+как
 \begin{equation*}
-    \label{eq-phi-3d}
-    \phi(x,y,z,t) = \InverseFourierY{
-        \frac{ \Sinh{\smash{2\pi \Kveclen (z+h)}} }{ 2\pi\Kveclen }
-        \frac{ \FourierY{ \zeta_t / \left( i f_1(x,y) + i f_2(x,y) - 1 \right)}{u,v} }
-        { \FourierY{\mathcal{D}_3\left( x,y,\zeta\left(x,y\right) \right)}{u,v} }
-    }{x,y},
+\phi(x,y,z,t) = \sum_n \frac{c_n g}{\omega_n}
+     e^{\sqrt{u_n^2+v_n^2} z}
+     \sin(u_n x + v_n y - \omega_n t + \epsilon_n).
 \end{equation*}
-где
-\(\FourierY{\mathcal{D}_3\left(x,y,z\right)}{u,v}=\Sinh{\smash{2\pi\Kveclen{}z}}\).
+Формула дифференцируется для получения производных потенциала, которые
+подставляются в динамическое граничное условие для определения поля давлений.
 
-*** Заключение
-Полученные в данном разделе формулы позволяют произвести численный расчет поля
-скоростей (а значит и поля давлений) вблизи дискретно или математически заданной
-взволнованной морской поверхности, минуя предположения линейной теории и теории
-волн малых амплитуд.
+** Определение поля давлений под дискретно заданной взволнованной поверхностью
+Аналитические решения граничных задач для классических уравнений часто
+используются для исследования различных свойств уравнений, и для таких
+исследований запись формулы общего решения неудобна ввиду своей сложности и
+наличия интегралов от неизвестных функций. Одним из методов нахождения
+аналитических решений ДУЧП является метод Фурье. Основой метода служит
+преобразование Фурье, применение которого к некоторым ДУЧП позволяет свести их к
+алгебраическому, а его решение записывается как обратное преобразование Фурье от
+некоторой функции (которая может содержать преобразования Фурье от других
+функций). Поскольку эти преобразования не всегда можно записать аналитически, то
+вместо этого ищутся частные решения задачи и анализируется их поведение в
+различных областях. В то же время, численный расчет дискретных преобразований
+Фурье возможен для любой дискретно заданной функции, используя алгоритмы БПФ.
+Эти алгоритмы используют симметрию комплексных экспонент для понижения
+асимптотической сложности с \(\mathcal{O}(n^2)\) до \(\mathcal{O}(n\log_{2}n)\).
+Таким образом, даже если общее решение содержит преобразования Фурье от
+неизвестных функций, они все равно могут быть взяты численно, а использование
+алгоритмов БПФ делает этот подход эффективным.
 
-** Моделирование нелинейности морских волн
-Модель АРСС позволяет учесть асимметричность распределения волновых аппликат,
-т.е.\nbsp{}генерировать морские волны, закон распределения аппликат которых
-имеет ненулевой эксцесс и асимметрию. Такой закон распределения характерен для
-реальных морских волн\nbsp{}cite:longuet1963nonlinear. Асимметричность волн
-моделируется с помощью нелинейного безынерционного преобразования (НБП)
-случайного процесса, однако, любое нелинейное преобразование случайного процесса
-приводит к преобразованию его АКФ. Для того чтобы подавить этот эффект,
-необходимо предварительно преобразовать АКФ, как показано
-в\nbsp{}cite:boukhanovsky1997thesis.
+Альтернативным подходом к решению ДУЧП является их сведение к разностным
+уравнениям, решаемым с помощью построения различных численных схем. При этом
+решение получается приближенным, а асимптотическая сложность соответствующих
+алгоритмов сопоставима со сложностью алгоритма БПФ. Например, для стационарного
+эллиптического уравнения в частных производных строится неявная численная схема;
+эта схема расчитывается итерационным методом, на каждом шаге которого ищется
+решение трехдиагональной или пятидиагональной СЛАУ методом прогонки.
+Асимптотическая сложность алгоритма составляет \(\mathcal{O}({n}{m})\), где
+\(n\)\nbsp{}--- количество точек на сетке взволнованной поверхности, а
+\(m\)\nbsp{}--- число итераций.
 
-**** Преобразование взволнованной поверхности.
-Формула \(z=f(y)\) преобразования взволнованной поверхности к необходимому
-одномерному закону распределения \(F(z)\) получается путем решения нелинейного
-трансцендентного уравнения \(F(z) = \Phi(y)\), где \(\Phi(y)\)\nbsp{}--- функция
-одномерного нормального закона распределения. Поскольку функция распределения
-аппликат морских волн часто задается некоторой аппроксимацией, основанной на
-натурных данных, то это уравнение целесообразно решать численно в каждой точке
-\(y_k|_{k=0}^N\) сетки сгенерированной поверхности относительно \(z_k\). Тогда
-уравнение запишется в виде
+Несмотря на широкое распространение, итеративные алгоритмы неэффективно
+отображаются на архитектуру параллельных машин ввиду неизбежной синхронизации
+процессов после каждой итерации; в частности, отображение на сопроцессоры может
+включать в себя копирование данных на сопроцессор и обратно на каждой итерации,
+что отрицательно сказывается на их производительности. В то же время, наличие
+большого количества преобразований Фурье в решении является скорее
+преимуществом, чем недостатком. Во-первых, решения, полученные с помощью метода
+Фурье, явные, а значит хорошо масштабируются на большое количество параллельно
+работающих вычислительных ядер с использованием простейших приемов параллельного
+программирования. Во-вторых, для алгоритмов БПФ существуют готовые
+оптимизированные реализации для различных архитектур процессоров и сопроцессоров
+(GPU, MIC). Эти преимущества обусловливают использование метода Фурье для
+получения явного решения задачи определения давлений под взволнованной морской
+поверхностью.
+
+*** Двухмерное поле скоростей
+:PROPERTIES:
+:CUSTOM_ID: sec-pressure-2d
+:END:
+
+**** Формула для жидкости бесконечной глубины.
+Задача Робена для уравнения Лапласа в двух измерениях записывается как
+\begin{align}
+    \label{eq-problem-2d}
+    & \phi_{xx}+\phi_{zz}=0,\\
+    & \zeta_t + \zeta_x\phi_x = \frac{\zeta_x}{\sqrt{1 + \zeta_x^2}} \phi_x - \phi_z, & \text{на }z=\zeta(x,t).\nonumber
+\end{align}
+Для ее решения воспользуемся методом Фурье. Возьмем преобразование Фурье от
+обоих частей уравнений Лапласа и получим
+\begin{equation*}
+    -4 \pi^2 \left( u^2 + v^2 \right)
+    \FourierY{\phi(x,z)}{u,v} = 0,
+\end{equation*}
+откуда имеем \(v = \pm i u\). Здесь и далее будет использоваться следующая
+симметричная форма преобразования Фурье:
+\begin{equation*}
+    \FourierY{f(x,y)}{u,v} =
+    \iint\limits_{-\infty}^{\phantom{--}\infty}
+    f(x,y)
+    e^{-2\pi i (x u + y v)}
+    dx dy.
+\end{equation*}
+Решение уравнения будем искать в виде обратного преобразования Фурье
+\(\phi(x,z)=\InverseFourierY{E(u,v)}{x,z}\). Подставляя[fn::Выражение \(v={-i}{u}\)
+не подходит в данной задаче, поскольку потенциал скорости должен стремиться к
+нулю с увеличением глубины до бесконечности.] \(v={i}{u}\) в формулу, решение
+перепишется как
 \begin{equation}
-    \label{eq-distribution-transformation}
-    F(z_k)
-    =
-    \frac{1}{\sqrt{2\pi}}
-    \int\limits_0^{y_k} \exp\left[ -\frac{t^2}{2} \right] dt
-    .
+    \label{eq-guessed-sol-2d}
+    \phi(x,z) = \InverseFourierY{e^{2\pi u z}E(u)}{x}.
 \end{equation}
-Поскольку функции распределения монотонны, для решения этого уравнения
-используется простейший численный метод половинного деления (метод бисекции).
-
-**** Предварительное преобразование АКФ.
-Для преобразования АКФ \(\gamma_z\) процесса ее необходимо разложить в ряд по
-полиномам Эрмита (ряд Грама---Шарлье)
+Для того чтобы подстановка \(z=\zeta(x,t)\) не помешала использованию
+преобразований Фурье в решении, перепишем\nbsp{}eqref:eq-guessed-sol-2d в виде
+свертки:
 \begin{equation*}
-    \gamma_z \left( \vec u \right)
+    \phi(x,z)
     =
-    \sum\limits_{m=0}^{\infty}
-    C_m^2 \frac{\gamma_y^m \left( \vec u \right)}{m!},
+    \Fun{z}
+    \ast
+    \InverseFourierY{E(u)}{x},
 \end{equation*}
-где
+где \(\Fun{z}\)\nbsp{}--- некоторая функция, вид которой будет определен в
+разделе\nbsp{}[[#sec-compute-delta]] и для которой выполняется соотношение
+\(\FourierY{\Fun{z}}{u}=e^{2\pi{u}{z}}\). Подставляя выражение для \(\phi\) в
+граничное условие, получим
 \begin{equation*}
-    C_m = \frac{1}{\sqrt{2\pi}}
-  \int\limits_{0}^\infty
-    f(y) H_m(y) \exp\left[ -\frac{y^2}{2} \right],
+    \zeta_t
+    =
+    \left( i f(x) - 1 \right)
+    \left[
+        \Fun{z}
+        \ast
+        \InverseFourierY{2\pi u E(u)}{x}
+    \right],
 \end{equation*}
-\(H_m\)\nbsp{}--- полином Эрмита, а \(f(y)\)\nbsp{}--- решение
-уравнения\nbsp{}eqref:eq-distribution-transformation. Воспользовавшись
-полиномиальной аппроксимацией \(f(y) \approx \sum\limits_i d_i y^i\) и
-аналитическими выражениями для полиномов Эрмита, формулу определения
-коэффициентов можно упростить, используя следующее равенство:
+где \(f(x) = {\zeta_x}/{\sqrt{1 + \zeta_x^2}} - \zeta_x\). Применяя преобразование
+Фурье к обеим частям, получаем выражение для коэффициентов \(E\):
 \begin{equation*}
-    \frac{1}{\sqrt{2\pi}}
-    \int\limits_\infty^\infty
-    y^k \exp\left[ -\frac{y^2}{2} \right]
-    =
-    \begin{cases}
-        (k-1)!! & \text{для четных }k,\\
-        0       & \text{для нечетных }k.
-    \end{cases}
+    E(u) =
+    \frac{1}{2\pi u}
+    \frac{
+    \FourierY{\zeta_t / \left(i f(x) - 1\right)}{u}
+    }{
+    \FourierY{\Fun{z}}{u}
+    }
 \end{equation*}
-Оптимальное количество коэффициентов \(C_m\) определяется путем вычисления их
-последовательно и критерий прекращения счета определяется совпадением дисперсий
-обоих полей с требуемой точностью \(\epsilon\):
-\begin{equation}
-    \label{eq-nit-error}
-    \left| \Var{z} - \sum\limits_{k=0}^m
-    \frac{C_k^2}{k!} \right| \leq \epsilon.
-\end{equation}
-
-В\nbsp{}cite:boukhanovsky1997thesis автор предлагает использовать полиномиальную
-аппроксимацию для \(f(y)\) также для преобразования поверхности, однако на
-практике в реализации взволнованной поверхности часто находятся точки,
-выпадающие за промежуток на котором построена аппроксимация, что приводит к
-резкому уменьшению ее точности. В этих точках
-уравнение\nbsp{}eqref:eq-distribution-transformation эффективнее решать методом
-бисекции. Использование полиномиальной аппроксимации в формулах для
-коэффициентов ряда Грама---Шарлье не приводит к аналогичным ошибкам.
-
-* Численные методы и результаты экспериментов
-** Форма АКФ для разных волновых профилей
-:PROPERTIES:
-:CUSTOM_ID: sec-wave-acfs
-:END:
-
-**** Аналитический метод.
-Прямой способ нахождения АКФ, соответствующей заданному профилю морской волны,
-заключается в применении теоремы Винера---Хинчина. Согласно этой теореме
-автокорреляционная функция \(K\) функции \(\zeta\) равна преобразованию Фурье от
-квадрата модуля этой функции:
+Выполняя подстановку \(z=\zeta(x,t)\) и подставляя полученное выражение
+в\nbsp{}eqref:eq-guessed-sol-2d, получаем окончательное выражение для
+\(\phi(x,z)\):
 \begin{equation}
-  K(t) = \Fourier{\left| \zeta(t) \right|^2}.
-  \label{eq-wiener-khinchin}
+    \label{eq-solution-2d}
+    \boxed{
+        \phi(x,z)
+        =
+        \InverseFourierY{
+            \frac{e^{2\pi u z}}{2\pi u}
+            \frac{
+            \FourierY{ \zeta_t / \left(i f(x) - 1\right) }{u}
+            }{
+            \FourierY{ \Fun{\zeta(x,t)} }{u}
+            }
+        }{x}.
+    }
 \end{equation}
-Если заменить \(\zeta\) на формулу для волнового профиля, то это выражение даст
-аналитическую формулу для соответствующей АКФ.
-
-Для трехмерного волнового профиля (два пространственных и одно временное
-измерение) аналитическая формула представляет собой многочлен высокой степени, и
-ее лучше всего вычислять с помощью программы для символьных вычислений. Затем,
-для практического применения она может быть аппроксимирована суперпозицией
-экспоненциально затухающих косинусов (именно так выглядит АКФ стационарного
-процесса АРСС\nbsp{}cite:box1976time).
 
-**** Эмпирический метод.
-Впрочем, для трехмерного случая существует более простой эмпирический метод
-нахождения формы АКФ, не требующий использования сложного программного
-обеспечения. Известно, что АКФ, представляющая собой суперпозицию
-экспоненциально затухающих косинусов, является решением уравнения Стокса для
-гравитационных волн\nbsp{}cite:boccotti1983wind. Значит, если в моделируемом
-морском волнении важна только форма волны, а не точные ее характеристики, то
-заданный волновой профиль можно просто домножить на затухающую экспоненту, чтобы
-получить подходящую АКФ. Эта АКФ не отражает параметры волн, такие как высота и
-период, зато это открывает возможность моделировать волны определенных форм,
-определяя профиль волны дискретно заданной функцией, домножая его на экспоненту
-и используя результирующую функцию в качестве АКФ. Таким образом, эмпирический
-метод неточен, но более простой по сравнению с применением теоремы
-Винера---Хинчина; он, в основном, полезен для тестирования модели АРСС.
+Множитель \(e^{2\pi u z}/(2\pi u)\) делает график функции, от которой берется
+обратное преобразования Фурье, несимметричным относительно оси \(OY\). Это
+затрудняет применение БПФ, поскольку оно требует периодичную функцию,
+принимающую нулевые значения на концах промежутка. В то же время, использование
+численного интегрирования вместо БПФ не позволит получить преимущество над
+решением всей системы уравнений с помощью разностных схем. Эту проблему можно
+обойти, используя формулу\nbsp{}eqref:eq-solution-2d-full для жидкости конечной
+глубины с заведомо большим значением глубины водоема \(h\). Вывод формулы дан в
+следующем разделе.
 
-**** АКФ стоячей волны.
-Профиль трехмерной плоской стоячей волны задается как
-\begin{equation}
-  \zeta(t, x, y) = A \sin (k_x x + k_y y) \sin (\sigma t).
-  \label{eq-standing-wave}
-\end{equation}
-Найдем АКФ с помощью аналитического метода. Домножив формулу на затухающую
-экспоненту (поскольку преобразование Фурье определено для функции \(f\), для
-которой справедливо \(f\underset{x\rightarrow\pm\infty}{\longrightarrow}0\)),
-получим
+**** Формула для жидкости конечной глубины.
+На дне водоема вертикальная составляющая скорости перемещения жидкости должна
+равняться нулю, т.е.\nbsp{}\(\phi_z=0\) на \(z=-h\), где \(h\)\nbsp{}--- глубина
+водоема. В этом случае пренебречь равенством \(v = -i u\), полученным из
+уравнения Лапласа, нельзя, и решение ищется в виде
 \begin{equation}
-  \zeta(t, x, y) =
-  A
-  \exp\left[-\alpha (|t|+|x|+|y|) \right]
-  \sin (k_x x + k_y y) \sin (\sigma t).
-  \label{eq-decaying-standing-wave}
+    \phi(x,z)
+    =
+    \InverseFourierY{
+        \left( C_1 e^{2\pi u z} + C_2 e^{-2\pi u z} \right)
+        E(u)
+    }{x}.
+    \label{eq-guessed-sol-2d-full}
 \end{equation}
-Затем, применяя трехмерное преобразование Фурье к обоим частям уравнения с
-помощью программы для символьных вычислений, получим многочлен высокой степени,
-который аппроксимируем выражением
+Подставляя \(\phi\) в условие на дне водоема, получим
+\begin{equation*}
+    C_1 e^{-2\pi u h} - C_2 e^{2\pi u h} = 0,
+\end{equation*}
+откуда имеем \(C_1=\frac{1}{2}C{e}^{2\pi{u}{h}}\) и
+\(C_2=-\frac{1}{2}C{e}^{-2\pi{u}{h}}\). Константа \(C\) здесь произвольна,
+поскольку при подстановке станет частью неизвестных коэффициентов \(E(u)\).
+Подставляя полученные выражения для \(C_1\) и \(C_2\)
+в\nbsp{}eqref:eq-guessed-sol-2d-full, получаем выражение
+\begin{equation*}
+    \phi(x,z) = \InverseFourierY{ \Sinh{2\pi u (z+h)} E(u) }{x}.
+\end{equation*}
+Подставляя \(\phi\) в граничное условие на свободной поверхности, получаем
+\begin{equation*}
+    \zeta_t = f(x) \InverseFourierY{ 2\pi i u \Sinh{2\pi u (z+h)} E(u) }{x}
+            - \InverseFourierY{ 2\pi u \SinhX{2\pi u (z+h)} E(u) }{x}.
+\end{equation*}
+Здесь \(\sinh\) и \(\cosh\) дают схожие результаты вблизи свободной поверхности, и,
+поскольку эта область является наиболее интересной с точки зрения практического
+применения, положим \(\Sinh{2\pi{u}(z+h)}\approx\SinhX{2\pi{u}(z+h)}\). Выполняя
+аналогичные предыдущему разделу операции, получаем окончательное выражение для
+\(\phi(x,z)\):
 \begin{equation}
-  K(t,x,y) =
-  \gamma
-  \exp\left[-\alpha (|t|+|x|+|y|) \right]
-  \cos \beta t
-  \cos \left[ \beta x + \beta y \right].
-  \label{eq-standing-wave-acf}
+\boxed{
+    \phi(x,z,t)
+    =
+  \InverseFourierY{
+        \frac{\Sinh{2\pi u (z+h)}}{2\pi u}
+        \frac{
+            \FourierY{ \zeta_t / \left(i f(x) - 1\right) }{u}
+        }{
+            \FourierY{ \FunSecond{\zeta(x,t)} }{u}
+        }
+    }{x},
+}
+    \label{eq-solution-2d-full}
 \end{equation}
-Таким образом, после применения теоремы Винера---Хинчина получаем исходную
-формулу, но с косинусами вместо синусов. Это различие важно, поскольку значение
-АКФ в точке \((0,0,0)\) равно дисперсии процесса АРСС, которое при использовании
-синусов было бы неверным.
-
-Если попытаться получить ту же самую формулу с помощью эмпирического метода, то
-выражение\nbsp{}eqref:eq-decaying-standing-wave необходимо адаптировать для
-соответствия\nbsp{}eqref:eq-standing-wave-acf. Это можно осуществить либо,
-изменяя фазу синуса, либо заменой синуса на косинус, чтобы сдвинуть максимум
-функции в начало координат.
+где \(\FunSecond{z}\)\nbsp{}--- некоторая функция, вид которой будет определен
+в\nbsp{}[[#sec-compute-delta]] и для которой выполняется соотношение
+\(\FourierY{\FunSecond{z}}{u}=\Sinh{2\pi{u}{z}}\).
 
-**** АКФ прогрессивной волны.
-Профиль трехмерной плоской прогрессивной волны задается как
-\begin{equation}
-  \zeta(t, x, y) = A \cos (\sigma t + k_x x + k_y y).
-  \label{eq-propagating-wave}
-\end{equation}
-Для аналитического метода повторение шагов из предыдущих двух параграфов дает
+**** Сведение к формулам линейной теории волн.
+Справедливость полученных формул проверим, подставив в качестве \(\zeta(x,t)\)
+известные аналитические выражения для плоских волн. Символьные вычисления
+преобразований Фурье в этом разделе производились с помощью пакета
+Mathematica\nbsp{}cite:mathematica10. В линейной теории широко используется
+предположение о малости амплитуд волн, что позволяет упростить исходную систему
+уравнений\nbsp{}eqref:eq-problem-2d до
+\begin{align*}
+    & \phi_{xx}+\phi_{zz}=0,\\
+    & \zeta_t = -\phi_z & \text{на }z=\zeta(x,t),
+\end{align*}
+решение которой запишется как
+\begin{equation*}
+    \phi(x,z,t)
+    =
+    -\InverseFourierY{
+        \frac{e^{2\pi u z}}{2\pi u}
+        \FourierY{\zeta_t}{u}
+    }{x}
+    .
+\end{equation*}
+Профиль прогрессивной волны описывается формулой
+\(\zeta(x,t)=A\cos(2\pi(kx-t))\). Подстановка этого выражения
+в\nbsp{}eqref:eq-solution-2d дает равенство
+\(\phi(x,z,t)=-\frac{A}{k}\sin(2\pi(kx-t))\Sinh{2\pi{k}{z}}\). Чтобы свести его
+к формуле линейной теории волн, представим гиперболический синус в
+экспоненциальной форме и отбросим член, содержащий \(e^{-2\pi{k}{z}}\), как
+противоречащий условию \(\phi\underset{z\rightarrow-\infty}{\longrightarrow}0\).
+После взятия действительной части выражения получится известная формула линейной
+теории \(\phi(x,z,t)=\frac{A}{k}e^{2\pi{k}{z}}\sin(2\pi(kx-t))\). Аналогично,
+предположение о малости амплитуд волн позволяет упростить
+формулу\nbsp{}eqref:eq-solution-2d-full до
+\begin{equation*}
+    \phi(x,z,t)
+    =
+    -\InverseFourierY{
+        \frac{\Sinh{2\pi u (z+h)}}{2\pi u \Sinh{2\pi u h}}
+        \FourierY{\zeta_t}{u}
+    }{x}.
+\end{equation*}
+Подстановка формулы для прогрессивной плоской волны вместо \(\zeta(x,t)\) дает
+равенство
 \begin{equation}
-  K(t,x,y) =
-  \gamma
-  \exp\left[-\alpha (|t|+|x|+|y|) \right]
-  \cos\left[\beta (t+x+y) \right].
-  \label{eq-propagating-wave-acf}
+    \label{eq-solution-2d-linear}
+    \phi(x,z,t)=\frac{A}{k}
+    \frac{\Sinh{2 \pi k (z+h)}}{ \Sinh{2 \pi k h} }
+    \sin(2 \pi (k x-t)),
 \end{equation}
-Для эмпирического метода профиль волны можно просто домножить на затухающую
-экспоненту, не изменяя положение максимума АКФ (как это требуется для стоячей
-волны).
-
-**** Сравнение изученных методов.
-Итого, аналитический метод нахождения АКФ морских волн сводится к следующим
-шагам.
-- Обеспечить затухание выражения для профиля волны на \(\pm\infty\), домножив
-  его на затухающую экспоненту.
-- Взять преобразование Фурье от квадрата модуля получившегося профиля,
-  воспользовавшись программой для символьных вычислений.
-- Аппроксимировать получившийся многочлен подходящим выражением для АКФ.
-
-Два примера этого раздела показывают, что затухающие профили стоячих и
-прогрессивных волн схожи по форме с соответствующими АКФ с тем лишь различием,
-что максимум АКФ должен быть перенесен в начало координат, чтобы сохранить
-дисперсию моделируемого процесса. Применение эмпирического метода нахождения АКФ
-сводится к следующим шагам.
-- Обеспечить затухание выражения для профиля волны на \(\pm\infty\), домножив его
-  на затухающую экспоненту.
-- Перенести максимум получившейся функции в начало координат, используя свойства
-  тригонометрических функций для сдвига фазы.
-
-** Дополнительные формулы, методы и алгоритмы для модели АРСС
-:PROPERTIES:
-:CUSTOM_ID: sec-arma-algorithms
-:END:
-
-*** Аппроксимация распределения аппликат
-Одним из параметров генератора взволнованной морской поверхности служит функция
-плотности распределения (ФПР) аппликат этой поверхности. Она задается либо
-полиномиальной аппроксимацией натурных данных, либо аналитически.
-
-**** Разложение в ряд Грама---Шарлье.
-В\nbsp{}cite:huang1980experimental было экспериментально показано, что
-распределение аппликат морской поверхности отличается от нормального ненулевым
-эксцессом и асимметрией. В\nbsp{}cite:rozhkov1996theory показано, что такое
-распределение раскладывается в ряд Грама---Шарлье (РГШ):
-\begin{align}
-    \label{eq-skew-normal-1}
-    & F(z; \mu=0, \sigma=1, \gamma_1, \gamma_2) \approx
-    \Phi(z; \mu, \sigma) \frac{2 + \gamma_2}{2}
-    - \frac{2}{3} \phi(z; \mu, \sigma)
-    \left(\gamma_2 z^3+\gamma_1
-    \left(2 z^2+1\right)\right) \nonumber
-    \\
-    & f(z; \gamma_1, \gamma_2) \approx
-    \frac{1}{\sigma\sqrt{2 \pi}}e^{-\frac{(z-\mu)^2}{2\sigma^2}}
-    \left[
-        1+
-        \frac{1}{6} \gamma_1 z \left(z^2-3\right)
-        + \frac{1}{24} \gamma_2 \left(z^4-6z^2+3\right)
-    \right],
-\end{align}
-где \(\Phi(z)\)\nbsp{}--- функция распределения (ФР) нормального закона,
-\(\phi\)\nbsp{}--- ФПР нормального закона, \(\gamma_1\)\nbsp{}--- асимметрия,
-\(\gamma_2\)\nbsp{}--- эксцесс, \(f\)\nbsp{}--- ФПР, \(F\)\nbsp{}--- функция
-распределения (ФР). Согласно\nbsp{}cite:rozhkov1990stochastic для аппликат
-морских волн значение асимметрии выбирается на интервале
-\(0,1\leq\gamma_1\leq{0,52}]\), а значение эксцесса на интервале
-\(0,1\leq\gamma_2\leq{0,7}\). Семейство плотностей распределения при различных
-параметрах показано на рис.\nbsp{}[[fig-skew-normal-1]].
+что соответствует формуле линейной теории для конечной глубины.
 
-#+name: fig-skew-normal-1
-#+begin_src R :file build/skew-normal-1-ru.pdf
-source(file.path("R", "common.R"))
-x <- seq(-3, 3, length.out=100)
-params <- data.frame(
-  skewness = c(0.00, 0.52, 0.00, 0.52),
-  kurtosis = c(0.00, 0.00, 0.70, 0.70),
-  linetypes = c("solid", "dashed", "dotdash", "dotted")
-)
-arma.skew_normal_1_plot(x, params)
-legend(
-  "topleft",
-  mapply(
-    function (s, k) {
-      as.expression(bquote(list(
-        gamma[1] == .(arma.fmt(s, 2)),
-        gamma[2] == .(arma.fmt(k, 2))
-      )))
-    },
-    params$skewness,
-    params$kurtosis
-  ),
-  lty = paste(params$linetypes)
-)
-#+end_src
+Различные записи решения уравнения Лапласа, в которых затухающая экспонента
+может встречаться как со знаком "+", так и со знаком "-", могут стать причиной
+разницы между формулами линейно теории и формулами, выведенными в данной работе,
+где вместо \(\sinh\) используется \(\cosh\). Выражение
+\(\frac{\Sinh{2\pi{k}(z+h)}}{\Sinh{2\pi{k}{h}}}\approx\frac{\sinh(2\pi{k}(z+h))}{\sinh(2\pi{k}{h})}\)
+превращается в строгое равенство на поверхности, и разница между правой левой
+частью увеличивается при приближении к дну водоема (для достаточно большой
+глубины ошибка вблизи поверхности жидкости незначительна). Поэтому для
+достаточно большой глубины можно использовать любую из функций (\(\cosh\) или
+\(\sinh\)) для вычисления потенциала скорости вблизи взволнованной поверхности.
 
-#+caption[График плотности распределения на основе РГШ]:
-#+caption: График плотности распределения закона, основанного на РГШ,
-#+caption: при различных значениях асимметрии \(\gamma_1\) и эксцесса
-#+caption: \(\gamma_2\).
-#+name: fig-skew-normal-1
-#+RESULTS: fig-skew-normal-1
-[[file:build/skew-normal-1-ru.pdf]]
+Сведение формул\nbsp{}eqref:eq-solution-2d и\nbsp{}eqref:eq-solution-2d-full к
+формулам линейной теории волн показывает, что формула\nbsp{}eqref:eq-solution-2d
+для жидкости бесконечной глубины не подходит для вычисления потенциала скорости
+с использованием метода Фурье, т.к.\nbsp{}не обладает необходимой для
+преобразования Фурье симметрией. Однако, для такого случая можно использовать
+формулу для конечной глубины, полагая \(h\) равным некоторому характерному
+значению глубины. Для стоячих волн сведение к формулам линейной теории
+происходит с аналогичными предположениями.
 
-**** Асимметричное нормальное распределение.
-Альтернативной аппроксимацией распределения волновых аппликат служит формула
-асимметричного нормального распределения (АНР):
+*** Трехмерное поле скоростей
+В трех измерениях исходная система уравнений\nbsp{}eqref:eq-problem
+переписывается как
 \begin{align}
-    \label{eq-skew-normal-2}
-    F(z; \alpha) & = \frac{1}{2}
-   \mathrm{erfc}\left[-\frac{z}{\sqrt{2}}\right]-2 T(z,\alpha ), \nonumber \\
-    f(z; \alpha) & = \frac{e^{-\frac{z^2}{2}}}{\sqrt{2 \pi }}
-   \mathrm{erfc}\left[-\frac{\alpha z}{\sqrt{2}}\right],
+    \label{eq-problem-3d}
+    & \phi_{xx} + \phi_{yy} + \phi_{zz} = 0,\\
+    & \zeta_t + \zeta_x\phi_x + \zeta_y\phi_y
+    =
+    \frac{\zeta_x}{\SqrtZeta{1 + \zeta_x^2 + \zeta_y^2}} \phi_x
+    +\frac{\zeta_y}{\SqrtZeta{1 + \zeta_x^2 + \zeta_y^2}} \phi_y
+    - \phi_z, & \text{на }z=\zeta(x,y,t).\nonumber
 \end{align}
-где \(T\)\nbsp{}--- функция Оуэна\nbsp{}cite:owen1956tables. Эта формула не
-позволяет задать значения асимметрии и эксцесса по отдельности\nbsp{}--- оба
-значения регулируются параметром \(\alpha\). Преимущество данной формулы лишь в
-относительной простоте вычисления: эта функция встроена в некоторые программы и
-библиотеки математических функций. График функции для разных значений \(\alpha\)
-представлен на рис.\nbsp{}[[fig-skew-normal-2]].
-
-#+name: fig-skew-normal-2
-#+begin_src R :file build/skew-normal-2-ru.pdf
-source(file.path("R", "common.R"))
-x <- seq(-3, 3, length.out=100)
-alpha <- c(0.00, 0.87, 2.25, 4.90)
-params <- data.frame(
-  alpha = alpha,
-  skewness = arma.bits.skewness_2(alpha),
-  kurtosis = arma.bits.kurtosis_2(alpha),
-  linetypes = c("solid", "dashed", "dotdash", "dotted")
-)
-arma.skew_normal_2_plot(x, params)
-legend(
-  "topleft",
-  mapply(
-    function (a, s, k) {
-      as.expression(bquote(list(
-        alpha == .(arma.fmt(a, 2)),
-        gamma[1] == .(arma.fmt(s, 2)),
-        gamma[2] == .(arma.fmt(k, 2))
-      )))
-    },
-    params$alpha,
-    params$skewness,
-    params$kurtosis
-  ),
-  lty = paste(params$linetypes)
-)
-#+end_src
-
-#+caption[График плотности АНР]:
-#+caption: График плотности АНР при различных значениях коэффициента
-#+caption: асимметрии \(\alpha\).
-#+name: fig-skew-normal-2
-#+RESULTS: fig-skew-normal-2
-[[file:build/skew-normal-2.pdf]]
-
-**** Тестирование.                                              :noexport:
-Решение уравнения\nbsp{}eqref:eq-distribution-transformation с выбранной
-функцией распределения можно произвести либо в каждой точке генерируемой
-поверхности, что даст наиболее точные результаты, либо в каждой точке
-фиксированной сетки, интерполировав решение методом наименьших квадратов (МНК).
-Во втором случае точность будет меньше. Например, интерполяция многочленом 12-го
-порядка на сетке из 500 узлов, построенной на промежутке
-\(-5\sigma_z\leq{z}\leq{5}\sigma_z\), дает погрешность
-\(\approx{0,43}\cdot10^{-3}\). Увеличение порядка многочлена приводит либо к
-переполнениям при интерполяции МНК, либо к дополнительным коэффициентам близким
-к нулю; увеличение размера сетки влияет на результат незначительно. В
-большинстве случаев трех коэффициентов ряда Грама---Шарлье было достаточно для
-преобразования АКФ; относительная погрешность без интерполяции составляет
-\(10^{-5}\).
-
-*** Алгоритм генерации белого шума
-Чтобы исключить периодичность из генерируемой моделью морского волнения
-реализации взволнованной поверхности, для генерации белого шума необходимо
-использовать ГПСЧ с достаточно большим периодом. В качестве такого генератора в
-работе используется параллельная реализация вихря
-Мерсенна\nbsp{}cite:matsumoto1998mersenne с периодом \(2^{19937}-1\). Это
-позволяет создавать апериодические реализации взволнованной морской поверхности
-для любых сценариев применения, встречаемых на практике.
-
-Запуск нескольких ГПСЧ с разными начальными состояниями в параллельных потоках
-не гарантирует некоррелированность генерируемых последовательностей
-псевдослучайных чисел, и для предоставления такой гарантии используется алгоритм
-динамического создания вихрей Мерсенна\nbsp{}cite:matsumoto1998dynamic. Суть
-алгоритма заключается в поиске таких матриц начальных состояний генераторов,
-которые бы дали максимально некоррелированные последовательности псевдослучайных
-чисел при параллельном запуске нескольких вихрей Мерсенна с этими начальными
-состояниями. Поскольку на поиск начальных состояний тратится значительное
-количество процессорного времени, то вектор состояний создается предварительно
-для заведомо большего количества параллельных потоков и сохраняется в файл,
-который впоследствии считывается основной программой перед началом генерации
-белого шума.
-
-*** Алгоритм генерации взволнованной поверхности
-В модели АРСС значение подъема взволнованной поверхности в каждой точке зависит
-от предыдущих по пространству и времени значений, из-за чего в начале реализации
-образуется так называемый /интервал разгона/
-(см.\nbsp{}рис.\nbsp{}[[fig-ramp-up-interval]])\nbsp{}--- промежуток, на котором
-реализация не соответствует заданной АКФ. Способ решения этой проблемы зависит
-от контекста, в котором производится имитационное моделирование.
-
-Если реализация используется в контексте расчета остойчивости судна без учета
-маневрирования, то интервал никак не повлияет результаты эксперимента, поскольку
-находится на границе (далеко от исследуемого морского объекта). Если изучается
-остойчивость судна в условиях маневрирования, то интервал проще всего исключить
-из реализации (размер интервала примерно равен числу коэффициентов АР по каждому
-из измерений). Однако, это приводит к потере большого числа точек, поскольку
-исключение происходит по каждому из трех измерений. Альтернативным подходом
-является генерация взволнованной поверхности на интервале разгона моделью ЛХ и
-генерация остальной реализации с помощью модели АРСС.
-
-В алгоритме генерации взволнованной поверхности используется параллелизм по
-данным: реализация делится на равные части, каждая из которых генерируется
-независимо. В случае модели СС в начале каждой из частей также присутствует
-интервал разгона. Для его исключения используется метод /сшивания/, часто
-применяемый в обработке цифровых
-сигналов\nbsp{}cite:oppenheim1989discrete,svoboda2011efficient,pavel2013algorithms.
-Суть метода заключается в добавлении интервала равного по размеру интервалу
-разгона в конец каждой из частей. Затем взволнованная поверхность генерируется в
-каждой точке каждой из частей (включая добавленный интервал), интервал в конце
-части \(N\) накладывается на интервал разгона в начале части \(N+1\), и значения
-в соответствующих точках складываются. В случае модели АР части генерируются с
-учетом авторегрессионных зависимостей между ними. Подробное описание
-параллельных и распределенных алгоритмов для каждой из моделей дано в
-разделе\nbsp{}[[#sec-parallel]] и\nbsp{}[[#sec-distributed]].
+Для ее решения также воспользуемся методом Фурье. Возьмем преобразование Фурье
+от обоих частей уравнений Лапласа и получим
+\begin{equation*}
+    -4 \pi^2 \left( u^2 + v^2 + w^2 \right)
+    \FourierY{\phi(x,y,z)}{u,v,w} = 0,
+\end{equation*}
+откуда имеем \(w=\pm{i}\sqrt{u^2+v^2}\). Решение уравнения будем искать в виде
+обратного преобразования Фурье
+\(\phi(x,y,z)=\InverseFourierY{E(u,v,w)}{x,y,z}\). Подставляя
+\(w=i\sqrt{u^2+v^2}=i\Kveclen\) в исходную формулу, получаем
+\begin{equation*}
+    \phi(x,y,z) = \InverseFourierY{
+        \left(
+            C_1 e^{2\pi \Kveclen z}
+            -C_2 e^{-2\pi \Kveclen z}
+        \right)
+        E(u,v)
+    }{x,y}.
+\end{equation*}
+Подставляя \(\phi\) в условие на дне водоема аналогично двухмерному случаю,
+получаем
+\begin{equation}
+    \label{eq-guessed-sol-3d}
+    \phi(x,y,z) = \InverseFourierY{
+        \Sinh{2\pi \Kveclen (z+h)} E(u,v)
+    }{x,y}.
+\end{equation}
+Подставляя выражение для \(\phi\) в граничное условие, получим
+\begin{equation*}
+    \arraycolsep=1.4pt
+    \begin{array}{rl}
+        \zeta_t = & i f_1(x,y) \InverseFourierY{2 \pi u \Sinh{2\pi \Kveclen (z+h)}E(u,v)}{x,y} \\
+        + & i f_2(x,y) \InverseFourierY{2 \pi v \Sinh{2\pi \Kveclen (z+h)}E(u,v)}{x,y} \\
+        - & \InverseFourierY{2 \pi \sqrt{u^2+v^2} \Sinh{2\pi \Kveclen (z+h)}E(u,v)}{x,y}
+    \end{array}
+\end{equation*}
+где \(f_1(x,y)={\zeta_x}/{\SqrtZeta{1+\zeta_x^2+\zeta_y^2}}-\zeta_x\) и
+\(f_2(x,y)={\zeta_y}/{\SqrtZeta{1+\zeta_x^2+\zeta_y^2}}-\zeta_y\).
 
-#+name: fig-ramp-up-interval
-#+begin_src R :file build/ramp-up-interval-ru.pdf
-source(file.path("R", "common.R"))
-arma.plot_ramp_up_interval(label="Интервал разгона")
-#+end_src
+Также как и в разделе\nbsp{}[[#sec-pressure-2d]] мы предполагаем, что
+\(\Sinh{2\pi{u}(z+h)}\approx\SinhX{2\pi{u}(z+h)}\) вблизи свободной поверхности,
+однако в трехмерном случае этого недостаточно для решения задачи. Для того чтобы
+получить явную формулу для коэффициентов \(E\), мы должны предположить, что
+преобразования Фурье в равенстве имеют радиально симметричные ядра,
+т.е.\nbsp{}заменить \(u\) и \(v\) на \(\Kveclen\). Есть два момента,
+поддерживающих это предположение. Во-первых, в численной реализации
+интегрирование ведется по положительным волновым числам, так что знак \(u\) и
+\(v\) не влияет на решение. Во-вторых, скорость роста \(\cosh\) в ядре интеграла
+значительно выше, чем скорость роста \(u\) или \(\Kveclen\), так что замена
+слабо влияет на величину решения. Несмотря на эти два момента, использование
+более математически строго подхода было бы предпочтительнее.
 
-#+caption[Интервал разгона в начале реализации]:
-#+caption: Интервал разгона в начале реализации.
-#+name: fig-ramp-up-interval
-#+RESULTS: fig-ramp-up-interval
-[[file:build/ramp-up-interval-ru.pdf]]
+Выполняя замену, применяя преобразование Фурье к обеим частям равенства, и,
+подставляя результат в\nbsp{}eqref:eq-guessed-sol-3d, получаем выражение для
+\(\phi\):
+\begin{equation*}
+    \label{eq-phi-3d}
+    \phi(x,y,z,t) = \InverseFourierY{
+        \frac{ \Sinh{\smash{2\pi \Kveclen (z+h)}} }{ 2\pi\Kveclen }
+        \frac{ \FourierY{ \zeta_t / \left( i f_1(x,y) + i f_2(x,y) - 1 \right)}{u,v} }
+        { \FourierY{\mathcal{D}_3\left( x,y,\zeta\left(x,y\right) \right)}{u,v} }
+    }{x,y},
+\end{equation*}
+где
+\(\FourierY{\mathcal{D}_3\left(x,y,z\right)}{u,v}=\Sinh{\smash{2\pi\Kveclen{}z}}\).
 
 *** Формулы нормировки для потенциалов скоростей
 :PROPERTIES:
@@ -1491,130 +1597,13 @@ arma.plot_ramp_up_interval(label="Интервал разгона")
 | \(\Fun{z}\)       | \(\delta (x+i z)\)                                           | \(\frac{1}{2 h}\mathrm{sech}\left(\frac{\pi  (x-i (h+z))}{2 h}\right)\)                                                                |
 | \(\FunSecond{z}\) | \(\frac{1}{2}\left[\delta (x-i z) + \delta (x+i z) \right]\) | \(\frac{1}{4 h}\left[\text{sech}\left(\frac{\pi  (x-i (h+z))}{2 h}\right)+\text{sech}\left(\frac{\pi  (x+i(h+z))}{2 h}\right)\right]\) |
 
-** Верификация модели АРСС
-:PROPERTIES:
-:CUSTOM_ID: sec-verification
-:END:
-
-Для модели АР в
-работах\nbsp{}cite:degtyarev2011modelling,degtyarev2013synoptic,boukhanovsky1997thesis
-экспериментальным путем были верифицированы
-- распределения различных характеристик волн (высоты волн, длины волн, длины
-  гребней, период волн, уклон волн, показатель трехмерности),
-- дисперсионное соотношение,
-- сохранение интегральных характеристик для случая смешанного волнения.
-В данной работе верифицируются как модель АР, так и СС путем сравнения
-распределений различных характеристик волн.
-
-*** Верификация интегральных характеристик взволнованной поверхности
-В\nbsp{}cite:rozhkov1990stochastic авторы показывают, что некоторые
-характеристики морских волн (перечисленные в табл.\nbsp{}[[tab-weibull-shape]])
-имеют распределение Вейбулла, а подъем взволнованной поверхности\nbsp{}---
-нормальное распределение. Для верификации генерируемых моделями АР и СС
-реализаций используются спрямленные диаграммы (графики, в которых по оси \(OX\)
-откладываются квантили функции распределения, вычисленные аналитически, а по оси
-\(OY\)\nbsp{}--- вычисленные экспериментально). Если экспериментально полученное
-распределение соответствует аналитическому, то график представляет собой прямую
-линию. Концы графика могут отклоняться от прямой линии, поскольку не могут быть
-надежно получены из реализации конечной длины.
-
-#+name: tab-weibull-shape
-#+caption[Значение коэффициента формы распределения Вейбулла]:
-#+caption: Значение коэффициента формы распределения Вейбулла для различных
-#+caption: характеристик волн.
-#+attr_latex: :booktabs t
-| Характеристика          | Коэффициент формы \(k\) |
-|-------------------------+-------------------------|
-|                         | <l>                     |
-| Высота волны            | 2                       |
-| Длина волны             | 2,3                     |
-| Длина гребня волны      | 2,3                     |
-| Период волны            | 3                       |
-| Уклон волны             | 2,5                     |
-| Показатель трехмерности | 2,5                     |
-
-Верификация производится для стоячих и прогрессивных волн. Соответствующие АКФ и
-спрямленные диаграммы распределений характеристик волн представлены на
-рис.\nbsp{}[[acf-slices]],\nbsp{}[[standing-wave-distributions]],\nbsp{}[[propagating-wave-distributions]].
-
-#+name: propagating-wave-distributions
-#+begin_src R :file build/propagating-wave-qqplots-ru.pdf
-source(file.path("R", "common.R"))
-par(pty="s", mfrow=c(2, 2))
-arma.qqplot_grid(
-  file.path("build", "propagating_wave"),
-  c("elevation", "heights_y", "lengths_y", "periods"),
-  c("подъем", "высота по Y", "длина по Y", "период"),
-  xlab="x",
-  ylab="y"
-)
-#+end_src
-
-#+caption[Спрямленные диаграммы для прогрессивных волн]:
-#+caption: Спрямленные диаграммы для прогрессивных волн.
-#+name: propagating-wave-distributions
-#+RESULTS: propagating-wave-distributions
-[[file:build/propagating-wave-qqplots.pdf]]
-
-#+name: standing-wave-distributions
-#+begin_src R :file build/standing-wave-qqplots-ru.pdf
-source(file.path("R", "common.R"))
-par(pty="s", mfrow=c(2, 2))
-arma.qqplot_grid(
-  file.path("build", "standing_wave"),
-  c("elevation", "heights_y", "lengths_y", "periods"),
-  c("подъем", "высота по Y", "длина по Y", "период"),
-  xlab="x",
-  ylab="y"
-)
-#+end_src
-
-#+caption[Спрямленные диаграммы для стоячих волн]:
-#+caption: Спрямленные диаграммы для стоячих волн.
-#+name: standing-wave-distributions
-#+RESULTS: standing-wave-distributions
-[[file:build/standing-wave-qqplots-ru.pdf]]
-
-#+name: acf-slices
-#+header: :width 6 :height 9
-#+begin_src R :file build/acf-slices-ru.pdf
-source(file.path("R", "common.R"))
-propagating_acf <- read.csv(file.path("build", "propagating_wave", "acf.csv"))
-standing_acf <- read.csv(file.path("build", "standing_wave", "acf.csv"))
-par(mfrow=c(5, 2), mar=c(0,0,0,0))
-for (i in seq(0, 4)) {
-  arma.wavy_plot(standing_acf, i, zlim=c(-5,5))
-  arma.wavy_plot(propagating_acf, i, zlim=c(-5,5))
-}
-#+end_src
-
-#+caption[Временные срезы АКФ для стоячих и прогрессивных волн]:
-#+caption: Временные срезы АКФ для стоячих (слева) и прогрессивных (справа) волн.
-#+name: acf-slices
-#+RESULTS: acf-slices
-[[file:build/acf-slices-ru.pdf]]
-
-Хвосты распределений на рис.\nbsp{}[[propagating-wave-distributions]] отклоняются от
-оригинального распределения для характеристик отдельных волн, поскольку каждую
-волну необходимо извлечь из полученной взволнованной поверхности, чтобы измерить
-ее длину, период и высоту. Алгоритм, который бы гарантировал безошибочное
-извлечение всех волн, не известен, поскольку волны могут и часто накладываются
-друг на друга. Правый хвост распределения Вейбулла отклоняется больше, поскольку
-он представляет редко возникающие волны.
-
-Степень соответствия для стоячих волн (рис.\nbsp{}[[standing-wave-distributions]])
-ниже для высот и длин, примерно одинакова для подъема поверхности и выше для
-периодов волн. Более низкая степень соответствия длин и высот может быть
-результатом того, что распределения были получены эмпирически для морских волн,
-которые, в основном, являются прогрессивными, и аналогичные распределения для
-стоячих волн могут отличаться. Более высокая степень соответствия периодов волн
-является следствием того, что периоды стоячих волн извлекаются более точно,
-поскольку волн не перемещаются вне моделируемой области взволнованной
-поверхности. Одинаковая степень соответствия для подъема поверхности получается
-из-за того, что это характеристика поверхности (и соответствующего процесса АР
-или СС), и она не зависит от типа волн.
+*** Заключение
+Полученные в данном разделе формулы позволяют произвести численный расчет поля
+скоростей (а значит и поля давлений) вблизи дискретно или математически заданной
+взволнованной морской поверхности, минуя предположения линейной теории и теории
+волн малых амплитуд.
 
-*** Верификация полей потенциалов скоростей
+** Верификация полей потенциалов скоростей
 :PROPERTIES:
 :CUSTOM_ID: sec-compare-formulae
 :END:
@@ -1627,7 +1616,7 @@ for (i in seq(0, 4)) {
 ввиду выводы раздела\nbsp{}[[#sec-pressure-2d]], сравниваются только формулы для
 случая конечной глубины.
 
-**** Отличие от формул линейной теории волн.
+*** Отличие от формул линейной теории волн.
 Для того чтобы получить поля потенциалов скоростей, взволнованная морская
 поверхность генерировалась с помощью модели АР с варьированием амплитуды волн. В
 численной реализации волновые числа в преобразованиях Фурье выбирались на
@@ -1700,7 +1689,7 @@ arma.plot_velocity_potential_field_legend(
 [[file:build/plain-wave-velocity-field-comparison-ru.pdf]]
 
 
-**** Отличие от формул теории волн малой амплитуды.
+*** Отличие от формул теории волн малой амплитуды.
 Эксперимент, в котором сравнивались поля потенциалов скоростей, полученные
 численно различными формулами, показал, что поля скоростей, полученные по
 формуле\nbsp{}eqref:eq-solution-2d-full и формуле для волн малой
@@ -1746,86 +1735,6 @@ arma.plot_velocity(
 #+RESULTS: fig-velocity-field-2d
 [[file:build/large-and-small-amplitude-velocity-field-comparison-ru.pdf]]
 
-*** Верификация нелинейного безынерционного преобразования
-Для того чтобы измерить влияние НБП на форму результирующей взволнованной
-поверхности, было сгенерировано три реализации:
-- реализация с Гауссовым распределением (без НБП),
-- реализация с распределением на основе ряда Грама---Шарлье (РГШ), и
-- реализация с асимметричным нормальным распределением (АНР).
-Начальные состояния ГПСЧ были заданы одинаковыми для всех запусков программы,
-чтобы модель АРСС выдавала одни и те же значения для каждой реализации. Было
-проведено два эксперимента: для стоячих и прогрессивных волн с АКФ, заданными
-формулами из раздела\nbsp{}[[#sec-wave-acfs]].
-
-В то время как эксперимент показал, что применение НБП с распределением РГШ
-увеличивает крутизну волн, то же самое нельзя сказать об асимметричном
-нормальном распределении (рис.\nbsp{}[[fig-nit]]). Использование этого распределения
-приводит к взволнованной поверхности, в которой аппликаты всегда больше или
-равны нулю. Таким образом, асимметричное нормальное распределение не подходит
-для НБП. НБП увеличивает высоту и крутизну как прогрессивных, так и стоячих
-волн. Увеличение асимметрии или эксцесса РГШ приводит в увеличению как высоты,
-так и крутизны волн. Ошибка аппроксимации АКФ (ур.\nbsp{}eqref:eq-nit-error)
-принимает значения от 0,20 для РГШ до 0,70 для АНР (табл.\nbsp{}[[tab-nit-error]]).
-
-#+name: fig-nit
-#+header: :width 7 :height 7
-#+begin_src R :file build/nit.pdf
-source(file.path("R", "nonlinear.R"))
-par(mfrow=c(2, 1), mar=c(4,4,4,0.5), family='serif')
-args <- list(
-  graphs=c('Гауссово', 'РГШ', 'АНР'),
-  linetypes=c('solid', 'dashed', 'dotted')
-)
-args$title <- 'Прогрессивные волны'
-arma.plot_nonlinear(file.path("build", "nit-propagating"), args)
-args$title <- 'Стоячие волны'
-arma.plot_nonlinear(file.path("build", "nit-standing"), args)
-#+end_src
-
-#+name: fig-nit
-#+caption[Срезы взволнованной поверхности с различными распределениями]:
-#+caption: Срезы взволнованной поверхности с различными распределениями
-#+caption: волновых аппликат (Гауссово, РГШ и асимметричное нормальное).
-#+RESULTS: fig-nit
-[[file:build/nit.pdf]]
-
-#+name: tab-nit-error
-#+caption[Ошибки аппроксимации АКФ]:
-#+caption: Ошибки аппроксимации АКФ (разность дисперсий) для различных
-#+caption: распределений волновых аппликат. \(N\)\nbsp{}--- количество
-#+caption: коэффициентов аппроксимации АКФ.
-#+attr_latex: :booktabs t
-| Тип волн      | Распред. | \(\gamma_1\) | \(\gamma_2\) | \(\alpha\) | Ошибка | \(N\) | Высота волн |
-|---------------+----------+--------------+--------------+------------+--------+-------+-------------|
-|               |          |          <r> |          <r> |        <r> |    <r> |       |         <r> |
-| прогрессивные | Гауссово |              |              |            |        |       |        2,41 |
-| прогрессивные | РГШ      |         2,25 |          0,4 |            |   0,20 |     2 |        2,75 |
-| прогрессивные | АНР      |              |              |          1 |   0,70 |     3 |        1,37 |
-| стоячие       | Гауссово |              |              |            |        |       |        1,73 |
-| стоячие       | РГШ      |         2,25 |          0,4 |            |   0,26 |     2 |        1,96 |
-| стоячие       | АНР      |              |              |          1 |   0,70 |     3 |        0,94 |
-
-Таким образом, единственный тестовый сценарий, который показал приемлемые
-результаты\nbsp{}--- это реализации с распределением на основе РГШ для
-прогрессивных и стоячих волн. АНР искажает взволнованную поверхность для обоих
-типов волн. Реализации с распределением на основе РГШ характеризуются большой
-ошибкой аппроксимации АКФ, что приводит к увеличению высоты волн. Причина
-большой ошибки заключается в неточности аппроксимации РГШ, которая не сходится
-для всевозможных функций\nbsp{}cite:wallace1958asymptotic. Несмотря на большую
-ошибку, изменение высоты волн невелико (табл.\nbsp{}[[tab-nit-error]]).
-
-*** Нефизическая природа модели
-В силу своей нефизической природы модель АРСС не включает в себя понятие морской
-волны; вместо этого она моделирует взволнованную поверхность как единое целое.
-Движения отдельных волн и их форма часто получаются грубыми, а точное количество
-генерируемых волн неизвестно. Несмотря на это, интегральные характеристики
-взволнованной поверхности соответствуют реальным морским волнам.
-
-Теоретически, профили самих морских волн могут быть использованы в качестве АКФ,
-если предварительно обеспечить их экспоненциальное затухание. Это может
-позволить генерировать волны произвольных профилей и является одной из тем
-дальнейших исследований.
-
 * Высокопроизводительный программный комплекс для моделирования морского волнения
 ** Реализация для систем с общей памятью (SMP)
 **** Параллельные алгоритмы для моделей АР, СС и ЛХ.
@@ -1913,9 +1822,9 @@ arma.plot_ar_cubes_2d(3, 3, xlabel="Индекс части (X)", ylabel="Инд
 вычисляется от результата. После этого, каждая часть записывается в выходной
 массив, а перекрывающие друг друга (из-за заполнения нулями) точки складываются
 друг с другом. Этот алгоритм известен в области обработки сигналов как
-"overlap-add"\nbsp{}cite:svoboda2011efficient. Заполнение нулями необходимо для
-предотвращения маскированных ошибок: без него результатом вычислений была бы
-циклическая свертка.
+"overlap-add"\nbsp{}cite:oppenheim1989discrete,svoboda2011efficient,pavel2013algorithms.
+Заполнение нулями необходимо для предотвращения маскированных ошибок: без него
+результатом вычислений была бы циклическая свертка.
 
 Несмотря на то что алгоритм модели СС разбивает поверхность на те же самые части
 (но, возможно, другого размера), что и алгоритм модели АР, отсутствие
@@ -2375,6 +2284,57 @@ OpenGL увеличивает производительность путем и
 артефактов изображения на экране, которые можно убрать, только перезагрузив
 компьютер.
 
+**** Параллельный алгоритм генерации белого шума
+Чтобы исключить периодичность из генерируемой моделью морского волнения
+реализации взволнованной поверхности, для генерации белого шума необходимо
+использовать ГПСЧ с достаточно большим периодом. В качестве такого генератора в
+работе используется параллельная реализация вихря
+Мерсенна\nbsp{}cite:matsumoto1998mersenne с периодом \(2^{19937}-1\). Это
+позволяет создавать апериодические реализации взволнованной морской поверхности
+для любых сценариев применения, встречаемых на практике.
+
+Запуск нескольких ГПСЧ с разными начальными состояниями в параллельных потоках
+не гарантирует некоррелированность генерируемых последовательностей
+псевдослучайных чисел, и для предоставления такой гарантии используется алгоритм
+динамического создания вихрей Мерсенна\nbsp{}cite:matsumoto1998dynamic. Суть
+алгоритма заключается в поиске таких матриц начальных состояний генераторов,
+которые бы дали максимально некоррелированные последовательности псевдослучайных
+чисел при параллельном запуске нескольких вихрей Мерсенна с этими начальными
+состояниями. Поскольку на поиск начальных состояний тратится значительное
+количество процессорного времени, то вектор состояний создается предварительно
+для заведомо большего количества параллельных потоков и сохраняется в файл,
+который впоследствии считывается основной программой перед началом генерации
+белого шума.
+
+**** Устранение интервала разгона.
+В модели АР значение подъема взволнованной поверхности в каждой точке зависит от
+предыдущих по пространству и времени значений, из-за чего в начале реализации
+образуется так называемый /интервал разгона/
+(см.\nbsp{}рис.\nbsp{}[[fig-ramp-up-interval]])\nbsp{}--- промежуток, на котором
+реализация не соответствует заданной АКФ. Способ решения этой проблемы зависит
+от контекста, в котором производится имитационное моделирование. Если реализация
+используется в контексте расчета остойчивости судна без учета маневрирования, то
+интервал никак не повлияет результаты эксперимента, поскольку находится на
+границе (далеко от исследуемого морского объекта). Альтернативным подходом
+является генерация взволнованной поверхности на интервале разгона моделью ЛХ и
+генерация остальной реализации с помощью модели АР. Если изучается остойчивость
+судна в условиях маневрирования, то интервал проще всего исключить из реализации
+(размер интервала примерно равен числу коэффициентов АР по каждому из
+измерений). Однако, это приводит к потере большого числа точек, поскольку
+исключение происходит по каждому из трех измерений.
+
+#+name: fig-ramp-up-interval
+#+begin_src R :file build/ramp-up-interval-ru.pdf
+source(file.path("R", "common.R"))
+arma.plot_ramp_up_interval(label="Интервал разгона")
+#+end_src
+
+#+caption[Интервал разгона в начале реализации процесса АР]:
+#+caption: Интервал разгона в начале реализации процесса АР.
+#+name: fig-ramp-up-interval
+#+RESULTS: fig-ramp-up-interval
+[[file:build/ramp-up-interval-ru.pdf]]
+
 **** Заключение.
 Тесты показали, что видеокарта превосходит центральный процессор по
 производительности в задачах с большим количеством арифметических операций,
diff --git a/arma-thesis.org b/arma-thesis.org
@@ -598,486 +598,195 @@ propagating waves. Using mixed ARMA process for both types of waves might
 increase model precision, given that coefficient correction formulae for three
 dimensions become available, which is the objective of the future research.
 
-** Known pressure field determination formulae
-**** Small amplitude waves theory.
-In\nbsp{}cite:stab2012,degtyarev1998modelling,degtyarev1997analysis the authors
-propose a solution for inverse problem of hydrodynamics of potential flow in the
-framework of small-amplitude wave theory (under assumption that wave length is
-much larger than height: \(\lambda \gg h\)). In that case inverse problem is
-linear and reduces to Laplace equation with mixed boundary conditions, and
-equation of motion is solely used to determine pressures for calculated velocity
-potential derivatives. The assumption of small amplitudes means the slow decay
-of wind wave coherence function, i.e. small change of local wave number in time
-and space compared to the wavy surface elevation (\(z\) coordinate). This
-assumption allows to calculate elevation \(z\) derivative as \(\zeta_z=k\zeta\),
-where \(k\) is wave number. In two-dimensional case the solution is written
-explicitly as
-\begin{align}
-    \left.\frac{\partial\phi}{\partial x}\right|_{x,t}= &
-        -\frac{1}{\sqrt{1+\alpha^{2}}}e^{-I(x)}
-            \int\limits_{0}^x\frac{\partial\dot{\zeta}/\partial
-                z+\alpha\dot{\alpha}}{\sqrt{1+\alpha^{2}}}e^{I(x)}dx,\label{eq-old-sol-2d}\\
-    I(x)= & \int\limits_{0}^x\frac{\partial\alpha/\partial z}{1+\alpha^{2}}dx,\nonumber
-\end{align}
-where \(\alpha\) is wave slope. In three-dimensional case the solution is
-written in the form of elliptic partial differential equation (PDE):
-\begin{align*}
-    & \frac{\partial^2 \phi}{\partial x^2} \left( 1 + \alpha_x^2 \right) +
-    \frac{\partial^2 \phi}{\partial y^2} \left( 1 + \alpha_y^2 \right) +
-    2\alpha_x\alpha_y \frac{\partial^2 \phi}{\partial x \partial y} + \\
-    & \left(
-        \frac{\partial \alpha_x}{\partial z} +
-        \alpha_x \frac{\partial \alpha_x}{\partial x} +
-        \alpha_y \frac{\partial \alpha_x}{\partial y}
-    \right) \frac{\partial \phi}{\partial x} + \\
-    & \left(
-        \frac{\partial \alpha_y}{\partial z} +
-        \alpha_x \frac{\partial \alpha_y}{\partial x} +
-        \alpha_y \frac{\partial \alpha_y}{\partial y}
-    \right) \frac{\partial \phi}{\partial y} + \\
-    & \frac{\partial \dot{\zeta}}{\partial z} +
-    \alpha_x \dot{\alpha_x} + \alpha_y \dot{\alpha_y} = 0.
-\end{align*}
-The authors suggest transforming this equation to finite differences and solve
-it numerically.
-
-As will be shown in section\nbsp{}[[#sec-compare-formulae]],
-formula\nbsp{}eqref:eq-old-sol-2d diverges when attempted to calculate velocity
-field for large-amplitude waves, and this is the reason that it can not be used
-in conjunction with a model, that generates arbitrary-amplitude waves.
-
-**** Linearisation of boundary condition.
-:PROPERTIES:
-:CUSTOM_ID: linearisation
-:END:
-
-LH model allows to derive an explicit formula for velocity field by linearising
-kinematic boundary condition. Velocity potential formula is written as
-\begin{equation*}
-\phi(x,y,z,t) = \sum_n \frac{c_n g}{\omega_n}
-     e^{\sqrt{u_n^2+v_n^2} z}
-     \sin(u_n x + v_n y - \omega_n t + \epsilon_n).
-\end{equation*}
-This formula is differentiated to obtain velocity potential derivatives, which
-are plugged to dynamic boundary condition to determine pressure field.
-
-** Determining wave pressures for discretely given wavy surface
-Analytic solutions to boundary problems in classical equations are often used to
-study different properties of the solution, and for that purpose general
-solution formula is too difficult to study, as it contains integrals of unknown
-functions. Fourier method is one of the methods to find analytic solutions to
-PDE. It is based on Fourier transform, applying of which to some PDEs reduces
-them to algebraic equations, and the solution is written as inverse Fourier
-transform of some function (which may contain Fourier transforms of other
-functions). Since it is not possible to write analytic forms of these Fourier
-transforms in some cases, unique solutions are found and their behaviour is
-studied in different domains instead. At the same time, computing discrete
-Fourier transforms numerically is possible for any discretely defined function
-using FFT algorithms. These algorithms use symmetry of complex exponentials to
-decrease asymptotic complexity from \(\mathcal{O}(n^2)\) to
-\(\mathcal{O}(n\log_{2}n)\). So, even if general solution contains Fourier
-transforms of unknown functions, they still can be computed numerically, and FFT
-family of algorithms makes this approach efficient.
-
-Alternative approach to solve PDEs is to reduce them to difference equations,
-which are solved by constructing various numerical schemes. This approach leads
-to approximate solution, and asymptotic complexity of corresponding algorithms
-is comparable to that of FFT. For example, for stationary elliptic PDE an
-implicit numerical scheme is constructed; the scheme is computed by iterative
-method on each step of which a tridiagonal or five-diagonal system of algebraic
-equations is solved by Thomas algorithm. Asymptotic complexity of this approach
-is \(\mathcal{O}({n}{m})\), where \(n\) is the number of wavy surface grid
-points and \(m\) is the number of iterations.
-
-Despite their wide spread, iterative algorithms are inefficient on parallel
-computer architectures due to inevitable process synchronisation after each
-iteration; in particular, their mapping to co-processors may involve copying
-data in and out of the co-processor on each iteration, which negatively affects
-their performance. At the same time, high number of Fourier transforms in the
-solution is an advantage, rather than a disadvantage. First, solutions obtained
-by Fourier method are explicit, hence their implementations scale with the large
-number of parallel computer cores. Second, there are implementations of FFT
-optimised for different processor architectures as well as co-processors (GPU,
-MIC) which makes it easy to get high performance on any computing platform.
-These advantages substantiate the choice of Fourier method to obtain explicit
-solution to the problem of determining pressures under wavy sea surface.
-
-*** Two-dimensional velocity field
-:PROPERTIES:
-:CUSTOM_ID: sec-pressure-2d
-:END:
+** Modelling non-linearity of sea waves
+ARMA model allows to model asymmetry of wave elevation distribution,
+i.e.\nbsp{}generate sea waves, distribution of \(z\)-coordinate of which has
+non-nought kurtosis and asymmetry. Such distribution is inherent to real sea
+waves\nbsp{}cite:longuet1963nonlinear and given by either polynomial
+approximation of /in situ/ data or analytic formula. Wave asymmetry is modelled
+by non-linear inertia-less transform (NIT) of stochastic process, however,
+transforming resulting wavy surface means transforming initial ACF. In order to
+alleviate this, ACF must be preliminary transformed as shown
+in\nbsp{}cite:boukhanovsky1997thesis.
 
-**** Formula for infinite depth fluid.
-Two-dimensional Laplace equation with Robin boundary condition is written as
-\begin{align}
-    \label{eq-problem-2d}
-    & \phi_{xx}+\phi_{zz}=0,\\
-    & \zeta_t + \zeta_x\phi_x = \frac{\zeta_x}{\sqrt{1 + \zeta_x^2}} \phi_x - \phi_z, & \text{на }z=\zeta(x,t).\nonumber
-\end{align}
-Use Fourier method to solve this problem. Applying Fourier transform to both
-sides of the equation yields
-\begin{equation*}
-    -4 \pi^2 \left( u^2 + v^2 \right)
-    \FourierY{\phi(x,z)}{u,v} = 0,
-\end{equation*}
-hence \(v = \pm i u\). Hereinafter we use the following symmetric form of
-Fourier transform:
-\begin{equation*}
-    \FourierY{f(x,y)}{u,v} =
-    \iint\limits_{-\infty}^{\phantom{--}\infty}
-    f(x,y)
-    e^{-2\pi i (x u + y v)}
-    dx dy.
-\end{equation*}
-We seek solution in the form of inverse Fourier transform
-\(\phi(x,z)=\InverseFourierY{E(u,v)}{x,z}\). Plugging[fn::\(v={-i}{u}\) is not
-applicable because velocity potential must go to nought when depth goes to
-infinity.] \(v={i}{u}\) into the formula yields
+**** Wavy surface transformation.
+Explicit formula \(z=f(y)\) that transforms wavy surface to desired
+one-dimensional distribution \(F(z)\) is the solution of non-linear
+transcendental equation \(F(z)=\Phi(y)\), where \(\Phi(y)\)\nbsp{}---
+one-dimensional Gaussian distribution. Since distribution of wave elevation is
+often given by some approximation based on /in situ/ data, this equation is
+solved numerically with respect to \(z_k\) in each grid point \(y_k|_{k=0}^N\)
+of generated wavy surface. In this case equation is rewritten as
 \begin{equation}
-    \label{eq-guessed-sol-2d}
-    \phi(x,z) = \InverseFourierY{e^{2\pi u z}E(u)}{x}.
+    \label{eq-distribution-transformation}
+    F(z_k)
+    =
+    \frac{1}{\sqrt{2\pi}}
+    \int\limits_0^{y_k} \exp\left[ -\frac{t^2}{2} \right] dt
+    .
 \end{equation}
-In order to make substitution \(z=\zeta(x,t)\) not interfere with Fourier
-transforms, we rewrite\nbsp{}eqref:eq-guessed-sol-2d as a convolution:
+Since, distribution functions are monotonic, the simplest interval halving
+(bisection) numerical method is used to solve this equation.
+
+**** Preliminary ACF transformation.
+In order to transform ACF \(\gamma_z\) of the process, it is expanded in series
+of Hermite polynomials (Gram---Charlier series)
 \begin{equation*}
-    \phi(x,z)
+    \gamma_z \left( \vec u \right)
     =
-    \Fun{z}
-    \ast
-    \InverseFourierY{E(u)}{x},
+    \sum\limits_{m=0}^{\infty}
+    C_m^2 \frac{\gamma_y^m \left( \vec u \right)}{m!},
 \end{equation*}
-where \(\Fun{z}\)\nbsp{}--- a function, form of which is defined in
-section\nbsp{}[[#sec-compute-delta]] and which satisfies equation
-\(\FourierY{\Fun{z}}{u}=e^{2\pi{u}{z}}\). Plugging formula \(\phi\) into the
-boundary condition yields
+where
 \begin{equation*}
-    \zeta_t
-    =
-    \left( i f(x) - 1 \right)
-    \left[
-        \Fun{z}
-        \ast
-        \InverseFourierY{2\pi u E(u)}{x}
-    \right],
+    C_m = \frac{1}{\sqrt{2\pi}}
+  \int\limits_{0}^\infty
+    f(y) H_m(y) \exp\left[ -\frac{y^2}{2} \right],
 \end{equation*}
-where \(f(x)={\zeta_x}/{\sqrt{1+\zeta_x^2}}-\zeta_x\). Applying Fourier transform
-to both sides of this equation yields formula for coefficients \(E\):
+\(H_m\)\nbsp{}--- Hermite polynomial, and \(f(y)\)\nbsp{}--- solution to
+equation\nbsp{}eqref:eq-distribution-transformation. Plugging polynomial
+approximation \(f(y)\approx\sum\limits_{i}d_{i}y^i\) and analytic formulae for
+Hermite polynomial yields
 \begin{equation*}
-    E(u) =
-    \frac{1}{2\pi u}
-    \frac{
-    \FourierY{\zeta_t / \left(i f(x) - 1\right)}{u}
-    }{
-    \FourierY{\Fun{z}}{u}
-    }
+    \frac{1}{\sqrt{2\pi}}
+    \int\limits_\infty^\infty
+    y^k \exp\left[ -\frac{y^2}{2} \right]
+    =
+    \begin{cases}
+        (k-1)!! & \text{if }k\text{ is even},\\
+        0       & \text{if }k\text{ is odd},
+    \end{cases}
 \end{equation*}
-Finally, substituting \(z\) for \(\zeta(x,t)\) and plugging resulting equation
-into\nbsp{}eqref:eq-guessed-sol-2d yields formula for \(\phi(x,z)\):
+which simplifies the former equation. Optimal number of coefficients \(C_m\) is
+determined by computing them sequentially and stopping when variances of both
+fields become equal with desired accuracy \(\epsilon\):
 \begin{equation}
-    \label{eq-solution-2d}
-    \boxed{
-        \phi(x,z)
-        =
-        \InverseFourierY{
-            \frac{e^{2\pi u z}}{2\pi u}
-            \frac{
-            \FourierY{ \zeta_t / \left(i f(x) - 1\right) }{u}
-            }{
-            \FourierY{ \Fun{\zeta(x,t)} }{u}
-            }
-        }{x}.
-    }
+    \label{eq-nit-error}
+    \left| \Var{z} - \sum\limits_{k=0}^m
+    \frac{C_k^2}{k!} \right| \leq \epsilon.
 \end{equation}
 
-Multiplier \(e^{2\pi{u}{z}}/(2\pi{u})\) makes graph of a function to which
-Fourier transform is applied asymmetric with respect to \(OY\) axis. This makes
-it difficult to use FFT, because it expects periodic function which takes nought
-values on both ends of the interval. At the same time, using numerical
-integration instead of FFT is not faster than solving the initial system of
-equations with numerical schemes. This problem is alleviated by using
-formula\nbsp{}eqref:eq-solution-2d-full for finite depth fluid with knowingly
-large depth \(h\). This formula is derived in the following section.
+In\nbsp{}cite:boukhanovsky1997thesis the author suggests using polynomial
+approximation \(f(y)\) also for wavy surface transformation, however, in
+practice sea surface realisation often contains points, in which
+\(z\)-coordinate is beyond the limits of the approximation, leading to sharp
+precision decrease. In these points it is more efficient to solve
+equation\nbsp{}eqref:eq-distribution-transformation by bisection method. Using
+the same approximation in Gram---Charlier series does not lead to such errors.
 
-**** Formula for finite depth fluid.
-On the sea bottom vertical fluid velocity component equals nought,
-i.e.\nbsp{}\(\phi_z=0\) on \(z=-h\), where \(h\)\nbsp{}--- water depth. In this
-case equation \(v=-{i}{u}\), which came from Laplace equation, can not be
-neglected, hence the solution is sought in the following form:
-\begin{equation}
-    \phi(x,z)
-    =
-    \InverseFourierY{
-        \left( C_1 e^{2\pi u z} + C_2 e^{-2\pi u z} \right)
-        E(u)
-    }{x}.
-    \label{eq-guessed-sol-2d-full}
-\end{equation}
-Plugging \(\phi\) into the boundary condition on the sea bottom yields
-\begin{equation*}
-    C_1 e^{-2\pi u h} - C_2 e^{2\pi u h} = 0,
-\end{equation*}
-hence \(C_1=\frac{1}{2}C{e}^{2\pi{u}{h}}\) and
-\(C_2=-\frac{1}{2}C{e}^{-2\pi{u}{h}}\). Constant \(C\) may take arbitrary values
-here, because after plugging it becomes part of unknown coefficients \(E(u)\).
-Plugging formulae for \(C_1\) and \(C_2\)
-into\nbsp{}eqref:eq-guessed-sol-2d-full yields
-\begin{equation*}
-    \phi(x,z) = \InverseFourierY{ \Sinh{2\pi u (z+h)} E(u) }{x}.
-\end{equation*}
-Plugging \(\phi\) into the boundary condition on the free surface yields
-\begin{equation*}
-    \zeta_t = f(x) \InverseFourierY{ 2\pi i u \Sinh{2\pi u (z+h)} E(u) }{x}
-            - \InverseFourierY{ 2\pi u \SinhX{2\pi u (z+h)} E(u) }{x}.
-\end{equation*}
-Here \(\sinh\) and \(\cosh\) give similar results near free surface, and since this
-is the main area of interest in practical applications, we assume that
-\(\Sinh{2\pi{u}(z+h)}\approx\SinhX{2\pi{u}(z+h)}\). Performing analogous to the
-previous section transformations yields final formula for \(\phi(x,z)\):
-\begin{equation}
-\boxed{
-    \phi(x,z,t)
-    =
-  \InverseFourierY{
-        \frac{\Sinh{2\pi u (z+h)}}{2\pi u}
-        \frac{
-            \FourierY{ \zeta_t / \left(i f(x) - 1\right) }{u}
-        }{
-            \FourierY{ \FunSecond{\zeta(x,t)} }{u}
-        }
-    }{x},
-}
-    \label{eq-solution-2d-full}
-\end{equation}
-where \(\FunSecond{z}\) is a function, a form of which is defined in
-section\nbsp{}[[#sec-compute-delta]] and which satisfies equation
-\(\FourierY{\FunSecond{z}}{u}=\Sinh{2\pi{u}{z}}\).
+**** Gram---Charlier series expansion.
+In\nbsp{}cite:huang1980experimental the authors experimentally show, that PDF of
+sea surface elevation is distinguished from normal distribution by non-nought
+kurtosis and skewness. In\nbsp{}cite:rozhkov1996theory the authors show, that
+this type of PDF expands in Gram---Charlier series (GCS):
+\begin{align}
+    \label{eq-skew-normal-1}
+    & F(z; \mu=0, \sigma=1, \gamma_1, \gamma_2) \approx
+    \Phi(z; \mu, \sigma) \frac{2 + \gamma_2}{2}
+    - \frac{2}{3} \phi(z; \mu, \sigma)
+    \left(\gamma_2 z^3+\gamma_1
+    \left(2 z^2+1\right)\right) \nonumber
+    \\
+    & f(z; \mu, \sigma, \gamma_1, \gamma_2) \approx
+    \phi(z; \mu, \sigma)
+    \left[
+        1+
+        \frac{1}{3!} \gamma_1 H_3 \left(\frac{z-\mu}{\sigma}\right)
+        + \frac{1}{4!} \gamma_2 H_4 \left(\frac{z-\mu}{\sigma}\right)
+    \right],
+\end{align}
+where \(\Phi(z)\)\nbsp{}--- cumulative density function (CDF) of normal
+distribution, \(\phi\)\nbsp{}--- PDF of normal distribution,
+\(\gamma_1\)\nbsp{}--- skewness, \(\gamma_2\)\nbsp{}--- kurtosis,
+\(f\)\nbsp{}--- PDF, \(F\)\nbsp{}--- cumulative distribution function (CDF).
+According to\nbsp{}cite:rozhkov1990stochastic for sea waves skewness is selected
+from interval \(0.1\leq\gamma_1\leq{0.52}]\) and kurtosis from interval
+\(0.1\leq\gamma_2\leq{0.7}\). Family of probability density functions for
+different parameters is shown in fig.\nbsp{}[[fig-skew-normal-1]].
 
-**** Reduction to the formulae from linear wave theory.
-Check the validity of derived formulae by substituting \(\zeta(x,t)\) with known
-analytic formula for plain waves. Symbolic computation of Fourier transforms in
-this section were performed in Mathematica\nbsp{}cite:mathematica10. In the
-framework of linear wave theory assume that waves have small amplitude compared
-to their lengths, which allows us to simplify initial system of
-equations\nbsp{}eqref:eq-problem-2d to
-\begin{align*}
-    & \phi_{xx}+\phi_{zz}=0,\\
-    & \zeta_t = -\phi_z & \text{на }z=\zeta(x,t),
-\end{align*}
-solution to which is written as
-\begin{equation*}
-    \phi(x,z,t)
-    =
-    -\InverseFourierY{
-        \frac{e^{2\pi u z}}{2\pi u}
-        \FourierY{\zeta_t}{u}
-    }{x}
-    .
-\end{equation*}
-Propagating wave profile is defined as \(\zeta(x,t)=A\cos(2\pi(kx-t))\).
-Plugging this formula into\nbsp{}eqref:eq-solution-2d yields
-\(\phi(x,z,t)=-\frac{A}{k}\sin(2\pi(kx-t))\Sinh{2\pi{k}{z}}\). In order to
-reduce it to the formula from linear wave theory, rewrite hyperbolic sine in
-exponential form and discard the term containing \(e^{-2\pi{k}{z}}\) as
-contradicting condition
-\(\phi\underset{z\rightarrow-\infty}{\longrightarrow}0\). Taking real part of
-the resulting formula yields
-\(\phi(x,z,t)=\frac{A}{k}e^{2\pi{k}{z}}\sin(2\pi(kx-t))\), which corresponds to
-the known formula from linear wave theory. Similarly, under small-amplitude
-waves assumption the formula for finite depth
-fluid\nbsp{}eqref:eq-solution-2d-full is reduced to
-\begin{equation*}
-    \phi(x,z,t)
-    =
-    -\InverseFourierY{
-        \frac{\Sinh{2\pi u (z+h)}}{2\pi u \Sinh{2\pi u h}}
-        \FourierY{\zeta_t}{u}
-    }{x}.
-\end{equation*}
-Substituting \(\zeta(x,t)\) with propagating plain wave profile formula yields
-\begin{equation}
-    \label{eq-solution-2d-linear}
-    \phi(x,z,t)=\frac{A}{k}
-    \frac{\Sinh{2 \pi k (z+h)}}{ \Sinh{2 \pi k h} }
-    \sin(2 \pi (k x-t)),
-\end{equation}
-which corresponds to the formula from linear wave theory for finite depth fluid.
+#+NAME: fig-skew-normal-1
+#+begin_src R :file build/skew-normal-1.pdf
+source(file.path("R", "common.R"))
+x <- seq(-3, 3, length.out=100)
+params <- data.frame(
+  skewness = c(0.00, 0.52, 0.00, 0.52),
+  kurtosis = c(0.00, 0.00, 0.70, 0.70),
+  linetypes = c("solid", "dashed", "dotdash", "dotted")
+)
+arma.skew_normal_1_plot(x, params)
+legend(
+  "topleft",
+  mapply(
+    function (s, k) {
+      as.expression(bquote(list(
+        gamma[1] == .(arma.fmt(s, 2)),
+        gamma[2] == .(arma.fmt(k, 2))
+      )))
+    },
+    params$skewness,
+    params$kurtosis
+  ),
+  lty = paste(params$linetypes)
+)
+#+end_src
 
-Different forms of Laplace equation solutions, in which decaying exponent is
-written with either "+" or "-" signs, may cause incompatibilities between
-formulae from linear wave theory and formulae derived in this work, where
-\(\sinh\) is used instead of \(\cosh\). Equality
-\(\frac{\Sinh{2\pi{k}(z+h)}}{\Sinh{2\pi{k}{h}}}\approx\frac{\sinh(2\pi{k}(z+h))}{\sinh(2\pi{k}{h})}\)
-becomes strict on the free surface, and difference between left-hand and
-right-hand sides increases when approaching sea bottom (for sufficiently large
-depth difference near free surface is negligible). So, for sufficiently large
-depth any function (\(\cosh\) or \(\sinh\)) may be used for velocity potential
-computation near free surface.
+#+caption[Graph of PDF of GCS-based distribution]: Graph of probability density function of GCS-based distribution law for different values of skewness \(\gamma_1\) and kurtosis \(\gamma_2\).
+#+name: fig-skew-normal-1
+#+RESULTS: fig-skew-normal-1
+[[file:build/skew-normal-1.pdf]]
 
-Reducing\nbsp{}eqref:eq-solution-2d and\nbsp{}eqref:eq-solution-2d-full to the
-known formulae from linear wave theory shows, that formula for infinite
-depth\nbsp{}eqref:eq-solution-2d is not suitable to compute velocity potentials
-with Fourier method, because it does not have symmetry, which is required for
-Fourier transform. However, formula for finite depth can be used instead by
-setting \(h\) to some characteristic water depth. For standing wave reducing to
-linear wave theory formulae is made under the same assumptions.
-*** Three-dimensional velocity field
-Three-dimensional version of\nbsp{}eqref:eq-problem is written as
+**** Skew-normal distribution.
+Alternative approach is to approximate distribution of sea wavy surface
+elevation by skew-normal (SN) distribution:
 \begin{align}
-    \label{eq-problem-3d}
-    & \phi_{xx} + \phi_{yy} + \phi_{zz} = 0,\\
-    & \zeta_t + \zeta_x\phi_x + \zeta_y\phi_y
-    =
-    \frac{\zeta_x}{\SqrtZeta{1 + \zeta_x^2 + \zeta_y^2}} \phi_x
-    +\frac{\zeta_y}{\SqrtZeta{1 + \zeta_x^2 + \zeta_y^2}} \phi_y
-    - \phi_z, & \text{на }z=\zeta(x,y,t).\nonumber
+    \label{eq-skew-normal-2}
+    F(z; \alpha) & = \frac{1}{2}
+   \mathrm{erfc}\left[-\frac{z}{\sqrt{2}}\right]-2 T(z,\alpha ), \nonumber \\
+    f(z; \alpha) & = \frac{e^{-\frac{z^2}{2}}}{\sqrt{2 \pi }}
+   \mathrm{erfc}\left[-\frac{\alpha z}{\sqrt{2}}\right],
 \end{align}
-Again, use Fourier method to solve it. Applying Fourier transform to both sides
-of Laplace equation yields
-\begin{equation*}
-    -4 \pi^2 \left( u^2 + v^2 + w^2 \right)
-    \FourierY{\phi(x,y,z)}{u,v,w} = 0,
-\end{equation*}
-hence \(w=\pm{i}\sqrt{u^2+v^2}\). We seek solution in the form of inverse Fourier
-transform \(\phi(x,y,z)=\InverseFourierY{E(u,v,w)}{x,y,z}\). Plugging
-\(w=i\sqrt{u^2+v^2}=i\Kveclen\) into the formula yields
-\begin{equation*}
-    \phi(x,y,z) = \InverseFourierY{
-        \left(
-            C_1 e^{2\pi \Kveclen z}
-            -C_2 e^{-2\pi \Kveclen z}
-        \right)
-        E(u,v)
-    }{x,y}.
-\end{equation*}
-Plugging \(\phi\) into the boundary condition on the sea bottom (analogous to
-two-dimensional case) yields
-\begin{equation}
-    \label{eq-guessed-sol-3d}
-    \phi(x,y,z) = \InverseFourierY{
-        \Sinh{2\pi \Kveclen (z+h)} E(u,v)
-    }{x,y}.
-\end{equation}
-Plugging \(\phi\) into the boundary condition on the free surface yields
-\begin{equation*}
-    \arraycolsep=1.4pt
-    \begin{array}{rl}
-        \zeta_t = & i f_1(x,y) \InverseFourierY{2 \pi u \Sinh{2\pi \Kveclen (z+h)}E(u,v)}{x,y} \\
-        + & i f_2(x,y) \InverseFourierY{2 \pi v \Sinh{2\pi \Kveclen (z+h)}E(u,v)}{x,y} \\
-        - & \InverseFourierY{2 \pi \Kveclen \SinhX{2\pi \Kveclen (z+h)}E(u,v)}{x,y}
-    \end{array}
-\end{equation*}
-where \(f_1(x,y)={\zeta_x}/{\SqrtZeta{1+\zeta_x^2+\zeta_y^2}}-\zeta_x\) and
-\(f_2(x,y)={\zeta_y}/{\SqrtZeta{1+\zeta_x^2+\zeta_y^2}}-\zeta_y\).
-
-Like in section\nbsp{}[[#sec-pressure-2d]] we assume that
-\(\Sinh{2\pi{u}(z+h)}\approx\SinhX{2\pi{u}(z+h)}\) near free surface, but in
-three-dimensional case this is not enough to solve the problem. In order to get
-explicit formula for coefficients \(E\) we need to assume, that all Fourier
-transforms in the equation have radially symmetric kernels, i.e.\nbsp{}replace
-\(u\) and \(v\) with \(\Kveclen\). There are two points supporting this
-assumption. First, in numerical implementation integration is done over positive
-wave numbers, so the sign of \(u\) and \(v\) does not affect the solution.
-Second, the growth rate of \(\cosh\) term of the integral kernel is much higher
-than the one of \(u\) or \(\Kveclen\), so the substitution has small effect on
-the magnitude of the solution. Despite these two points, a use of more
-mathematically rigorous approach would be preferable.
-
-Making the replacement, applying Fourier transform to both sides of the equation
-and plugging the result into\nbsp{}eqref:eq-guessed-sol-3d yields formula for
-\(\phi\):
-\begin{equation}
-    \label{eq-phi-3d}
-    \phi(x,y,z,t) = \InverseFourierY{
-        \frac{ \Sinh{\smash{2\pi \Kveclen (z+h)}} }{ 2\pi\Kveclen }
-        \frac{ \FourierY{ \zeta_t / \left( i f_1(x,y) + i f_2(x,y) - 1 \right)}{u,v} }
-        { \FourierY{\mathcal{D}_3\left( x,y,\zeta\left(x,y\right) \right)}{u,v} }
-    }{x,y},
-\end{equation}
-where
-\(\FourierY{\mathcal{D}_3\left(x,y,z\right)}{u,v}=\Sinh{\smash{2\pi\Kveclen{}z}}\).
+where \(T\)\nbsp{}--- Owen \(T\)-function\nbsp{}cite:owen1956tables. Using this
+formula it is impossible to specify skewness and kurtosis separately\nbsp{}---
+both values are adjusted via \(\alpha\) parameter. The only advantage of the
+formula is its relative computational simplicity: this function is available in
+some programmes and mathematical libraries. Its graph for different values of
+\(\alpha\) is shown in fig.\nbsp{}[[fig-skew-normal-2]].
 
-*** Conclusion
-Formulae derived in this section allow to numerically compute velocity potential
-field (and hence pressure field) near discretely or mathematically given wavy
-sea surface, bypassing assumptions of linear wavy theory and theory of small
-amplitude waves.
-
-** Modelling non-linearity of sea waves
-ARMA model allows to model asymmetry of wave elevation distribution,
-i.e.\nbsp{}generate sea waves, distribution of \(z\)-coordinate of which has
-non-nought kurtosis and asymmetry. Such distribution is inherent to real sea
-waves\nbsp{}cite:longuet1963nonlinear. Wave asymmetry is modelled by non-linear
-inertia-less transform (NIT) of stochastic process, however, transforming
-resulting wavy surface means transforming initial ACF. In order to alleviate
-this, ACF must be preliminary transformed as shown
-in\nbsp{}cite:boukhanovsky1997thesis.
-
-**** Wavy surface transformation.
-Explicit formula \(z=f(y)\) that transforms wavy surface to desired
-one-dimensional distribution \(F(z)\) is the solution of non-linear
-transcendental equation \(F(z)=\Phi(y)\), where \(\Phi(y)\)\nbsp{}---
-one-dimensional Gaussian distribution. Since distribution of wave elevation is
-often given by some approximation based on /in situ/ data, this equation is
-solved numerically with respect to \(z_k\) in each grid point \(y_k|_{k=0}^N\)
-of generated wavy surface. In this case equation is rewritten as
-\begin{equation}
-    \label{eq-distribution-transformation}
-    F(z_k)
-    =
-    \frac{1}{\sqrt{2\pi}}
-    \int\limits_0^{y_k} \exp\left[ -\frac{t^2}{2} \right] dt
-    .
-\end{equation}
-Since, distribution functions are monotonic, the simplest interval halving
-(bisection) numerical method is used to solve this equation.
-
-**** Preliminary ACF transformation.
-In order to transform ACF \(\gamma_z\) of the process, it is expanded in series
-of Hermite polynomials (Gram---Charlier series)
-\begin{equation*}
-    \gamma_z \left( \vec u \right)
-    =
-    \sum\limits_{m=0}^{\infty}
-    C_m^2 \frac{\gamma_y^m \left( \vec u \right)}{m!},
-\end{equation*}
-where
-\begin{equation*}
-    C_m = \frac{1}{\sqrt{2\pi}}
-  \int\limits_{0}^\infty
-    f(y) H_m(y) \exp\left[ -\frac{y^2}{2} \right],
-\end{equation*}
-\(H_m\)\nbsp{}--- Hermite polynomial, and \(f(y)\)\nbsp{}--- solution to
-equation\nbsp{}eqref:eq-distribution-transformation. Plugging polynomial
-approximation \(f(y)\approx\sum\limits_{i}d_{i}y^i\) and analytic formulae for
-Hermite polynomial yields
-\begin{equation*}
-    \frac{1}{\sqrt{2\pi}}
-    \int\limits_\infty^\infty
-    y^k \exp\left[ -\frac{y^2}{2} \right]
-    =
-    \begin{cases}
-        (k-1)!! & \text{if }k\text{ is even},\\
-        0       & \text{if }k\text{ is odd},
-    \end{cases}
-\end{equation*}
-which simplifies the former equation. Optimal number of coefficients \(C_m\) is
-determined by computing them sequentially and stopping when variances of both
-fields become equal with desired accuracy \(\epsilon\):
-\begin{equation}
-    \label{eq-nit-error}
-    \left| \Var{z} - \sum\limits_{k=0}^m
-    \frac{C_k^2}{k!} \right| \leq \epsilon.
-\end{equation}
+#+name: fig-skew-normal-2
+#+begin_src R :file build/skew-normal-2.pdf
+source(file.path("R", "common.R"))
+x <- seq(-3, 3, length.out=100)
+alpha <- c(0.00, 0.87, 2.25, 4.90)
+params <- data.frame(
+  alpha = alpha,
+  skewness = arma.bits.skewness_2(alpha),
+  kurtosis = arma.bits.kurtosis_2(alpha),
+  linetypes = c("solid", "dashed", "dotdash", "dotted")
+)
+arma.skew_normal_2_plot(x, params)
+legend(
+  "topleft",
+  mapply(
+    function (a, s, k) {
+      as.expression(bquote(list(
+        alpha == .(arma.fmt(a, 2)),
+        gamma[1] == .(arma.fmt(s, 2)),
+        gamma[2] == .(arma.fmt(k, 2))
+      )))
+    },
+    params$alpha,
+    params$skewness,
+    params$kurtosis
+  ),
+  lty = paste(params$linetypes)
+)
+#+end_src
 
-In\nbsp{}cite:boukhanovsky1997thesis the author suggests using polynomial
-approximation \(f(y)\) also for wavy surface transformation, however, in
-practice sea surface realisation often contains points, in which
-\(z\)-coordinate is beyond the limits of the approximation, leading to sharp
-precision decrease. In these points it is more efficient to solve
-equation\nbsp{}eqref:eq-distribution-transformation by bisection method. Using
-the same approximation in Gram---Charlier series does not lead to such errors.
+#+caption[Graph of PDF of SN distribution]: Graph of PDF of skew-normal distribution for different values of skewness coefficient \(\alpha\).
+#+name: fig-skew-normal-2
+#+RESULTS: fig-skew-normal-2
+[[file:build/skew-normal-2.pdf]]
 
-* Numerical methods and experimental results
 ** The shape of ACF for different types of waves
 :PROPERTIES:
 :CUSTOM_ID: sec-wave-acfs
@@ -1190,361 +899,675 @@ steps.
 - Move maximum value of the resulting function to the origin by using
   trigonometric identities to shift the phase.
 
-** Additional formulae, methods and algorithms for ARMA model
+** ARMA model verification
 :PROPERTIES:
-:CUSTOM_ID: sec-arma-algorithms
+:CUSTOM_ID: sec-verification
 :END:
 
-*** Wave elevation distribution approximation
-One of the parameters of sea wavy surface generator is probability density
-function (PDF) of the surface elevation. This distribution is given by either
-polynomial approximation of /in situ/ data or analytic formula.
+In\nbsp{}cite:degtyarev2011modelling,degtyarev2013synoptic,boukhanovsky1997thesis
+AR model the following items are verified experimentally:
+- probability distributions of different wave characteristics (wave heights,
+  lengths, crests, periods, slopes, three-dimensionality),
+- dispersion relation,
+- retention of integral characteristics for mixed wave sea state.
+In this work both AR and MA model are verified by comparing probability
+distributions of different wave characteristics.
 
-**** Gram---Charlier series expansion.
-In\nbsp{}cite:huang1980experimental the authors experimentally show, that PDF of
-sea surface elevation is distinguished from normal distribution by non-nought
-kurtosis and skewness. In\nbsp{}cite:rozhkov1996theory the authors show, that
-this type of PDF expands in Gram---Charlier series (GCS):
-\begin{align}
-    \label{eq-skew-normal-1}
-    & F(z; \mu=0, \sigma=1, \gamma_1, \gamma_2) \approx
-    \Phi(z; \mu, \sigma) \frac{2 + \gamma_2}{2}
-    - \frac{2}{3} \phi(z; \mu, \sigma)
-    \left(\gamma_2 z^3+\gamma_1
-    \left(2 z^2+1\right)\right) \nonumber
-    \\
-    & f(z; \mu, \sigma, \gamma_1, \gamma_2) \approx
-    \phi(z; \mu, \sigma)
-    \left[
-        1+
-        \frac{1}{3!} \gamma_1 H_3 \left(\frac{z-\mu}{\sigma}\right)
-        + \frac{1}{4!} \gamma_2 H_4 \left(\frac{z-\mu}{\sigma}\right)
-    \right],
-\end{align}
-where \(\Phi(z)\)\nbsp{}--- cumulative density function (CDF) of normal
-distribution, \(\phi\)\nbsp{}--- PDF of normal distribution,
-\(\gamma_1\)\nbsp{}--- skewness, \(\gamma_2\)\nbsp{}--- kurtosis,
-\(f\)\nbsp{}--- PDF, \(F\)\nbsp{}--- cumulative distribution function (CDF).
-According to\nbsp{}cite:rozhkov1990stochastic for sea waves skewness is selected
-from interval \(0.1\leq\gamma_1\leq{0.52}]\) and kurtosis from interval
-\(0.1\leq\gamma_2\leq{0.7}\). Family of probability density functions for
-different parameters is shown in fig.\nbsp{}[[fig-skew-normal-1]].
+*** Verification of wavy surface integral characteristics
+In\nbsp{}cite:rozhkov1990stochastic the authors show that several sea wave
+characteristics (listed in table\nbsp{}[[tab-weibull-shape]]) have Weibull
+distribution, and wavy surface elevation has Gaussian distribution. In order to
+verify that distributions corresponding to generated realisation are correct,
+quantile-quantile plots are used (plots where analytic quantile values are used
+for \(OX\) axis and estimated quantile values for \(OY\) axis). If the estimated
+distribution matches analytic then the graph has the form of the straight line.
+Tails of the graph may diverge from the straight line, because they can not be
+reliably estimated from the finite-size realisation.
 
-#+NAME: fig-skew-normal-1
-#+begin_src R :file build/skew-normal-1.pdf
+#+name: tab-weibull-shape
+#+caption[Values of Weibull shape parameter]:
+#+caption: Values of Weibull shape parameter for different wave characteristics.
+#+attr_latex: :booktabs t
+| Characteristic       | Weibull shape (\(k\)) |
+|----------------------+-----------------------|
+|                      | <l>                   |
+| Wave height          | 2                     |
+| Wave length          | 2.3                   |
+| Crest length         | 2.3                   |
+| Wave period          | 3                     |
+| Wave slope           | 2.5                   |
+| Three-dimensionality | 2.5                   |
+
+Verification was performed for standing and propagating waves. The corresponding
+ACFs and quantile-quantile plots of wave characteristics distributions are shown
+in
+fig.\nbsp{}[[propagating-wave-distributions]],\nbsp{}[[standing-wave-distributions]],\nbsp{}[[acf-slices]].
+
+#+name: propagating-wave-distributions
+#+begin_src R :file build/propagating-wave-qqplots.pdf
 source(file.path("R", "common.R"))
-x <- seq(-3, 3, length.out=100)
-params <- data.frame(
-  skewness = c(0.00, 0.52, 0.00, 0.52),
-  kurtosis = c(0.00, 0.00, 0.70, 0.70),
-  linetypes = c("solid", "dashed", "dotdash", "dotted")
-)
-arma.skew_normal_1_plot(x, params)
-legend(
-  "topleft",
-  mapply(
-    function (s, k) {
-      as.expression(bquote(list(
-        gamma[1] == .(arma.fmt(s, 2)),
-        gamma[2] == .(arma.fmt(k, 2))
-      )))
-    },
-    params$skewness,
-    params$kurtosis
-  ),
-  lty = paste(params$linetypes)
+par(pty="s", mfrow=c(2, 2))
+arma.qqplot_grid(
+  file.path("build", "propagating_wave"),
+  c("elevation", "heights_y", "lengths_y", "periods"),
+  c("elevation", "height Y", "length Y", "period"),
+  xlab="x",
+  ylab="y"
 )
 #+end_src
 
-#+caption[Graph of PDF of GCS-based distribution]: Graph of probability density function of GCS-based distribution law for different values of skewness \(\gamma_1\) and kurtosis \(\gamma_2\).
-#+name: fig-skew-normal-1
-#+RESULTS: fig-skew-normal-1
-[[file:build/skew-normal-1.pdf]]
-
-**** Skew-normal distribution.
-Alternative approach is to approximate distribution of sea wavy surface
-elevation by skew-normal (SN) distribution:
-\begin{align}
-    \label{eq-skew-normal-2}
-    F(z; \alpha) & = \frac{1}{2}
-   \mathrm{erfc}\left[-\frac{z}{\sqrt{2}}\right]-2 T(z,\alpha ), \nonumber \\
-    f(z; \alpha) & = \frac{e^{-\frac{z^2}{2}}}{\sqrt{2 \pi }}
-   \mathrm{erfc}\left[-\frac{\alpha z}{\sqrt{2}}\right],
-\end{align}
-where \(T\)\nbsp{}--- Owen \(T\)-function\nbsp{}cite:owen1956tables. Using this
-formula it is impossible to specify skewness and kurtosis separately\nbsp{}---
-both values are adjusted via \(\alpha\) parameter. The only advantage of the
-formula is its relative computational simplicity: this function is available in
-some programmes and mathematical libraries. Its graph for different values of
-\(\alpha\) is shown in fig.\nbsp{}[[fig-skew-normal-2]].
+#+caption[Quantile-quantile plots for propagating waves]:
+#+caption: Quantile-quantile plots for propagating waves.
+#+name: propagating-wave-distributions
+#+RESULTS: propagating-wave-distributions
+[[file:build/propagating-wave-qqplots.pdf]]
 
-#+name: fig-skew-normal-2
-#+begin_src R :file build/skew-normal-2.pdf
+#+name: standing-wave-distributions
+#+begin_src R :file build/standing-wave-qqplots.pdf
 source(file.path("R", "common.R"))
-x <- seq(-3, 3, length.out=100)
-alpha <- c(0.00, 0.87, 2.25, 4.90)
-params <- data.frame(
-  alpha = alpha,
-  skewness = arma.bits.skewness_2(alpha),
-  kurtosis = arma.bits.kurtosis_2(alpha),
-  linetypes = c("solid", "dashed", "dotdash", "dotted")
-)
-arma.skew_normal_2_plot(x, params)
-legend(
-  "topleft",
-  mapply(
-    function (a, s, k) {
-      as.expression(bquote(list(
-        alpha == .(arma.fmt(a, 2)),
-        gamma[1] == .(arma.fmt(s, 2)),
-        gamma[2] == .(arma.fmt(k, 2))
-      )))
-    },
-    params$alpha,
-    params$skewness,
-    params$kurtosis
-  ),
-  lty = paste(params$linetypes)
+par(pty="s", mfrow=c(2, 2))
+arma.qqplot_grid(
+  file.path("build", "standing_wave"),
+  c("elevation", "heights_y", "lengths_y", "periods"),
+  c("elevation", "height Y", "length Y", "period"),
+  xlab="x",
+  ylab="y"
 )
 #+end_src
 
-#+caption[Graph of PDF of SN distribution]: Graph of PDF of skew-normal distribution for different values of skewness coefficient \(\alpha\).
-#+name: fig-skew-normal-2
-#+RESULTS: fig-skew-normal-2
-[[file:build/skew-normal-2.pdf]]
-
-**** Evaluation.                                                :noexport:
-Equation\nbsp{}eqref:eq-distribution-transformation with selected wave elevation
-distribution may be solved either in every point of generated wavy surface,
-which gives the most accurate results, or in every fixed grid point
-interpolating result via least-squares (LS) polynomial. In the second case
-precision is lower. For example, interpolating 12^th order polynomial on a fixed
-grid of 500 points on interval \(-5\sigma_z\leq{z}\leq{5}\sigma_z\) gives error
-of \(\approx{0.43}\cdot10^{-3}\). Increasing polynomial order leads to either
-numeric overflows during LS interpolation, or more coefficient close to nought;
-increasing the size of the grid has insignificant effect on the result. In the
-majority of cases three Gram---Charlier series coefficients is enough to
-transform ACF; relative error without interpolation is \(10^{-5}\).
-
-*** White noise generation algorithm
-In order to eliminate periodicity from wavy surface generated by sea wave model,
-it is imperative to use PRNG with sufficiently large period to generate white
-noise. Parallel Mersenne Twister\nbsp{}cite:matsumoto1998mersenne with a period
-of \(2^{19937}-1\) is used as a generator in this work. It allows to produce
-aperiodic sea wavy surface realisations in any practical usage scenarios.
-
-There is no guarantee that multiple PRNGs executed in parallel threads with
-distinct initial states produce uncorrelated pseudo-random number sequences, and
-algorithm of dynamic creation of Mersenne
-Twisters\nbsp{}cite:matsumoto1998dynamic is used to provide such guarantee. The
-essence of the algorithm is to find matrices of initial generator states, that
-give maximally uncorrelated pseudo-random number sequences when Mersenne
-Twisters are executed in parallel with these initial states. Since finding such
-initial states consumes considerable amount of processor time, vector of initial
-states is created preliminary with knowingly larger number of parallel threads
-and saved to a file, which is then read before starting white noise generation.
-
-*** Wavy surface generation algorithm
-In ARMA model value of wavy surface elevation at a particular point depends on
-previous in space and time points, as a result the so called /ramp-up interval/
-(see fig.\nbsp{}[[fig-ramp-up-interval]]), in which realisation does not correspond
-to specified ACF, forms in the beginning of the realisation. There are several
-solutions to this problem which depend on the simulation context.
-
-If realisation is used in the context of ship stability simulation without
-manoeuvring, ramp-up interval will not affect results of the simulation, because
-it is located on the border (too far away from the studied marine object). If
-ship stability with manoeuvring is studied, then the interval may be simply
-discarded from the realisation (the size of the interval approximately equals
-the number of AR coefficients in each dimension). However, this may lead to loss
-of a very large number of points, because discarding occurs for each of three
-dimensions. Alternative approach is to generate sea wavy surface on ramp-up
-interval with LH model and generate the rest of the realisation with ARMA model.
-
-Algorithm of wavy surface generation is data-parallel: realisation is divided
-into equal parts each of which is generated independently. For MA model there is
-ramp-up interval at the beginning of each realisation. To eliminate it
-/overlap-add/
-method\nbsp{}cite:oppenheim1989discrete,svoboda2011efficient,pavel2013algorithms
-(a popular method in signal processing) is used. The essence of the method is to
-add another interval, size of which is equal to the ramp-up interval size, to
-the end of each part. Then wavy surface is generated in each point of each part
-(including points from the added interval), the interval at the end of part
-\(N\) is superimposed on the ramp-up interval at the beginning of the part
-\(N+1\), and values in corresponding points are added. For AR model, parts are
-generated taking into account autoregressive dependencies between them. Detailed
-explanation of parallel and distributed algorithms is given in
-sections\nbsp{}[[#sec-parallel]] and\nbsp{}[[#sec-distributed]].
+#+caption[Quantile-quantile plots for standing waves]:
+#+caption: Quantile-quantile plots for standing waves.
+#+name: standing-wave-distributions
+#+RESULTS: standing-wave-distributions
+[[file:build/standing-wave-qqplots.pdf]]
 
-#+name: fig-ramp-up-interval
-#+begin_src R :file build/ramp-up-interval.pdf
+#+name: acf-slices
+#+header: :width 6 :height 9
+#+begin_src R :file build/acf-slices.pdf
 source(file.path("R", "common.R"))
-arma.plot_ramp_up_interval()
+propagating_acf <- read.csv(file.path("build", "propagating_wave", "acf.csv"))
+standing_acf <- read.csv(file.path("build", "standing_wave", "acf.csv"))
+par(mfrow=c(5, 2), mar=c(0,0,0,0))
+for (i in seq(0, 4)) {
+  arma.wavy_plot(standing_acf, i, zlim=c(-5,5))
+  arma.wavy_plot(propagating_acf, i, zlim=c(-5,5))
+}
 #+end_src
 
-#+caption[Ramp-up interval at the beginning of the realisation]:
-#+caption: Ramp-up interval at the beginning of the realisation.
-#+name: fig-ramp-up-interval
-#+RESULTS: fig-ramp-up-interval
-[[file:build/ramp-up-interval.pdf]]
+#+caption[ACF time slices for standing and propagating waves]:
+#+caption: ACF time slices for standing (left column) and propagating waves
+#+caption: (right column).
+#+name: acf-slices
+#+RESULTS: acf-slices
+[[file:build/acf-slices.pdf]]
 
-*** Velocity potential normalisation formulae
-:PROPERTIES:
-:CUSTOM_ID: sec-compute-delta
-:END:
+Graph tails in fig.\nbsp{}[[propagating-wave-distributions]] deviate from original
+distribution for individual wave characteristics, because every wave have to be
+extracted from the resulting wavy surface to measure its length, period and
+height. There is no algorithm that guarantees correct extraction of all waves,
+because they may and often overlap each other. Weibull distribution right tail
+represents infrequently occurring waves, so it deviates more than left tail.
 
-In solutions\nbsp{}eqref:eq-solution-2d and\nbsp{}eqref:eq-solution-2d-full to
-two-dimensional pressure determination problem there are functions
-\(\Fun{z}=\InverseFourierY{e^{2\pi{u}{z}}}{x}\) and
-\(\FunSecond{z}=\InverseFourierY{\Sinh{2\pi{u}{z}}}{x}\) which has multiple
-analytic representations and are difficult to compute. Each function is a
-Fourier transform of linear combination of exponents which reduces to poorly
-defined Dirac delta function of a complex argument (see
-table\nbsp{}[[tab-delta-functions]]). The usual way of handling this type of
-functions is to write them as multiplication of Dirac delta functions of real
-and imaginary part, however, this approach does not work here, because applying
-inverse Fourier transform to this representation does not produce exponent,
-which severely warp resulting velocity field. In order to get unique analytic
-definition normalisation factor \(1/\Sinh{2\pi{u}{h}}\) (which is also included
-in formula for \(E(u)\)) may be used. Despite the fact that normalisation allows
-to obtain adequate velocity potential field, numerical experiments show that
-there is little difference between this field and the one produced by formulae
-from linear wave theory, in which terms with \(\zeta\) are omitted. As a result,
-the formula for three-dimensional case was not derived.
+Correspondence rate for standing waves (fig.\nbsp{}[[standing-wave-distributions]])
+is lower for height and length, roughly the same for surface elevation and
+higher for wave period distribution tails. Lower correspondence degree for
+length and height may be attributed to the fact that Weibull distributions were
+obtained empirically for sea waves which are typically propagating, and
+distributions may be different for standings waves. Higher correspondence degree
+for wave periods is attributed to the fact that wave periods of standing waves
+are extracted more precisely as the waves do not move outside simulated wavy
+surface region. The same correspondence degree for wave elevation is obtained,
+because this is the characteristic of the wavy surface (and corresponding AR or
+MA process) and is not affected by the type of waves.
 
-#+name: tab-delta-functions
-#+caption[Formulae for computing \(\Fun{z}\) and \(\FunSecond{z}\)]:
-#+caption: Formulae for computing \(\Fun{z}\) and \(\FunSecond{z}\) from
-#+caption: sec.\nbsp{}[[#sec-pressure-2d]], that use normalisation to eliminate
-#+caption: uncertainty from definition of Dirac delta function of complex argument.
-#+attr_latex: :booktabs t
-| Function          | Without normalisation                                        | Normalised                                                                                                                             |
-|-------------------+--------------------------------------------------------------+----------------------------------------------------------------------------------------------------------------------------------------|
-| \(\Fun{z}\)       | \(\delta (x+i z)\)                                           | \(\frac{1}{2 h}\mathrm{sech}\left(\frac{\pi  (x-i (h+z))}{2 h}\right)\)                                                                |
-| \(\FunSecond{z}\) | \(\frac{1}{2}\left[\delta (x-i z) + \delta (x+i z) \right]\) | \(\frac{1}{4 h}\left[\text{sech}\left(\frac{\pi  (x-i (h+z))}{2 h}\right)+\text{sech}\left(\frac{\pi  (x+i(h+z))}{2 h}\right)\right]\) |
+*** Verification of nonlinear inertialess transformation
+In order to measure the effect of NIT on the shape of the resulting wavy
+surface, three realisations were generated:
+- realisation with Gaussian distribution (without NIT),
+- realisation with Gram---Charlier series (GCS) based distribution, and
+- realisation with skew normal distribution.
+The seed of PRNG was set to be the same for all programme executions to make ARMA
+model produce the same values for each realisation. There we two experiments:
+for standing and propagating waves with ACFs given by formulae from
+section\nbsp{}[[#sec-wave-acfs]].
 
-** ARMA model verification
-:PROPERTIES:
-:CUSTOM_ID: sec-verification
-:END:
+While the experiment showed that applying NIT with GCS-based distribution
+increases wave steepness, the same is not true for skew normal distribution
+(fig.\nbsp{}[[fig-nit]]). Using this distribution results in wavy surface each
+\(z\)-coordinate of which is always greater or equal to nought. So, skew normal
+distribution is unsuitable for NIT. NIT increases the wave height and steepness
+of both standing and propagating waves. Increasing either skewness or kurtosis
+of GCS-based distribution increases both wave height and steepness. The error of
+ACF approximation (eq.\nbsp{}eqref:eq-nit-error) ranges from 0.20 for GCS-based
+distribution to 0.70 for skew normal distribution (table\nbsp{}[[tab-nit-error]]).
 
-In\nbsp{}cite:degtyarev2011modelling,degtyarev2013synoptic,boukhanovsky1997thesis
-AR model the following items are verified experimentally:
-- probability distributions of different wave characteristics (wave heights,
-  lengths, crests, periods, slopes, three-dimensionality),
-- dispersion relation,
-- retention of integral characteristics for mixed wave sea state.
-In this work both AR and MA model are verified by comparing probability
-distributions of different wave characteristics.
+#+name: fig-nit
+#+header: :width 7 :height 7
+#+begin_src R :file build/nit.pdf
+source(file.path("R", "nonlinear.R"))
+par(mfrow=c(2, 1), mar=c(4,4,4,0.5), family='serif')
+args <- list(
+  graphs=c('Gaussian', 'Gram—Charlier', 'Skew normal'),
+  linetypes=c('solid', 'dashed', 'dotted')
+)
+args$title <- 'Propagating waves'
+arma.plot_nonlinear(file.path("build", "nit-propagating"), args)
+args$title <- 'Standing waves'
+arma.plot_nonlinear(file.path("build", "nit-standing"), args)
+#+end_src
 
-*** Verification of wavy surface integral characteristics
-In\nbsp{}cite:rozhkov1990stochastic the authors show that several sea wave
-characteristics (listed in table\nbsp{}[[tab-weibull-shape]]) have Weibull
-distribution, and wavy surface elevation has Gaussian distribution. In order to
-verify that distributions corresponding to generated realisation are correct,
-quantile-quantile plots are used (plots where analytic quantile values are used
-for \(OX\) axis and estimated quantile values for \(OY\) axis). If the estimated
-distribution matches analytic then the graph has the form of the straight line.
-Tails of the graph may diverge from the straight line, because they can not be
-reliably estimated from the finite-size realisation.
+#+name: fig-nit
+#+caption[Wavy surface slices with different distributions]:
+#+caption: Wavy surface slices with different distributions
+#+caption: of wave elevation (Gaussian, GCS-based and SN).
+#+RESULTS: fig-nit
+[[file:build/nit.pdf]]
 
-#+name: tab-weibull-shape
-#+caption[Values of Weibull shape parameter]:
-#+caption: Values of Weibull shape parameter for different wave characteristics.
+#+name: tab-nit-error
+#+caption[Errors of ACF approximations]:
+#+caption: Errors of ACF approximations (the difference of variances) for
+#+caption: different wave elevation distributions. \(N\) is the number of
+#+caption: coefficients of ACF approximation.
 #+attr_latex: :booktabs t
-| Characteristic       | Weibull shape (\(k\)) |
-|----------------------+-----------------------|
-|                      | <l>                   |
-| Wave height          | 2                     |
-| Wave length          | 2.3                   |
-| Crest length         | 2.3                   |
-| Wave period          | 3                     |
-| Wave slope           | 2.5                   |
-| Three-dimensionality | 2.5                   |
+| Wave type   | Distribution | \(\gamma_1\) | \(\gamma_2\) | \(\alpha\) | Error | \(N\) | Wave height |
+|-------------+--------------+--------------+--------------+------------+-------+-------+-------------|
+| propagating | Gaussian     |              |              |            |       |       |        2.41 |
+| propagating | GCS-based    |         2.25 |          0.4 |            |  0.20 |     2 |        2.75 |
+| propagating | skew normal  |              |              |          1 |  0.70 |     3 |        1.37 |
+| standing    | Gaussian     |              |              |            |       |       |        1.73 |
+| standing    | GCS-based    |         2.25 |          0.4 |            |  0.26 |     2 |        1.96 |
+| standing    | skew normal  |              |              |          1 |  0.70 |     3 |        0.94 |
 
-Verification was performed for standing and propagating waves. The corresponding
-ACFs and quantile-quantile plots of wave characteristics distributions are shown
-in
-fig.\nbsp{}[[propagating-wave-distributions]],\nbsp{}[[standing-wave-distributions]],\nbsp{}[[acf-slices]].
+To summarise, the only test case that showed acceptable results is realisation
+with GCS-based distribution for both standing and propagating waves. Skew normal
+distribution warps wavy surface for both types of waves. GCS-based realisations
+have large error of ACF approximation, which results in increase of wave height.
+The reason for the large error is that GCS approximations are not accurate as
+they do not converge for all possible
+functions\nbsp{}cite:wallace1958asymptotic. Despite the large error, the change
+in wave height is small (table\nbsp{}[[tab-nit-error]]).
 
-#+name: propagating-wave-distributions
-#+begin_src R :file build/propagating-wave-qqplots.pdf
-source(file.path("R", "common.R"))
-par(pty="s", mfrow=c(2, 2))
-arma.qqplot_grid(
-  file.path("build", "propagating_wave"),
-  c("elevation", "heights_y", "lengths_y", "periods"),
-  c("elevation", "height Y", "length Y", "period"),
-  xlab="x",
-  ylab="y"
+**** Wave height                                                :noexport:
+:PROPERTIES:
+:header-args:R: :results output org
+:END:
+
+#+header:
+#+begin_src R :results output org
+source(file.path("R", "nonlinear.R"))
+propagating <- arma.print_wave_height(file.path("build", "nit-propagating"))
+standing <- arma.print_wave_height(file.path("build", "nit-standing"))
+result <- data.frame(
+  h1=c(propagating$h1, standing$h1),
+  h2=c(propagating$h2, standing$h2),
+  h3=c(propagating$h3, standing$h3)
 )
+rownames(result) <- c('propagating', 'standing')
+colnames(result) <- c('none', 'gcs', 'sn')
+ascii(result)
 #+end_src
 
-#+caption[Quantile-quantile plots for propagating waves]:
-#+caption: Quantile-quantile plots for propagating waves.
-#+name: propagating-wave-distributions
-#+RESULTS: propagating-wave-distributions
-[[file:build/propagating-wave-qqplots.pdf]]
+#+RESULTS:
+#+BEGIN_SRC org
+|             | none |  gcs |   sn |
+|-------------+------+------+------|
+| propagating | 2.41 | 2.75 | 1.37 |
+| standing    | 1.73 | 1.96 | 0.94 |
+#+END_SRC
+
+*** Non-physical nature of ARMA model
+ARMA model, owing to its non-physical nature, does not have the notion of sea
+wave; it simulates wavy surface as a whole instead. Motions of individual waves
+and their shape are often rough, and the total number of waves can not be
+determined precisely. However, integral characteristics of wavy surface match
+the ones of real sea waves.
+
+Theoretically, sea waves themselves can be chosen as ACFs, the only
+pre-processing step is to make them decay exponentially. This may allow
+to generate waves of arbitrary profiles, and is one of the directions of future
+work.
+
+* Pressure field under discretely given wavy surface
+** Known pressure field determination formulae
+**** Small amplitude waves theory.
+In\nbsp{}cite:stab2012,degtyarev1998modelling,degtyarev1997analysis the authors
+propose a solution for inverse problem of hydrodynamics of potential flow in the
+framework of small-amplitude wave theory (under assumption that wave length is
+much larger than height: \(\lambda \gg h\)). In that case inverse problem is
+linear and reduces to Laplace equation with mixed boundary conditions, and
+equation of motion is solely used to determine pressures for calculated velocity
+potential derivatives. The assumption of small amplitudes means the slow decay
+of wind wave coherence function, i.e. small change of local wave number in time
+and space compared to the wavy surface elevation (\(z\) coordinate). This
+assumption allows to calculate elevation \(z\) derivative as \(\zeta_z=k\zeta\),
+where \(k\) is wave number. In two-dimensional case the solution is written
+explicitly as
+\begin{align}
+    \left.\frac{\partial\phi}{\partial x}\right|_{x,t}= &
+        -\frac{1}{\sqrt{1+\alpha^{2}}}e^{-I(x)}
+            \int\limits_{0}^x\frac{\partial\dot{\zeta}/\partial
+                z+\alpha\dot{\alpha}}{\sqrt{1+\alpha^{2}}}e^{I(x)}dx,\label{eq-old-sol-2d}\\
+    I(x)= & \int\limits_{0}^x\frac{\partial\alpha/\partial z}{1+\alpha^{2}}dx,\nonumber
+\end{align}
+where \(\alpha\) is wave slope. In three-dimensional case the solution is
+written in the form of elliptic partial differential equation (PDE):
+\begin{align*}
+    & \frac{\partial^2 \phi}{\partial x^2} \left( 1 + \alpha_x^2 \right) +
+    \frac{\partial^2 \phi}{\partial y^2} \left( 1 + \alpha_y^2 \right) +
+    2\alpha_x\alpha_y \frac{\partial^2 \phi}{\partial x \partial y} + \\
+    & \left(
+        \frac{\partial \alpha_x}{\partial z} +
+        \alpha_x \frac{\partial \alpha_x}{\partial x} +
+        \alpha_y \frac{\partial \alpha_x}{\partial y}
+    \right) \frac{\partial \phi}{\partial x} + \\
+    & \left(
+        \frac{\partial \alpha_y}{\partial z} +
+        \alpha_x \frac{\partial \alpha_y}{\partial x} +
+        \alpha_y \frac{\partial \alpha_y}{\partial y}
+    \right) \frac{\partial \phi}{\partial y} + \\
+    & \frac{\partial \dot{\zeta}}{\partial z} +
+    \alpha_x \dot{\alpha_x} + \alpha_y \dot{\alpha_y} = 0.
+\end{align*}
+The authors suggest transforming this equation to finite differences and solve
+it numerically.
+
+As will be shown in section\nbsp{}[[#sec-compare-formulae]],
+formula\nbsp{}eqref:eq-old-sol-2d diverges when attempted to calculate velocity
+field for large-amplitude waves, and this is the reason that it can not be used
+in conjunction with a model, that generates arbitrary-amplitude waves.
+
+**** Linearisation of boundary condition.
+:PROPERTIES:
+:CUSTOM_ID: linearisation
+:END:
+
+LH model allows to derive an explicit formula for velocity field by linearising
+kinematic boundary condition. Velocity potential formula is written as
+\begin{equation*}
+\phi(x,y,z,t) = \sum_n \frac{c_n g}{\omega_n}
+     e^{\sqrt{u_n^2+v_n^2} z}
+     \sin(u_n x + v_n y - \omega_n t + \epsilon_n).
+\end{equation*}
+This formula is differentiated to obtain velocity potential derivatives, which
+are plugged to dynamic boundary condition to determine pressure field.
+
+** Determining wave pressures for discretely given wavy surface
+Analytic solutions to boundary problems in classical equations are often used to
+study different properties of the solution, and for that purpose general
+solution formula is too difficult to study, as it contains integrals of unknown
+functions. Fourier method is one of the methods to find analytic solutions to
+PDE. It is based on Fourier transform, applying of which to some PDEs reduces
+them to algebraic equations, and the solution is written as inverse Fourier
+transform of some function (which may contain Fourier transforms of other
+functions). Since it is not possible to write analytic forms of these Fourier
+transforms in some cases, unique solutions are found and their behaviour is
+studied in different domains instead. At the same time, computing discrete
+Fourier transforms numerically is possible for any discretely defined function
+using FFT algorithms. These algorithms use symmetry of complex exponentials to
+decrease asymptotic complexity from \(\mathcal{O}(n^2)\) to
+\(\mathcal{O}(n\log_{2}n)\). So, even if general solution contains Fourier
+transforms of unknown functions, they still can be computed numerically, and FFT
+family of algorithms makes this approach efficient.
+
+Alternative approach to solve PDEs is to reduce them to difference equations,
+which are solved by constructing various numerical schemes. This approach leads
+to approximate solution, and asymptotic complexity of corresponding algorithms
+is comparable to that of FFT. For example, for stationary elliptic PDE an
+implicit numerical scheme is constructed; the scheme is computed by iterative
+method on each step of which a tridiagonal or five-diagonal system of algebraic
+equations is solved by Thomas algorithm. Asymptotic complexity of this approach
+is \(\mathcal{O}({n}{m})\), where \(n\) is the number of wavy surface grid
+points and \(m\) is the number of iterations.
+
+Despite their wide spread, iterative algorithms are inefficient on parallel
+computer architectures due to inevitable process synchronisation after each
+iteration; in particular, their mapping to co-processors may involve copying
+data in and out of the co-processor on each iteration, which negatively affects
+their performance. At the same time, high number of Fourier transforms in the
+solution is an advantage, rather than a disadvantage. First, solutions obtained
+by Fourier method are explicit, hence their implementations scale with the large
+number of parallel computer cores. Second, there are implementations of FFT
+optimised for different processor architectures as well as co-processors (GPU,
+MIC) which makes it easy to get high performance on any computing platform.
+These advantages substantiate the choice of Fourier method to obtain explicit
+solution to the problem of determining pressures under wavy sea surface.
+
+*** Two-dimensional velocity field
+:PROPERTIES:
+:CUSTOM_ID: sec-pressure-2d
+:END:
+
+**** Formula for infinite depth fluid.
+Two-dimensional Laplace equation with Robin boundary condition is written as
+\begin{align}
+    \label{eq-problem-2d}
+    & \phi_{xx}+\phi_{zz}=0,\\
+    & \zeta_t + \zeta_x\phi_x = \frac{\zeta_x}{\sqrt{1 + \zeta_x^2}} \phi_x - \phi_z, & \text{на }z=\zeta(x,t).\nonumber
+\end{align}
+Use Fourier method to solve this problem. Applying Fourier transform to both
+sides of the equation yields
+\begin{equation*}
+    -4 \pi^2 \left( u^2 + v^2 \right)
+    \FourierY{\phi(x,z)}{u,v} = 0,
+\end{equation*}
+hence \(v = \pm i u\). Hereinafter we use the following symmetric form of
+Fourier transform:
+\begin{equation*}
+    \FourierY{f(x,y)}{u,v} =
+    \iint\limits_{-\infty}^{\phantom{--}\infty}
+    f(x,y)
+    e^{-2\pi i (x u + y v)}
+    dx dy.
+\end{equation*}
+We seek solution in the form of inverse Fourier transform
+\(\phi(x,z)=\InverseFourierY{E(u,v)}{x,z}\). Plugging[fn::\(v={-i}{u}\) is not
+applicable because velocity potential must go to nought when depth goes to
+infinity.] \(v={i}{u}\) into the formula yields
+\begin{equation}
+    \label{eq-guessed-sol-2d}
+    \phi(x,z) = \InverseFourierY{e^{2\pi u z}E(u)}{x}.
+\end{equation}
+In order to make substitution \(z=\zeta(x,t)\) not interfere with Fourier
+transforms, we rewrite\nbsp{}eqref:eq-guessed-sol-2d as a convolution:
+\begin{equation*}
+    \phi(x,z)
+    =
+    \Fun{z}
+    \ast
+    \InverseFourierY{E(u)}{x},
+\end{equation*}
+where \(\Fun{z}\)\nbsp{}--- a function, form of which is defined in
+section\nbsp{}[[#sec-compute-delta]] and which satisfies equation
+\(\FourierY{\Fun{z}}{u}=e^{2\pi{u}{z}}\). Plugging formula \(\phi\) into the
+boundary condition yields
+\begin{equation*}
+    \zeta_t
+    =
+    \left( i f(x) - 1 \right)
+    \left[
+        \Fun{z}
+        \ast
+        \InverseFourierY{2\pi u E(u)}{x}
+    \right],
+\end{equation*}
+where \(f(x)={\zeta_x}/{\sqrt{1+\zeta_x^2}}-\zeta_x\). Applying Fourier transform
+to both sides of this equation yields formula for coefficients \(E\):
+\begin{equation*}
+    E(u) =
+    \frac{1}{2\pi u}
+    \frac{
+    \FourierY{\zeta_t / \left(i f(x) - 1\right)}{u}
+    }{
+    \FourierY{\Fun{z}}{u}
+    }
+\end{equation*}
+Finally, substituting \(z\) for \(\zeta(x,t)\) and plugging resulting equation
+into\nbsp{}eqref:eq-guessed-sol-2d yields formula for \(\phi(x,z)\):
+\begin{equation}
+    \label{eq-solution-2d}
+    \boxed{
+        \phi(x,z)
+        =
+        \InverseFourierY{
+            \frac{e^{2\pi u z}}{2\pi u}
+            \frac{
+            \FourierY{ \zeta_t / \left(i f(x) - 1\right) }{u}
+            }{
+            \FourierY{ \Fun{\zeta(x,t)} }{u}
+            }
+        }{x}.
+    }
+\end{equation}
+
+Multiplier \(e^{2\pi{u}{z}}/(2\pi{u})\) makes graph of a function to which
+Fourier transform is applied asymmetric with respect to \(OY\) axis. This makes
+it difficult to use FFT, because it expects periodic function which takes nought
+values on both ends of the interval. At the same time, using numerical
+integration instead of FFT is not faster than solving the initial system of
+equations with numerical schemes. This problem is alleviated by using
+formula\nbsp{}eqref:eq-solution-2d-full for finite depth fluid with knowingly
+large depth \(h\). This formula is derived in the following section.
+
+**** Formula for finite depth fluid.
+On the sea bottom vertical fluid velocity component equals nought,
+i.e.\nbsp{}\(\phi_z=0\) on \(z=-h\), where \(h\)\nbsp{}--- water depth. In this
+case equation \(v=-{i}{u}\), which came from Laplace equation, can not be
+neglected, hence the solution is sought in the following form:
+\begin{equation}
+    \phi(x,z)
+    =
+    \InverseFourierY{
+        \left( C_1 e^{2\pi u z} + C_2 e^{-2\pi u z} \right)
+        E(u)
+    }{x}.
+    \label{eq-guessed-sol-2d-full}
+\end{equation}
+Plugging \(\phi\) into the boundary condition on the sea bottom yields
+\begin{equation*}
+    C_1 e^{-2\pi u h} - C_2 e^{2\pi u h} = 0,
+\end{equation*}
+hence \(C_1=\frac{1}{2}C{e}^{2\pi{u}{h}}\) and
+\(C_2=-\frac{1}{2}C{e}^{-2\pi{u}{h}}\). Constant \(C\) may take arbitrary values
+here, because after plugging it becomes part of unknown coefficients \(E(u)\).
+Plugging formulae for \(C_1\) and \(C_2\)
+into\nbsp{}eqref:eq-guessed-sol-2d-full yields
+\begin{equation*}
+    \phi(x,z) = \InverseFourierY{ \Sinh{2\pi u (z+h)} E(u) }{x}.
+\end{equation*}
+Plugging \(\phi\) into the boundary condition on the free surface yields
+\begin{equation*}
+    \zeta_t = f(x) \InverseFourierY{ 2\pi i u \Sinh{2\pi u (z+h)} E(u) }{x}
+            - \InverseFourierY{ 2\pi u \SinhX{2\pi u (z+h)} E(u) }{x}.
+\end{equation*}
+Here \(\sinh\) and \(\cosh\) give similar results near free surface, and since this
+is the main area of interest in practical applications, we assume that
+\(\Sinh{2\pi{u}(z+h)}\approx\SinhX{2\pi{u}(z+h)}\). Performing analogous to the
+previous section transformations yields final formula for \(\phi(x,z)\):
+\begin{equation}
+\boxed{
+    \phi(x,z,t)
+    =
+  \InverseFourierY{
+        \frac{\Sinh{2\pi u (z+h)}}{2\pi u}
+        \frac{
+            \FourierY{ \zeta_t / \left(i f(x) - 1\right) }{u}
+        }{
+            \FourierY{ \FunSecond{\zeta(x,t)} }{u}
+        }
+    }{x},
+}
+    \label{eq-solution-2d-full}
+\end{equation}
+where \(\FunSecond{z}\) is a function, a form of which is defined in
+section\nbsp{}[[#sec-compute-delta]] and which satisfies equation
+\(\FourierY{\FunSecond{z}}{u}=\Sinh{2\pi{u}{z}}\).
+
+**** Reduction to the formulae from linear wave theory.
+Check the validity of derived formulae by substituting \(\zeta(x,t)\) with known
+analytic formula for plain waves. Symbolic computation of Fourier transforms in
+this section were performed in Mathematica\nbsp{}cite:mathematica10. In the
+framework of linear wave theory assume that waves have small amplitude compared
+to their lengths, which allows us to simplify initial system of
+equations\nbsp{}eqref:eq-problem-2d to
+\begin{align*}
+    & \phi_{xx}+\phi_{zz}=0,\\
+    & \zeta_t = -\phi_z & \text{на }z=\zeta(x,t),
+\end{align*}
+solution to which is written as
+\begin{equation*}
+    \phi(x,z,t)
+    =
+    -\InverseFourierY{
+        \frac{e^{2\pi u z}}{2\pi u}
+        \FourierY{\zeta_t}{u}
+    }{x}
+    .
+\end{equation*}
+Propagating wave profile is defined as \(\zeta(x,t)=A\cos(2\pi(kx-t))\).
+Plugging this formula into\nbsp{}eqref:eq-solution-2d yields
+\(\phi(x,z,t)=-\frac{A}{k}\sin(2\pi(kx-t))\Sinh{2\pi{k}{z}}\). In order to
+reduce it to the formula from linear wave theory, rewrite hyperbolic sine in
+exponential form and discard the term containing \(e^{-2\pi{k}{z}}\) as
+contradicting condition
+\(\phi\underset{z\rightarrow-\infty}{\longrightarrow}0\). Taking real part of
+the resulting formula yields
+\(\phi(x,z,t)=\frac{A}{k}e^{2\pi{k}{z}}\sin(2\pi(kx-t))\), which corresponds to
+the known formula from linear wave theory. Similarly, under small-amplitude
+waves assumption the formula for finite depth
+fluid\nbsp{}eqref:eq-solution-2d-full is reduced to
+\begin{equation*}
+    \phi(x,z,t)
+    =
+    -\InverseFourierY{
+        \frac{\Sinh{2\pi u (z+h)}}{2\pi u \Sinh{2\pi u h}}
+        \FourierY{\zeta_t}{u}
+    }{x}.
+\end{equation*}
+Substituting \(\zeta(x,t)\) with propagating plain wave profile formula yields
+\begin{equation}
+    \label{eq-solution-2d-linear}
+    \phi(x,z,t)=\frac{A}{k}
+    \frac{\Sinh{2 \pi k (z+h)}}{ \Sinh{2 \pi k h} }
+    \sin(2 \pi (k x-t)),
+\end{equation}
+which corresponds to the formula from linear wave theory for finite depth fluid.
+
+Different forms of Laplace equation solutions, in which decaying exponent is
+written with either "+" or "-" signs, may cause incompatibilities between
+formulae from linear wave theory and formulae derived in this work, where
+\(\sinh\) is used instead of \(\cosh\). Equality
+\(\frac{\Sinh{2\pi{k}(z+h)}}{\Sinh{2\pi{k}{h}}}\approx\frac{\sinh(2\pi{k}(z+h))}{\sinh(2\pi{k}{h})}\)
+becomes strict on the free surface, and difference between left-hand and
+right-hand sides increases when approaching sea bottom (for sufficiently large
+depth difference near free surface is negligible). So, for sufficiently large
+depth any function (\(\cosh\) or \(\sinh\)) may be used for velocity potential
+computation near free surface.
+
+Reducing\nbsp{}eqref:eq-solution-2d and\nbsp{}eqref:eq-solution-2d-full to the
+known formulae from linear wave theory shows, that formula for infinite
+depth\nbsp{}eqref:eq-solution-2d is not suitable to compute velocity potentials
+with Fourier method, because it does not have symmetry, which is required for
+Fourier transform. However, formula for finite depth can be used instead by
+setting \(h\) to some characteristic water depth. For standing wave reducing to
+linear wave theory formulae is made under the same assumptions.
+*** Three-dimensional velocity field
+Three-dimensional version of\nbsp{}eqref:eq-problem is written as
+\begin{align}
+    \label{eq-problem-3d}
+    & \phi_{xx} + \phi_{yy} + \phi_{zz} = 0,\\
+    & \zeta_t + \zeta_x\phi_x + \zeta_y\phi_y
+    =
+    \frac{\zeta_x}{\SqrtZeta{1 + \zeta_x^2 + \zeta_y^2}} \phi_x
+    +\frac{\zeta_y}{\SqrtZeta{1 + \zeta_x^2 + \zeta_y^2}} \phi_y
+    - \phi_z, & \text{на }z=\zeta(x,y,t).\nonumber
+\end{align}
+Again, use Fourier method to solve it. Applying Fourier transform to both sides
+of Laplace equation yields
+\begin{equation*}
+    -4 \pi^2 \left( u^2 + v^2 + w^2 \right)
+    \FourierY{\phi(x,y,z)}{u,v,w} = 0,
+\end{equation*}
+hence \(w=\pm{i}\sqrt{u^2+v^2}\). We seek solution in the form of inverse Fourier
+transform \(\phi(x,y,z)=\InverseFourierY{E(u,v,w)}{x,y,z}\). Plugging
+\(w=i\sqrt{u^2+v^2}=i\Kveclen\) into the formula yields
+\begin{equation*}
+    \phi(x,y,z) = \InverseFourierY{
+        \left(
+            C_1 e^{2\pi \Kveclen z}
+            -C_2 e^{-2\pi \Kveclen z}
+        \right)
+        E(u,v)
+    }{x,y}.
+\end{equation*}
+Plugging \(\phi\) into the boundary condition on the sea bottom (analogous to
+two-dimensional case) yields
+\begin{equation}
+    \label{eq-guessed-sol-3d}
+    \phi(x,y,z) = \InverseFourierY{
+        \Sinh{2\pi \Kveclen (z+h)} E(u,v)
+    }{x,y}.
+\end{equation}
+Plugging \(\phi\) into the boundary condition on the free surface yields
+\begin{equation*}
+    \arraycolsep=1.4pt
+    \begin{array}{rl}
+        \zeta_t = & i f_1(x,y) \InverseFourierY{2 \pi u \Sinh{2\pi \Kveclen (z+h)}E(u,v)}{x,y} \\
+        + & i f_2(x,y) \InverseFourierY{2 \pi v \Sinh{2\pi \Kveclen (z+h)}E(u,v)}{x,y} \\
+        - & \InverseFourierY{2 \pi \Kveclen \SinhX{2\pi \Kveclen (z+h)}E(u,v)}{x,y}
+    \end{array}
+\end{equation*}
+where \(f_1(x,y)={\zeta_x}/{\SqrtZeta{1+\zeta_x^2+\zeta_y^2}}-\zeta_x\) and
+\(f_2(x,y)={\zeta_y}/{\SqrtZeta{1+\zeta_x^2+\zeta_y^2}}-\zeta_y\).
 
-#+name: standing-wave-distributions
-#+begin_src R :file build/standing-wave-qqplots.pdf
-source(file.path("R", "common.R"))
-par(pty="s", mfrow=c(2, 2))
-arma.qqplot_grid(
-  file.path("build", "standing_wave"),
-  c("elevation", "heights_y", "lengths_y", "periods"),
-  c("elevation", "height Y", "length Y", "period"),
-  xlab="x",
-  ylab="y"
-)
-#+end_src
+Like in section\nbsp{}[[#sec-pressure-2d]] we assume that
+\(\Sinh{2\pi{u}(z+h)}\approx\SinhX{2\pi{u}(z+h)}\) near free surface, but in
+three-dimensional case this is not enough to solve the problem. In order to get
+explicit formula for coefficients \(E\) we need to assume, that all Fourier
+transforms in the equation have radially symmetric kernels, i.e.\nbsp{}replace
+\(u\) and \(v\) with \(\Kveclen\). There are two points supporting this
+assumption. First, in numerical implementation integration is done over positive
+wave numbers, so the sign of \(u\) and \(v\) does not affect the solution.
+Second, the growth rate of \(\cosh\) term of the integral kernel is much higher
+than the one of \(u\) or \(\Kveclen\), so the substitution has small effect on
+the magnitude of the solution. Despite these two points, a use of more
+mathematically rigorous approach would be preferable.
 
-#+caption[Quantile-quantile plots for standing waves]:
-#+caption: Quantile-quantile plots for standing waves.
-#+name: standing-wave-distributions
-#+RESULTS: standing-wave-distributions
-[[file:build/standing-wave-qqplots.pdf]]
+Making the replacement, applying Fourier transform to both sides of the equation
+and plugging the result into\nbsp{}eqref:eq-guessed-sol-3d yields formula for
+\(\phi\):
+\begin{equation}
+    \label{eq-phi-3d}
+    \phi(x,y,z,t) = \InverseFourierY{
+        \frac{ \Sinh{\smash{2\pi \Kveclen (z+h)}} }{ 2\pi\Kveclen }
+        \frac{ \FourierY{ \zeta_t / \left( i f_1(x,y) + i f_2(x,y) - 1 \right)}{u,v} }
+        { \FourierY{\mathcal{D}_3\left( x,y,\zeta\left(x,y\right) \right)}{u,v} }
+    }{x,y},
+\end{equation}
+where
+\(\FourierY{\mathcal{D}_3\left(x,y,z\right)}{u,v}=\Sinh{\smash{2\pi\Kveclen{}z}}\).
 
-#+name: acf-slices
-#+header: :width 6 :height 9
-#+begin_src R :file build/acf-slices.pdf
-source(file.path("R", "common.R"))
-propagating_acf <- read.csv(file.path("build", "propagating_wave", "acf.csv"))
-standing_acf <- read.csv(file.path("build", "standing_wave", "acf.csv"))
-par(mfrow=c(5, 2), mar=c(0,0,0,0))
-for (i in seq(0, 4)) {
-  arma.wavy_plot(standing_acf, i, zlim=c(-5,5))
-  arma.wavy_plot(propagating_acf, i, zlim=c(-5,5))
-}
-#+end_src
+*** Velocity potential normalisation formulae
+:PROPERTIES:
+:CUSTOM_ID: sec-compute-delta
+:END:
 
-#+caption[ACF time slices for standing and propagating waves]:
-#+caption: ACF time slices for standing (left column) and propagating waves
-#+caption: (right column).
-#+name: acf-slices
-#+RESULTS: acf-slices
-[[file:build/acf-slices.pdf]]
+In solutions\nbsp{}eqref:eq-solution-2d and\nbsp{}eqref:eq-solution-2d-full to
+two-dimensional pressure determination problem there are functions
+\(\Fun{z}=\InverseFourierY{e^{2\pi{u}{z}}}{x}\) and
+\(\FunSecond{z}=\InverseFourierY{\Sinh{2\pi{u}{z}}}{x}\) which has multiple
+analytic representations and are difficult to compute. Each function is a
+Fourier transform of linear combination of exponents which reduces to poorly
+defined Dirac delta function of a complex argument (see
+table\nbsp{}[[tab-delta-functions]]). The usual way of handling this type of
+functions is to write them as multiplication of Dirac delta functions of real
+and imaginary part, however, this approach does not work here, because applying
+inverse Fourier transform to this representation does not produce exponent,
+which severely warp resulting velocity field. In order to get unique analytic
+definition normalisation factor \(1/\Sinh{2\pi{u}{h}}\) (which is also included
+in formula for \(E(u)\)) may be used. Despite the fact that normalisation allows
+to obtain adequate velocity potential field, numerical experiments show that
+there is little difference between this field and the one produced by formulae
+from linear wave theory, in which terms with \(\zeta\) are omitted. As a result,
+the formula for three-dimensional case was not derived.
 
-Graph tails in fig.\nbsp{}[[propagating-wave-distributions]] deviate from original
-distribution for individual wave characteristics, because every wave have to be
-extracted from the resulting wavy surface to measure its length, period and
-height. There is no algorithm that guarantees correct extraction of all waves,
-because they may and often overlap each other. Weibull distribution right tail
-represents infrequently occurring waves, so it deviates more than left tail.
+#+name: tab-delta-functions
+#+caption[Formulae for computing \(\Fun{z}\) and \(\FunSecond{z}\)]:
+#+caption: Formulae for computing \(\Fun{z}\) and \(\FunSecond{z}\) from
+#+caption: sec.\nbsp{}[[#sec-pressure-2d]], that use normalisation to eliminate
+#+caption: uncertainty from definition of Dirac delta function of complex argument.
+#+attr_latex: :booktabs t
+| Function          | Without normalisation                                        | Normalised                                                                                                                             |
+|-------------------+--------------------------------------------------------------+----------------------------------------------------------------------------------------------------------------------------------------|
+| \(\Fun{z}\)       | \(\delta (x+i z)\)                                           | \(\frac{1}{2 h}\mathrm{sech}\left(\frac{\pi  (x-i (h+z))}{2 h}\right)\)                                                                |
+| \(\FunSecond{z}\) | \(\frac{1}{2}\left[\delta (x-i z) + \delta (x+i z) \right]\) | \(\frac{1}{4 h}\left[\text{sech}\left(\frac{\pi  (x-i (h+z))}{2 h}\right)+\text{sech}\left(\frac{\pi  (x+i(h+z))}{2 h}\right)\right]\) |
 
-Correspondence rate for standing waves (fig.\nbsp{}[[standing-wave-distributions]])
-is lower for height and length, roughly the same for surface elevation and
-higher for wave period distribution tails. Lower correspondence degree for
-length and height may be attributed to the fact that Weibull distributions were
-obtained empirically for sea waves which are typically propagating, and
-distributions may be different for standings waves. Higher correspondence degree
-for wave periods is attributed to the fact that wave periods of standing waves
-are extracted more precisely as the waves do not move outside simulated wavy
-surface region. The same correspondence degree for wave elevation is obtained,
-because this is the characteristic of the wavy surface (and corresponding AR or
-MA process) and is not affected by the type of waves.
+*** Conclusion
+Formulae derived in this section allow to numerically compute velocity potential
+field (and hence pressure field) near discretely or mathematically given wavy
+sea surface, bypassing assumptions of linear wavy theory and theory of small
+amplitude waves.
 
-*** Verification of velocity potential fields
+** Verification of velocity potential fields
 :PROPERTIES:
 :CUSTOM_ID: sec-compare-formulae
 :END:
@@ -1557,7 +1580,7 @@ known, even for plain waves, so comparison is done numerically. Taking into
 account conclusions of section\nbsp{}[[#sec-pressure-2d]], only finite depth
 formulae are compared.
 
-**** The difference with linear wave theory formulae.
+*** The difference with linear wave theory formulae.
 In order to obtain velocity potential fields, sea wavy surface was generated
 by AR model with varying wave amplitude. In numerical implementation wave
 numbers in Fourier transforms were chosen on the interval from \(0\) to the
@@ -1627,7 +1650,7 @@ arma.plot_velocity_potential_field_legend(
 #+RESULTS: fig-potential-field-nonlinear
 [[file:build/plain-wave-velocity-field-comparison.pdf]]
 
-**** The difference with small-amplitude wave theory.
+*** The difference with small-amplitude wave theory.
 The experiment, in which velocity fields produced numerically by different
 formulae were compared, shows that velocity fields produced by
 formula\nbsp{}eqref:eq-solution-2d-full and\nbsp{}eqref:eq-old-sol-2d correspond
@@ -1672,113 +1695,6 @@ arma.plot_velocity(
 #+RESULTS: fig-velocity-field-2d
 [[file:build/large-and-small-amplitude-velocity-field-comparison.pdf]]
 
-*** Verification of nonlinear inertialess transformation
-In order to measure the effect of NIT on the shape of the resulting wavy
-surface, three realisations were generated:
-- realisation with Gaussian distribution (without NIT),
-- realisation with Gram---Charlier series (GCS) based distribution, and
-- realisation with skew normal distribution.
-The seed of PRNG was set to be the same for all programme executions to make ARMA
-model produce the same values for each realisation. There we two experiments:
-for standing and propagating waves with ACFs given by formulae from
-section\nbsp{}[[#sec-wave-acfs]].
-
-While the experiment showed that applying NIT with GCS-based distribution
-increases wave steepness, the same is not true for skew normal distribution
-(fig.\nbsp{}[[fig-nit]]). Using this distribution results in wavy surface each
-\(z\)-coordinate of which is always greater or equal to nought. So, skew normal
-distribution is unsuitable for NIT. NIT increases the wave height and steepness
-of both standing and propagating waves. Increasing either skewness or kurtosis
-of GCS-based distribution increases both wave height and steepness. The error of
-ACF approximation (eq.\nbsp{}eqref:eq-nit-error) ranges from 0.20 for GCS-based
-distribution to 0.70 for skew normal distribution (table\nbsp{}[[tab-nit-error]]).
-
-#+name: fig-nit
-#+header: :width 7 :height 7
-#+begin_src R :file build/nit.pdf
-source(file.path("R", "nonlinear.R"))
-par(mfrow=c(2, 1), mar=c(4,4,4,0.5), family='serif')
-args <- list(
-  graphs=c('Gaussian', 'Gram—Charlier', 'Skew normal'),
-  linetypes=c('solid', 'dashed', 'dotted')
-)
-args$title <- 'Propagating waves'
-arma.plot_nonlinear(file.path("build", "nit-propagating"), args)
-args$title <- 'Standing waves'
-arma.plot_nonlinear(file.path("build", "nit-standing"), args)
-#+end_src
-
-#+name: fig-nit
-#+caption[Wavy surface slices with different distributions]:
-#+caption: Wavy surface slices with different distributions
-#+caption: of wave elevation (Gaussian, GCS-based and SN).
-#+RESULTS: fig-nit
-[[file:build/nit.pdf]]
-
-#+name: tab-nit-error
-#+caption[Errors of ACF approximations]:
-#+caption: Errors of ACF approximations (the difference of variances) for
-#+caption: different wave elevation distributions. \(N\) is the number of
-#+caption: coefficients of ACF approximation.
-#+attr_latex: :booktabs t
-| Wave type   | Distribution | \(\gamma_1\) | \(\gamma_2\) | \(\alpha\) | Error | \(N\) | Wave height |
-|-------------+--------------+--------------+--------------+------------+-------+-------+-------------|
-| propagating | Gaussian     |              |              |            |       |       |        2.41 |
-| propagating | GCS-based    |         2.25 |          0.4 |            |  0.20 |     2 |        2.75 |
-| propagating | skew normal  |              |              |          1 |  0.70 |     3 |        1.37 |
-| standing    | Gaussian     |              |              |            |       |       |        1.73 |
-| standing    | GCS-based    |         2.25 |          0.4 |            |  0.26 |     2 |        1.96 |
-| standing    | skew normal  |              |              |          1 |  0.70 |     3 |        0.94 |
-
-To summarise, the only test case that showed acceptable results is realisation
-with GCS-based distribution for both standing and propagating waves. Skew normal
-distribution warps wavy surface for both types of waves. GCS-based realisations
-have large error of ACF approximation, which results in increase of wave height.
-The reason for the large error is that GCS approximations are not accurate as
-they do not converge for all possible
-functions\nbsp{}cite:wallace1958asymptotic. Despite the large error, the change
-in wave height is small (table\nbsp{}[[tab-nit-error]]).
-
-**** Wave height                                                :noexport:
-:PROPERTIES:
-:header-args:R: :results output org
-:END:
-
-#+header: 
-#+begin_src R :results output org
-source(file.path("R", "nonlinear.R"))
-propagating <- arma.print_wave_height(file.path("build", "nit-propagating"))
-standing <- arma.print_wave_height(file.path("build", "nit-standing"))
-result <- data.frame(
-  h1=c(propagating$h1, standing$h1),
-  h2=c(propagating$h2, standing$h2),
-  h3=c(propagating$h3, standing$h3)
-)
-rownames(result) <- c('propagating', 'standing')
-colnames(result) <- c('none', 'gcs', 'sn')
-ascii(result)
-#+end_src
-
-#+RESULTS:
-#+BEGIN_SRC org
-|             | none |  gcs |   sn |
-|-------------+------+------+------|
-| propagating | 2.41 | 2.75 | 1.37 |
-| standing    | 1.73 | 1.96 | 0.94 |
-#+END_SRC
-
-*** Non-physical nature of ARMA model
-ARMA model, owing to its non-physical nature, does not have the notion of sea
-wave; it simulates wavy surface as a whole instead. Motions of individual waves
-and their shape are often rough, and the total number of waves can not be
-determined precisely. However, integral characteristics of wavy surface match
-the ones of real sea waves.
-
-Theoretically, sea waves themselves can be chosen as ACFs, the only
-pre-processing step is to make them decay exponentially. This may allow
-to generate waves of arbitrary profiles, and is one of the directions of future
-work.
-
 * High-performance software implementation of sea wave simulation
 ** SMP implementation
 **** Parallel AR, MA and LH model algorithms.
@@ -1862,8 +1778,9 @@ previously computed Fourier transform of the coefficients, and inverse Fourier
 transform of the result is computed. After that, each part is written to the
 output array with overlapping (due to padding) points added to each other. This
 algorithm is commonly known in signal processing as
-"overlap-add"\nbsp{}cite:svoboda2011efficient. Padding with noughts is needed to
-prevent aliasing errors: without it the result would be circular convolution.
+"overlap-add"\nbsp{}cite:oppenheim1989discrete,svoboda2011efficient,pavel2013algorithms.
+Padding with noughts is needed to prevent aliasing errors: without it the result
+would be circular convolution.
 
 Despite the fact that MA model algorithm partitions the surface into the same
 parts (but possibly of different sizes) as AR model algorithm, the vicinity of
@@ -2302,6 +2219,52 @@ disadvantage of using OpenCL and OpenGL together is the requirement for manual
 locking of shared buffer: failure to do so results in appearance of screen image
 artefacts which can be removed only by rebooting the computer.
 
+**** Parallel white noise generation algorithm.
+In order to eliminate periodicity from wavy surface generated by sea wave model,
+it is imperative to use PRNG with sufficiently large period to generate white
+noise. Parallel Mersenne Twister\nbsp{}cite:matsumoto1998mersenne with a period
+of \(2^{19937}-1\) is used as a generator in this work. It allows to produce
+aperiodic sea wavy surface realisations in any practical usage scenarios.
+
+There is no guarantee that multiple PRNGs executed in parallel threads with
+distinct initial states produce uncorrelated pseudo-random number sequences, and
+algorithm of dynamic creation of Mersenne
+Twisters\nbsp{}cite:matsumoto1998dynamic is used to provide such guarantee. The
+essence of the algorithm is to find matrices of initial generator states, that
+give maximally uncorrelated pseudo-random number sequences when Mersenne
+Twisters are executed in parallel with these initial states. Since finding such
+initial states consumes considerable amount of processor time, vector of initial
+states is created beforehand with knowingly larger number of parallel threads
+and saved to a file, which is then read before starting white noise generation.
+
+**** Ramp-up interval elimination.
+In AR model value of wavy surface elevation at a particular point depends on
+previous in space and time points, as a result the so called /ramp-up interval/
+(see fig.\nbsp{}[[fig-ramp-up-interval]]), in which realisation does not correspond
+to specified ACF, forms in the beginning of the realisation. There are several
+solutions to this problem which depend on the simulation context. If realisation
+is used in the context of ship stability simulation without manoeuvring, ramp-up
+interval will not affect results of the simulation, because it is located on the
+border (too far away from the studied marine object). Alternative approach is to
+generate sea wavy surface on ramp-up interval with LH model and generate the
+rest of the realisation with AR model. If ship stability with manoeuvring is
+studied, then the interval may be simply discarded from the realisation (the
+size of the interval approximately equals the number of AR coefficients in each
+dimension). However, this may lead to loss of a very large number of points,
+because discarding occurs for each of three dimensions.
+
+#+name: fig-ramp-up-interval
+#+begin_src R :file build/ramp-up-interval.pdf
+source(file.path("R", "common.R"))
+arma.plot_ramp_up_interval()
+#+end_src
+
+#+caption[Ramp-up interval at the beginning of AR process realisation]:
+#+caption: Ramp-up interval at the beginning of AR process realisation.
+#+name: fig-ramp-up-interval
+#+RESULTS: fig-ramp-up-interval
+[[file:build/ramp-up-interval.pdf]]
+
 **** Conclusions.
 Benchmarks showed that GPU outperforms CPU in arithmetic intensive tasks,
 i.e.\nbsp{}tasks requiring high number of floating point operations per second,