Browse Source

Added text about velocity fields interpolation

master
Michael Uleysky 5 years ago
parent
commit
7a7e9d30ad
  1. 217
      Interpolation/inter.tex
  2. 139
      Interpolation/net.eps

217
Interpolation/inter.tex

@ -0,0 +1,217 @@
\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}

139
Interpolation/net.eps

@ -0,0 +1,139 @@
%!PS-Adobe-2.0 EPSF-2.0
%%Title: net.eps
%%BoundingBox: 85 85 210 215
%%EndComments
gsave
100 100 translate
10 dict begin
/l 100 def
/r 3 def
/mlw 0.5 def
/r0 2 def
/rs 1.5 def
/x 0.45 def
/y 0.7 def
/dlw 0.5 def
/dash [ 2 2 ] def
/fs 15 def
/sfs fs 0.5 mul def
/dn fs -0.1 mul def
/shb -4 def
/shu 6 def
/strparams {20 dict begin /inp exch def gsave newpath 0 0 moveto
/stringtype {inp false charpath} def
/arraytype
{
inp {dup dup 0 get findfont exch 1 get scalefont setfont
dup 2 get /x exch def dup 3 get /y exch def
dup 5 get {x y rmoveto 4 get false charpath}{pop} ifelse
} forall
} def
inp type exec pathbbox
/ury exch def /urx exch def /lly exch def /llx exch def
% Return values
urx llx sub ury lly sub llx neg lly neg
grestore end} def
/GetParam { 3 dict begin
/default exch def /key exch def /dict exch def
dict key known {dict key get}{default} ifelse
end
} def
/Show {20 dict begin /Params exch def /String exch def
/HAlign Params /HAlign -1 GetParam def
/VAlign Params /VAlign -1 GetParam def
/UseHShift Params /UseHShift false GetParam def
/UseVShift Params /UseVShift false GetParam def
String strparams /VShift exch def /HShift exch def /Height exch def /Width exch def
% Correct HAlign and VAlign
HAlign 0 lt {/HAlign -1 def}{} ifelse HAlign 0 gt {/HAlign 1 def}{} ifelse
VAlign 0 lt {/VAlign -1 def}{} ifelse VAlign 0 gt {/VAlign 1 def}{} ifelse
% Move initial point
UseHShift {HShift 0 rmoveto}{} ifelse UseVShift {0 VShift rmoveto}{} ifelse
HAlign 0 eq {Width -0.5 mul 0 rmoveto}{} ifelse HAlign 1 eq {Width neg 0 rmoveto}{} ifelse
VAlign 0 eq {0 Height -0.5 mul rmoveto}{} ifelse VAlign 1 eq {0 Height neg rmoveto}{} ifelse
/stringtype {show} def
/arraytype {{
dup dup 0 get findfont exch 1 get scalefont setfont currentpoint 3 -1 roll dup dup 2 get exch 3 get rmoveto dup 4 get show
dup length 5 eq {pop pop pop}{
dup length 6 eq {5 get {pop pop}{moveto} ifelse}{
dup 5 get exch 6 get currentpoint 4 -2 roll % x0 y0 x y ux uy
{1}{3}ifelse index exch % x0 y0 x y ye ux
{2}{4}ifelse index exch %x0 y0 x y xe ye
moveto pop pop pop pop
}ifelse}ifelse
}forall} def
String dup type exec
end} def
0 setgray
newpath
0 0 r 0 360 arc
fill
newpath
0 l r 0 360 arc
fill
newpath
l l r 0 360 arc
fill
newpath
l 0 r 0 360 arc
fill
newpath
x l mul y l mul r0 0 360 arc
fill
newpath
x l mul 0 rs 0 360 arc
fill
newpath
0 y l mul rs 0 360 arc
fill
mlw setlinewidth
newpath
0 0 moveto
0 l lineto
l l lineto
l 0 lineto
0 0 lineto
stroke
dlw setlinewidth
dash 0 setdash
newpath
x l mul 0 moveto
x l mul y l mul lineto
0 y l mul lineto
stroke
0 0 moveto [
[ /Times-Italic fs 0 shb (z) true]
[ /Times-Italic sfs 0 dn (00) true]
] << /HAlign 0 /VAlign 1>> Show
l 0 moveto [
[ /Times-Italic fs 0 shb (z) true]
[ /Times-Italic sfs 0 dn (10) true]
] << /HAlign 0 /VAlign 1>> Show
0 l moveto [
[ /Times-Italic fs 0 shu (z) true]
[ /Times-Italic sfs 0 dn (01) true]
] << /HAlign 0 /VAlign -1>> Show
l l moveto [
[ /Times-Italic fs 0 shu (z) true]
[ /Times-Italic sfs 0 dn (11) true]
] << /HAlign 0 /VAlign -1>> Show
l x mul 0 moveto [
[ /Times-Italic fs 0 shb (x) true]
] << /HAlign 0 /VAlign 1>> Show
0 l y mul moveto [
[ /Times-Italic fs -3 0 (y) true]
] << /HAlign 1 /VAlign 0>> Show
l x mul l y mul moveto [
[ /Times-Italic fs 3 3 (z) true]
] << /HAlign -1 /VAlign -1>> Show
end
grestore
%%EOF
Loading…
Cancel
Save