\documentclass[12pt]{article}
\usepackage{graphicx}
\usepackage{color}
\usepackage[utf8]{inputenc}
\usepackage[T2A]{fontenc}
\usepackage[russian]{babel}
\usepackage{amsmath}
\usepackage{amssymb}

\topmargin=-1.8cm
\oddsidemargin=-15mm
\evensidemargin=-15mm
\textheight=24.5cm
\textwidth=19cm
\tolerance=1000

\parskip=5pt plus 4pt minus 2pt
\tolerance=9000
%
\title{Интерполяция полей скорости}
\begin{document}
\maketitle
\section{Координаты и уравнения движения}
Пусть у нас задано поле скорости $(u_{i,j},v_{i,j})$ на некоторой двумерной сетке $(\lambda_{i,j},\varphi_{i,j})$, $i=0,1,\dots,N_x$, $j=0,1,\dots,N_y$, $\lambda_{i,j}$~--- долгота, $\varphi_{i,j}$~--- широта. Мы предполагаем, что эта сетка может быть преобразована к прямоугольному виду, то есть существуют такие функции $X(\lambda,\varphi)$ и $Y(\lambda,\varphi)$, что $x_{i,j}=X(\lambda_{i,j},\varphi_{i,j})=i\Delta x$ и $y_{i,j}=Y(\lambda_{i,j},\varphi_{i,j})=j\Delta y$, где $\Delta x$, $\Delta y$~--- некоторые постоянные числа. Проще говоря, существуют координаты, в которых наша сетка имеет простой прямоугольный вид. В самом простом случае сетка $(\lambda_{i,j},\varphi_{i,j})$ уже имеет прямоугольный вид. Но, для глобальных полей такую сетку не используют в силу сингулярности на Северном полюсе и зависимости разрешения от широты. Сами скорости задаются обычно в см/сек. Однако, для описания крупномасштабных процессов это неудобно и лучше использовать скорости, выраженные в км/сут
%
\begin{equation}
1\,\frac{\text{см}}{\text{сек}}=\frac{10^{-5}\,\text{км}}{\frac{1}{24*3600}\,\text{сут}}=0{,}864\frac{\text{км}}{\text{сут}}.
\end{equation}
%

Есть несколько возможных подходов к написанию лагранжевых уравнений адвекции.
%
\begin{enumerate}
\item Не преобразовывать ни скорости, ни координаты. Уравнения адвекции выглядят следующим образом:
%
\begin{equation}
\dot\lambda=\frac{10800\cos{\varphi}}{\pi R}u,\quad \dot\varphi=\frac{10800}{\pi R}v,
\end{equation}
%
где $R$~--- радиус Земли. Широта и долготы выражены в минутах, время~--- в сутках, $u$ и $v$~--- в км/сут. Коэффициент $\frac{\pi R}{10800}$ примерно равен морской миле~--- $1{,}852$ км.\\
{\em Плюсы:} нет необходимости в каких-либо преобразованиях координат, входные и выходные данные задаются непосредственно в географических координатах.\\
{\em Минусы:} вычисление косинуса на каждом шаге~--- затратная операция, в случае непрямоугольной сетки определение ячейки, в которой находятся текущие координаты (необходимо для нахождения интерполированых скоростей), является непростой операцией.
\item Перейдя от линейных скоростей в правой части к угловым, мы получаем максимально простые уравнения адвекции
%
\begin{equation}
\dot\lambda=u,\quad \dot\varphi=v.
\label{maineq}
\end{equation}
%
Здесь $u$ и $v$ выражены в мин/сут.\\
{\em Плюсы:} те же, что и у первого способа, быстрое численное интегрирование уравнений.\\
{\em Минусы:} всё те же проблемы с непрямоугольными сетками.
\item Преобразовать координаты и скорости в прямоугольную сетку.
Уравнения движения будут иметь тот же вид, что и в предыдущем случае, но для некоторых абстрактных координат и скоростей.
%
\begin{equation}
\dot x=u,\quad \dot y=v.
\end{equation}
%
{\em Плюсы:} максимально простое и быстрое численное интегрирование.\\
{\em Минусы:} мы считаем в абстрактных координатах, следовательно, все входные и выходные данные должны преобразовываться. Кроме этого, при данном подходе неудобно рассчитывать различные интегральные функции от траектории (например, её длину) из-за необходимости преобразования метрики на каждом шаге, что может сильно увеличить время счёта.
\end{enumerate}
%
Мы считаем последним способом ещё и потому, что в этом случае програмный код получается универсальным и не зависит от реальной сетки $(\lambda_{ij},\varphi_{ij})$.

\section{Интерполяция}
Поле скорости у нас задано только на сетке $(x_{i,j},y_{i,j},t_k)$, однако, для интегрирования уравнений адвекции необходимы значения скорости в произвольной точке пространства $(x,y)$ и в произвольный момент времени $t$. Для нахождения этих значений используются различные методы интерполяции. Мы используем два: линейную и бикубическую интерполяцию в пространстве и линейную и кубическими полиномами третьего порядка во времени.
\subsection{Интерполяция в пространстве}
Рассмотрим отдельную ячейку нашей пространственной сетки единичной ширины и высоты (это всегда можно сделать соответствующей нормировкой). В четырёх углах углах ячейки нам известны значения некоторой функции $z(x,y)$, а именно $z(0,0)=z_{00}$, $z(1,0)=z_{10}$, $z(0,1)=z_{01}$, $z(1,1)=z_{11}$ (см. Рис.~\ref{Fig-cell})
%
\begin{figure}[!htb]
\begin{center}
\includegraphics[width=0.4\textwidth]{net.eps}
\caption{Элементарная ячейка сетки, в которой производится интерполяция.}
\label{Fig-cell}
\end{center}
\end{figure}
%
Рассмотрим сперва линейную интерполяцию. В этом случае мы считаем, что функция $z(x,y)$ является линейной вдоль любой вертикальной или горизонтальной линии. Соответственно, мы можем найти $z(x,y)$, найдя сперва $z(x,0)$, $z(x,1)$ (или $z(0,y)$, $z(1,y)$, результат будет тот же)
%
\begin{equation}
\begin{gathered}
z(x,0)=(z(1,0)-z(0,0))x+z(0,0)=(z_{10}-z_{00})x+z_{00},\\
z(x,1)=(z(1,1)-z(0,1))x+z(0,1)=(z_{11}-z_{01})x+z_{01},\\
z(x,y)=(z(x,1)-z(x,0))y+z(x,0)=axy+bx+cy+d,\\
a=z_{11}-z_{01}-z_{10}+z_{00},\quad b=z_{10}-z_{00},\quad c=z_{01}-z_{00},\quad d=z_{00}.
\end{gathered}
\label{linint}
\end{equation}
%
Таким образом, несмотря на название метода, интерполирующая функция является нелинейной. Линейная интерполяция обеспечивает непрерывность функции $z(x,y)$ при пересечении границ ячейки, но её производные разрывны на этих границах.

Интерполирующая функция для бикубической интерполяции имеет вид
%
\begin{multline}
z(x,y)=a_{00}+a_{10}x+a_{20}x^2+a_{30}x^3+a_{01}y+a_{11}xy+a_{21}x^2y+a_{31}x^3y+\\
+a_{02}y^2+a_{12}xy^2+a_{22}x^2y^2+a_{32}x^3y^2+a_{03}y^3+a_{13}xy^3+a_{23}x^2y^3+a_{33}x^3y^3
\label{bicubic}
\end{multline}
%
и содержит $16$ коэффициентов. Для их нахождения необходимо знать не только значения функции на сетке, но и значения производных $\partial z/\partial x$, $\partial z/\partial y$ и $\partial^2 z/\partial x \partial y$. Так как мы имеем только сами значения функции, нам необходимо найти эти производные путём численного дифференцирования
%
\begin{equation}
\begin{aligned}
\frac{\partial z}{\partial x}&=z^x_{i,j}=\frac{z_{i+1,j}-z_{i-1,j}}{2},\\
\frac{\partial z}{\partial y}&=z^y_{i,j}=\frac{z_{i,j+1}-z_{i,j-1}}{2},\\
\frac{\partial^2 z}{\partial x\partial y}&=z^{xy}_{i,j}=\frac{z^y_{i+1,j}-z^y_{i-1,j}}{2}.
\end{aligned}
\end{equation}
%
В точках, где использование центральных разностей невозможно (вблизи границ сетки или около берега), используются правые или левые разности. Таким образом, для каждой ячейки сетки мы получаем $16$  значений функции, её первых производных и смешаной второй производной в углах ячейки. Коэффициенты интерполирующей функции получаются следующим образом
%
\begin{equation}
\begin{gathered}
a_{00}=z_{00},\quad
a_{01}=z^y_{00},\quad
a_{02}=3(z_{01}-z_{00})-2z^y_{00}-z^y_{01},\quad
a_{03}=2(z_{00}-z_{01})+z^y_{00}+z^y_{01},\\
a_{10}=z^x_{00},\quad
a_{11}=z^{xy}_{00},\quad
a_{12}=3(z^x_{01}-z^x_{00})-2z^{xy}_{00}-z^{xy}_{01},\quad
a_{13}=2(z^x_{00}-z^x_{01})+z^{xy}_{00}+z^{xy}_{01},\\
a_{20}=3(z_{10}-z_{00})-2z^x_{00}-z^x_{10},\quad
a_{21}=3(z^y_{10}-z^y_{00})-2z^{xy}_{00}-z^{xy}_{10},\\
\begin{aligned}
a_{22}=9(z_{00}-z_{01}-z_{10}+z_{11})&+6(z^x_{00}-z^x_{01}+z^y_{00}-z^y_{10})+{}\\
{}&+3(z^x_{10}-z^x_{11}+z^y_{01}-z^y_{11})+2(z^{xy}_{01}+z^{xy}_{10})+4z^{xy}_{00}+z^{xy}_{11},
\end{aligned}\\
\begin{aligned}
a_{23}=6(z_{01}+z_{10}-z_{00}-z_{11})+4(z^x_{01}-z^x_{00})&+3(z^y_{10}+z^y_{11}-z^y_{00}-z^y_{01})+{}\\
{}&+2(z^x_{11}-z^x_{10}-z^{xy}_{00}-z^{xy}_{01})-z^{xy}_{10}-z^{xy}_{11},
\end{aligned}\\
a_{30}=2(z_{00}-z_{10})+z^x_{00}+z^x_{10},\quad
a_{31}=2(z^y_{00}-z^y_{10})+z^{xy}_{00}+z^{xy}_{10},\\
\begin{aligned}
a_{32}=6(z_{01}+z_{10}-z_{00}-z_{11})+4(z^y_{10}-z^y_{00})&+3(z^x_{01}+z^x_{11}-z^x_{00}-z^x_{10})+{}\\
{}&+2(z^y_{11}-z^y_{01}-z^{xy}_{00}-z^{xy}_{10})-z^{xy}_{01}-z^{xy}_{11},\\
\end{aligned}\\
\begin{aligned}
a_{33}=4(z_{00}-z_{01}-z_{10}+z_{11})+2(z^x_{00}-z^x_{01}+z^x_{10}-z^x_{11}+z^y_{00}&-z^y_{10}+z^y_{01}-z^y_{11})+{}\\
{}&+z^{xy}_{00}+z^{xy}_{01}+z^{xy}_{10}+z^{xy}_{11}.
\end{aligned}
\end{gathered}
\label{bicub_coeff}
\end{equation}
%
Из-за более сложных формул, бикубическая интерполяция работает существенно медленнее линейной.

\subsection{Интерполяция во времени}
Так как нам нужно получать значения скоростей не только в произвольном месте, но и в произвольное время, необходимо делать дополнительно интерполяцию во времени. Эта задача существенно проще пространственной интерполяции, так как в данном случае мы имеем дело с одномерной функцией. Простейший вид интерполяции~--- линейный
%
\begin{equation}
z(x,y,t)=\frac{z(x,y,t_{k+1})-z(x,y,t_k)}{\Delta t} (t-t_k)+z(x,y,t_k).
\label{lint}
\end{equation}
%
где $\Delta t$~--- шаг сетки по времени.

Мы также используем интерполяцию кубическими полиномами Лагранжа\footnote{Полином Лагранжа степени $n$~--- это такой полином, который в $n+1$ точке имеет заданные значения.}. Для $t\in[t_k:t_{k+1}]$ мы имеем
%
\begin{equation}
z(x,y,t)=a_0(x,y)+a_1(x,y)t+a_2(x,y)t^2+a_3(x,y)t^3.
\label{cubt}
\end{equation}
%
Коэффициенты $a_i(x,y)$ рассчитываются из значений $z_i=z(x,y,t_i)$, которые, в свою очередь, находятся пространственной интерполяцией
%
\begin{equation}
\begin{gathered}
a_0=z_k,\qquad
a_1=\left\{
 \begin{aligned}
 -&\frac{11}{6}z_k+3z_{k+1}-\frac{3}{2}z_{k+2}+\frac{1}{3}z_{k+3}, &k&=0,\\
 -&\frac{1}{3}z_{k-1}-\frac{1}{2}z_k+z_{k+1}-\frac{1}{6}z_{k+2},   &k&\in[1:N_t-2],\\
  &\frac{1}{6}z_{k-2}-z_{k-1}+\frac{1}{2}z_k+\frac{1}{3}z_{k+1},   &k&=N_t-1,
 \end{aligned}\right.\\
a_2=\left\{
 \begin{aligned}
  &z_k-\frac{5}{2}z_{k+1}+2z_{k+2}-\frac{1}{2}z_{k+3}, &k&=0,\\
  &\frac{1}{2}z_{k-1}-z_k+\frac{1}{2}z_{k+1},          &k&\geqslant 1,
 \end{aligned}\right.\quad
a_3=\left\{
 \begin{aligned}
  -&\frac{1}{6}z_k+\frac{1}{2}z_{k+1}-\frac{1}{2}z_{k+2}+\frac{1}{6}z_{k+3}, &k&=0,\\
  -&\frac{1}{6}z_{k-1}+\frac{1}{2}z_k-\frac{1}{2}z_{k+1}+\frac{1}{6}z_{k+2}, &k&\in[1:N_t-2],\\
  -&\frac{1}{6}z_{k-2}+\frac{1}{2}z_{k-1}-\frac{1}{2}z_k+\frac{1}{6}z_{k+1}, &k&=N_t-1.
 \end{aligned}\right.
\end{gathered}
\end{equation}
%

\section{Общий алгоритм}
\begin{enumerate}
\item Скорости и координаты преобразуются в прямоугольную сетку. Создаётся файл для интерполяции, содержащий информацию о сетке (число шагов $N_x$, $N_y$, $N_t$ и величина шага $\Delta x$, $\Delta y$, $\Delta t$ по каждой координате и времени) и скорости. Можно представить себе сетку как стопку плоскостей, расчерченых на прямоугольники. Каждая плоскость~--- это карта скоростей в фиксированый момент времени.
\item Для расчёта мы используем либо линейную интерполяцию в пространстве и времени, либо бикубическую в пространстве и полиномиальную во времени. Каждая компонента скорости интерполируется независимо от другой.
 \begin{enumerate}
  \item Определяем ячейку $(i,j,k)$, в которой находится точка $(X,Y,T)$ и нормированые координаты внутри этой ячейки
  \begin{equation}
   \begin{aligned}
   i&=[X/\Delta x], &j&=[Y/\Delta y], &k&=[T/\Delta t],\\ x&=\frac{X-i\Delta x}{\Delta x}, &y&=\frac{Y-j\Delta y}{\Delta y}, &t&=\frac{T-k\Delta t}{\Delta t},
  \end{aligned}
  \end{equation}
  где $[r]$~--- целая часть числа $r$.
  \item По формуле \eqref{linint} находим значения $u$ и $v$ на верхней и нижней гранях ячейки в точке $(x,y)$.
  \item По формуле \eqref{lint} находим значения $u$ и $v$ в точке $(x,y,t)$.
 \end{enumerate}
Для второго метода алгоритм чуть сложнее.
 \begin{enumerate}
  \item Ячейка и нормированые координаты находятся так же, как и в первом случае.
  \item Для четырёх временных плоскостей, ближайших к $T$, находятся по формуле \eqref{bicub_coeff} коэффициенты бикубического полинома \eqref{bicubic}\footnote{На самом деле они рассчитываются только один раз для каждого квадратика на плоскости, а потом хранятся в памяти. Пересчёт их при каждом определении скоростей~--- слишком дорогое удовольствие.}.
  \item Рассчитываются по формуле \eqref{bicubic} значения скоростей в точке $(x,y)$ на четырёх временных плоскостях.
  \item Скорости окончательно рассчитываются по формуле \eqref{cubt}.
 \end{enumerate}
\item Полученные скорости подставляются в уравнение \eqref{maineq} и интегрируются.
\item Выходные данные либо используются как есть (для анализа), либо трансформируются обратно в географические координаты (для статей).
\end{enumerate}
\end{document}