\documentclass[11pt,reqno]{amsart}
\usepackage[cp1251]{inputenc}
\usepackage[english,russian]{babel}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{graphicx}
\usepackage{xcolor}
\usepackage[colorlinks]{hyperref}


\usepackage{titlesec}
\titlespacing{\section}{0pt}{2.5ex}{1.5ex}
\titlespacing{\subsection}{0pt}{1.5ex}{1ex}
\titlespacing{\subsubsection}{0pt}{1.5ex}{1ex}
\titleformat{\section}{\large\bfseries\centering}{\thesection}{1em}{}
\titleformat{\subsection}[runin]{\bfseries}{\thesubsection.}{0.5em}{}[.\mbox{\ }]
\titleformat{\subsubsection}[runin]{\bfseries}{\thesubsubsection.}{0.4em}{}[.\mbox{\ }]



%%%%%---------- ORCiD
\newcommand\orcidicon[1]{\href{https://orcid.org/#1}{\includegraphics[scale=0.02]{orcid.pdf}}}

\def\udcs{519.633} %Здесь автор определяет УДК своей работы
\def\mscs{65M25} %Здесь автор определяет классификаторы AMS своей работы
\setcounter{page}{144}

\newtheorem{lemma}{Лемма}
\newtheorem{theorem}{Теорема}
\newtheorem{definition}{Определение}
\newtheorem{corollary}{Следствие}
\newtheorem{proposition}{Предложение}


\definecolor{myblue}{rgb}{0.03, 0.27, 0.79}

\hypersetup{linkcolor=blue,filecolor=black,urlcolor=blue, citecolor=blue}

\def\logo{\textcolor{blue}{{\bf\Huge S\raisebox{0.2ex}{\hspace{0.55ex}\raisebox{0.05ex}e\hspace{-1.65ex}$\bigcirc$}MR}}}

\def\semrtop
     {
  \vbox{
     \noindent\logo \hspace{-3mm}  \parbox{12cm}{     \begin{center}
     {\LARGE \textsf{СИБИРСКИЕ ЭЛЕКТРОННЫЕ}} \\[2mm]
     {\LARGE \textsf{МАТЕМАТИЧЕСКИЕ ИЗВЕСТИЯ}} \\[2mm]
     {\Large \textsf{Siberian Electronic Mathematical Reports}} \\[1mm]
     {\Large\tt{\textcolor{blue}{http://semr.math.nsc.ru}}}\\[1mm]
     \raisebox{1ex}{ISSN 1813-3304 }
     \vspace*{-3mm}
%     {\small Институт математики им. С.Л.Соболева СО РАН}\\[-1mm]
%     {\small Sobolev Institute of Mathematics SB RAS}
     \end{center}
}
     \vspace{-1mm}
     \noindent
      \begin{tabular}{c}
     \hphantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa!} \\
     \hline\hline
     \end{tabular}

     \vspace{3mm}
     {\flushleft\it Том 20, № 2, стр. 144--145 (2023) \hspace{\fill}{\rm\small УДК \udcs \hspace{-4mm}}}
     \newline
     	{\rm\small \textcolor{blue}{https://doi.org/10.33048/semi.2023.20.???}}\hphantom{aaaaaaaaaaaaaaaaaaaaaaaa\;}{\rm\small MSC\ \ \mscs }
  }%\vbox
}



\begin{document}
\renewcommand{\refname}{References}


\thispagestyle{empty}

\title[\sc Разные шаги по времени для уравнения неразрывности]{\large Вычислительная технология использования разных шагов по времени в разных частях вычислительной области при решении двумерного уравнения неразрывности консервативным полулугранжевым методом}
\author[А.В. Вяткин, Е.В. Кучунова]{{\bf \large \textsc{А.В. Вяткин}}\protect\orcidicon{0000-0002-7609-5398}, {\bf \large \textsc{Е.В. Кучунова}}\protect\orcidicon{0000-0009-9800-4184}\\ \\
%\textit{Представлено} \textsc{П.П. Петровым}
%\textit{(заполняется редактором)}
}%
\address{Vyatkin Aleksandr Vladimirovich  % Контактные данные всех авторов и место работы указываются только на английском языке
\newline\hphantom{iii} Institute of Computational Modelling SB RAS,
\newline\hphantom{iii} str. Akademgorodok, 50/44,
\newline\hphantom{iii} 660036, Krasnoyarsk, Russia}%
\email{\textcolor{blue}{vyatkin@icm.krasn.ru}}%

\address{Kuchunova Elena Vladimirovna  % Контактные данные всех авторов и место работы указываются только на английском языке
\newline\hphantom{iii} Siberian Federal University,
\newline\hphantom{iii} pr. Svobodny, 79,
\newline\hphantom{iii} 660041, Krasnoyarsk, Russia}%
\email{\textcolor{blue}{HKuchunova@sfu-kras.ru}}%

\thanks{\sc Vyatkin A.V., Kuchunova E.V.
Algorithm of different time step usage in different parts of computational domain based on semi-Lagrangian method for continuity equation}
\thanks{\copyright \ 2024 Вяткин А.В, Кучунова Е.В.}
\thanks{\rm Работа поддержана Красноярским математическим центром, финансируемым Минобрнауки РФ в рамках мероприятий по созданию и развитию региональных НОМЦ (Соглашение 075-02-2025-1606)}
\thanks{\it Поступила 1 января 2023 г., опубликована 31 декабря 2023 г.}%


\semrtop \vspace{1cm}
\maketitle {\small
\begin{quote}
\noindent{\bf Abstract:} In the paper semi-Lagrangian algorithm for two-dimensional continuity equation is considered. Computational domain is decomposed into two parts. Various time steps in different parts of computational domain are used. Theoretical investigation and numerical experiments confirm first order of convergence and conservation low validity.\medskip

\noindent{\bf Keywords:} continuity equation, semi-Lagrangian method, conservation low, irregular grid, Eulerian-Lagrangian method.
 \end{quote}
}

\bigskip

\section{Введение}

Свое развитие полулагранжевые методы, ранее называемые методом характеристик или методом траекторий, начали в середине XX века при моделировании движении газа \cite{1}-\cite{4} в задачах моделирования обтекания тела, а также в задаче прогнозирования погоды \cite{5}. Главным свойством этих методов является то, что они основаны на моделировании перемещения субстанции посредством вычисления траекторий движения. Со временем эти методы начали применяться для решения системы уравнений Навье-Стокса \cite{6}-\cite{9}, а также более широкого класса задач, в которых, в качестве одного из процессов, присутствует адвекция \cite{10}-\cite{14}. 

У этих методов можно выделить три главных преимущества. Первое преимущество состоит в том, что для этих методов условие устойчивости Куранта-Фридрихса-Леви выполняется автоматически (с точностью до погрешности аппроксимации траекторий движения субстанции) \cite{10}, \cite{15}, \cite{16}, что позволяет проводить расчеты с большими шагами по времени. Второе преимущество состоит в возможности использовать разные пространственные сетки на разных слоях по времени. Это преимущество привело к активному использованию адаптивных сеток \cite{17}-\cite{19}. Третье преимущество состоит в выполнении закона сохранения для численного решения в дискретном аналоге нормы пространства $L_{1} $ \cite{11}, \cite{15}, \cite{17}, \cite{20}. Подробное описание преимуществ полулагранжевого подхода можно найти в статьях \cite{10}, \cite{19}, \cite{21}-\cite{23} вместе с обширным обзором литературы.

В данной работе показано наличие еще одного преимущества, которое состоит в возможности использовать разные шаги по времени в разных частях вычислительной области. При этом, при моделировании осуществляется сквозное течение субстанции из одной части вычислительной области в другую и обратно. Теоретически и на основе вычислительного эксперимента подтверждена сходимость численного решения к точному решению с первым порядком точности, а также выполнение закона сохранения без использования каких-либо поправочных (весовых) коэффициентов, искажающих численное решение. В работе вычислительная область разбита на две части, в каждой из которых используются разные шаги по времени. Основная сложность при описании метода состоит в построении численного решения в узлах сетки, лежащих на границе между двумя подобластями. В данной работе теоретически обоснована и экспериментально показана принципиальная возможность создания алгоритма с указанными свойствами. %При этом в работе не обсуждается экономия во времени проведения расчетов и увеличение погрешности в численном решении от использования данного алгоритма. 

\section{Постановка задачи}

Пусть $\Omega = \left[0,1\right]\times \left[0,1\right]$ -- единичный квадрат с границей $\Gamma $. На множестве $D =\left[0, T \right] \times \Omega $, $T > 0$, рассмотрим двумерное уравнение неразрывности с ненулевой правой частью
%
\begin{equation} \label{GrindEQ__1_} 
\frac{\partial \rho }{\partial t} +\frac{\partial (\rho u)}{\partial x} +\frac{\partial (\rho v)}{\partial y} =f(t,x,y). 
\end{equation} 
%
Здесь $\rho \left(t, x, y \right)$ -- искомая функция, а $u\left(t,x,y\right)$, $v\left(t,x,y\right)$ и $f\left(t,x,y\right)$ -- известные непрерывно-дифференцируемые функции из класса $C^1(D)$. Для простоты будем полагать, что на границе $\Gamma $ выполняется условие 
%
\begin{equation} \label{GrindEQ__2_} 
u(t,x,y)=v(t,x,y)=0  \quad \forall {\rm \; }\left(t,x,y\right)\in \Gamma .  
\end{equation} 
%
Для искомой функции $\rho \left(t,x,y\right)$ зададим начальное условие 
%
\begin{equation} \label{GrindEQ__3_} 
\rho \left(0,x,y\right)=\rho _{{\rm init}} (x,y)  \quad \forall \left(x,y\right)\in \Omega, \end{equation} 
%
где $\rho_{{\rm init}} (x,y)$ -- известная функция. На множестве $\Omega $ введем сетку $$G_{h}^{\Omega } =\left\{(x_{i} ,y_{j}) : x_{i} =ih, y_{j} =ih; i,j=0,...,N\right\}$$ с шагом $h=1/N$, $N=10\cdot 2^{n} $, $n\in {\mathbb N}$. Пусть числа $N_{1} $ и $N_{2} $ такие, что $N_{1} \cdot h=0,2$ и $N_{2} \cdot h=0,8$. Разделим множество узлов $G_{h}^{\Omega } $ на два подмножества:
%
\begin{equation*}
G_{2\tau }^{\Omega } = \left\{(x_{l} ,y_{n} ) : l,n = N_{1} +1, \dots, N_{2} -1\right\} \; \text{и} \; G_{\tau }^{\Omega } = G_{h}^{\Omega } \backslash G_{2\tau }^{\Omega }.    
\end{equation*}
%
 Множество узлов сетки $G_{\tau }^{\Omega } $, которые принадлежат внутренней границе множества $G_{\tau }^{\Omega } $, обозначим через $\Gamma _{\tau }^{\Omega } $. Множество узлов сетки $G_{2\tau }^{\Omega } $, которые принадлежат внешней границе множества $G_{2\tau }^{\Omega } $, обозначим через $\Gamma _{2\tau }^{\Omega } $. Определим множество узлов, принадлежащих границам смежных сеток $G_{\tau }^{\Omega } $ и $G_{2\tau }^{\Omega } $ в виде $\Gamma _{h}^{\Omega } =\Gamma _{\tau }^{\Omega } \cup \Gamma _{2\tau }^{\Omega } $, рис.~\ref{Pic0} (а). Обозначим $\tau =T/K$, где $K$~--~четное число. Положим $G_{h}^{D} =G_{\tau }^{D} \cup G_{2\tau }^{D} $, $\Gamma _{h}^{D} =\Gamma _{\tau }^{D} \cup \Gamma _{2\tau }^{D} $ где 
%
\begin{equation*}
    \begin{gathered}
        G_{\tau }^{D} =\left\{(t_{k} ,x_{i} ,y_{j} ):t_{k} =k\tau ,(x_{i} ,y_{j} )\in G_{\tau }^{\Omega }; \; k=0,...,K\right\}, \\   
        G_{2\tau }^{D} =\left\{(t_{2m} ,x_{l} ,y_{n} ):t_{2m} =2m\tau ,(x_{l} ,y_{n} )\in G_{2\tau }^{\Omega }; \; m=0,...,K/2\right\}, \\
        \Gamma _{\tau }^{D} =\left\{(t_{k} ,x_{i} ,y_{j} ):t_{k} =k\tau ,(x_{i} ,y_{j} )\in \Gamma _{\tau }^{\Omega }; \; k=0,...,K\right\}, \\
        \Gamma _{2\tau }^{D} =\left\{(t_{2m} ,x_{l} ,y_{n} ):t_{2m} =2m\tau ,(x_{l} ,y_{n} )\in \Gamma _{2\tau }^{\Omega }; \; m=0,...,K/2\right\}.        
    \end{gathered}    
\end{equation*}
%
Для каждого узла $\left(t_{k} ,x_{p} ,y_{q} \right)$ сетки $G_{h}^{D} $ рассмотрим его окрестность 
%
\begin{equation*}
    \Omega _{p,q}^{k} =\left\{t=t_{k} \right\}\times \left(\left[x_{p} -h/2,x_{p} +h/2\right)\times \left[y_{q} -h/2,y_{q} +h/2\right)\cap \Omega\right). 
\end{equation*}
%
Решение задачи \eqref{GrindEQ__1_}--\eqref{GrindEQ__3_} будем искать в окрестностях $\Omega _{p,q}^{k} $ всех узлов сетки $G_{h}^{D} $ в виде кусочно-постоянной функции $\rho ^{h} \left(t,x,y\right)$ равной константе в каждой окрестности $\Omega _{p,q}^{k} $ узла $\left(t_{k} ,x_{p} ,y_{q} \right)$
%
\begin{equation*}
    \rho ^{h} \left(t_{k} ,x,y\right)=\rho ^{h} \left(t_{k} ,x_{p} ,y_{q} \right)   \forall {\rm \; }\left(t_{k} ,x,y\right)\in \Omega _{p,q}^{k}.  
\end{equation*}
%
Введем обозначения $g\left(t_{k} ,x_{p} ,y_{q} \right)=g_{p,q}^{k} $ и $\rho ^{h} \left(t_{k} ,x_{p} ,y_{q} \right)=\rho _{p,q}^{h,k} $. Рассмотрим узел $\left(x_{p} ,y_{q} \right)\in G_{h}^{\Omega } $. Пусть
%
\begin{equation}
\nonumber 
    \Omega _{p,q} =\left[x_{p} -h/2,x_{p} +h/2\right)\times \left[y_{q} -h/2,y_{q} +h/2\right)\cap D.
\end{equation}
%
Введем обозначения
%
\begin{equation*}
    \Omega _{\tau }^{} = \bigcup _{\left(x_{p} ,y_{q} \right)\in G_{\tau }^{\Omega } }\Omega _{p,q},
    \qquad
    \Omega _{2\tau } =\bigcup _{\left(x_{p} ,y_{q} \right)\in G_{2\tau }^{\Omega } }\Omega _{p,q},
    \qquad 
    \Omega _{\Gamma } =\bigcup _{\left(x_{p} ,y_{q} \right)\in \Gamma _{h}^{\Omega } }\Omega _{p,q},
\end{equation*}
%
\begin{equation*}
    D_{\tau } =\left[0,T\right]\times \Omega _{\tau }, 
    \quad 
    D_{2\tau } =\left[0,T\right]\times \Omega _{2\tau }, 
    \quad
    D_{\Gamma } =\left[0,T\right]\times \Omega _{\Gamma }.
\end{equation*}
%
Понятно, что $\Omega _{\Gamma } \cap \Omega _{\tau } \ne \emptyset $, $\Omega _{\Gamma } \cap \Omega _{2\tau } \ne \emptyset $. Внутреннюю границу множества $D_{\tau } $, находящуюся внутри множества $D_{\Gamma } $, обозначим $S$, рис.~\ref{Pic0} (б). 
%
\begin{figure}
    \centering
    а) \includegraphics[scale=0.9]{Pic0a.png} \\
    б) \includegraphics[scale=0.9]{Pic0b.png}
    \caption{а) Множества узлов смежных сеток $G_{\tau }^{\Omega }$, $G_{2\tau }^{\Omega } $ и $\Gamma _{h}^{\Omega }$ б) Множества $\Omega_{\tau}$, $\Omega_{2\tau}$ и граница $S$}
    \label{Pic0}
\end{figure}
%
Для упрощения теоретических выкладок положим, что 
%
\begin{equation} \label{GrindEQ__4_} 
\tau = c_{\tau ,h} h,  
\end{equation} 
%
где $c_{\tau ,h}$ -- некоторая константа такая, что
%
\begin{equation} \label{GrindEQ__5_} 
c_{\tau ,h} \leq\frac{1}{4\mathop{\max }\limits_{\left(t,x,y\right)\in D} \left\{\left|u\left(t,x,y\right)\right|,\left|v\left(t,x,y\right)\right|\right\}}  
\end{equation} 
%
В общем случае представленного метода ограничение \eqref{GrindEQ__5_} может быть ослаблено так, что может использоваться константа $c_{\tau ,h}$ большего значения. В данной работе ограничение \eqref{GrindEQ__5_} используется для упрощения программной реализации описываемого метода.

\section{Интегральное тождество}

Рассмотрим узел $\left(t_{k} ,x_{p} ,y_{q} \right)\in G_{h}^{D} $, который не принадлежит множеству $\Gamma _{h}^{D}$. Пусть точка $A^{k} =\left(t_{k} ,A^{k,x} ,A^{k,y} \right)$ принадлежит границе множества $\bar{\Omega }_{p,q}^{k} $, где $\bar{\Omega }_{p,q}^{k} $ -- замыкание области $\Omega _{p,q}^{k} $. Построим траекторию движения $C_{A^{k} } \left(t,\tilde{x}\left(t;A^{k} \right),\tilde{y}\left(t;A^{k} \right)\right)$ точки $A^{k} $ обратно по времени со слоя $t=t_{k} $ на слой $t=t_{k-s} $. Здесь $s=1$, если $\left(t_{k} ,x_{p} ,y_{q} \right)\in G_{\tau }^{D} $, и $s=2$, если $\left(t_{k} ,x_{p} ,y_{q} \right)\in G_{2\tau }^{D} $. Траектория $C_{A^{k} } \left(t,\tilde{x}\left(t;A^{k} \right),\tilde{y}\left(t;A^{k} \right)\right)$ описывается решением задачи Коши для следующей системы обыкновенных дифференциальных уравнений
%
\begin{equation} \label{GrindEQ__6_} 
\frac{d\tilde{x}}{dt} \left(t; A^{k} \right)=u\left(t,\tilde{x},\tilde{y}\right),\, \, \, \frac{d\tilde{y}}{dt} \left(t; A^{k} \right)=v\left(t,\tilde{x},\tilde{y}\right),\, \, \, \, \left[t\in t_{k-s} ,t_{k} \right], 
\end{equation} 
%
с начальными данными
%
\begin{equation} \label{GrindEQ__7_} 
\tilde{x}\left(t_{k}; A^{k} \right)=A^{k,x}, \quad  \tilde{y}\left(t_{k}; A^{k} \right)=A^{k,y} .  
\end{equation} 
В силу ограничения \eqref{GrindEQ__5_} траектория 
$C_{A^{k} } \left(t,\tilde{x}\left(t;A^{k} \right),\tilde{y}\left(t;A^{k} \right)\right)$ не пересекает поверхность $S$. Множество $C^{{\rm set}} \left(\bar{\Omega }_{p,q}^{k} \right)$ траекторий, выпущенных из всех точек границы множества $\bar{\Omega }_{p,q}^{k} $, формируют в плоскости $t=t^{*} $ замкнутую кривую, которую рассмотрим в виде границы четырехугольника $Q_{p,q} \left(t^{*} \right)$. Обозначим $Q_{p,q}^{k-s} =Q_{p,q} \left(t_{k-s} \right)$, рис.~\ref{Pic1}. Пусть $V_{p,q}^{k} $ -- множество, ограниченное сверху четырехугольником $\left\{t=t_{k} \right\}\times \bar{\Omega }_{p,q}^{k} $, снизу -- четырехугольником $\left\{t=t_{k-s} \right\}\times Q_{p,q}^{k-s} $, по бокам -- множеством траекторий $C^{{\rm set}} \left(\bar{\Omega }_{p,q}^{k} \right)$. 
%
\begin{figure}
    \centering
    \includegraphics{Pic1-new.png}
    \caption{Формирование четырехугольников $Q_{N_{1} -1,q}^{k-1} $ и $Q_{N_{1} +2,q}^{k-2}$}
    \label{Pic1}
\end{figure}
\noindent 
%
Учитывая теорему Гаусса-Остроградского, несложно доказать следующую теорему.

\begin{theorem}
Решение задачи \eqref{GrindEQ__1_}--\eqref{GrindEQ__3_} удовлетворяет уравнению
%
\begin{equation} \label{GrindEQ__8_} 
\int\limits_{\Omega _{p,q}^{k} }\rho  \left(t_{k} ,x,y\right)dxdy=\int\limits_{Q_{p,q}^{k-s} }\rho  \left(t_{k-s} ,x,y\right)dxdy
+
\int\limits_{V_{p,q}^{k} }f\left(t,x,y\right)dtdxdy .  
\end{equation} 
\end{theorem}

\section{Построение численного решения в узлах, не принадлежащих границам сеток на стыке сеток}

Чтобы построить алгоритм, аппроксимируем каждое слагаемое из равенства \eqref{GrindEQ__8_}. Обозначим ${\rm meas}\left(\Omega _{p,q}^{k} \right)$ -- площадь окрестности $\Omega _{p,q}^{k} $. Раскладывая функцию $\rho \left(t_{k} ,x,y\right)$ по формуле Тейлора в окрестности точки $\left(t_{k} ,x_{p} ,y_{q} \right)$ по переменным $x$ и $y$, несложно показать, что 

\begin{equation} \label{GrindEQ__9_} 
\begin{split} {\forall \left(x_{p} ,y_{q} \right)\notin \Gamma {\rm \; \; \; }\int\limits_{\Omega _{p,q}^{k} }\rho  \left(t_{k} ,x,y\right)dxdy=\rho^{k}_{p,q} {\rm meas}\left(\Omega _{p,q}^{k} \right)+O\left(h^{4} \right),} 
\\
{\forall \left(x_{p} ,y_{q} \right)\in \Gamma {\rm \; \; \; }\int\limits _{\Omega _{p,q}^{k} }\rho  \left(t_{k} ,x,y\right)dxdy=\rho^{k}_{p,q} {\rm meas}\left(\Omega _{p,q}^{k} \right)+O\left(h^{3} \right).} 
\end{split} 
\end{equation} 
Для вычисления интеграла по множеству $Q_{p,q}^{k-s} $, мы аппроксимируем $Q_{p,q}^{k-s} $ четырехугольником $P_{p,q}^{k-s} $ с прямыми сторонами. Для этого рассмотрим четыре вершины $A_{d}^{k} =\left(t_{k} ,A_{d}^{k,x} ,A_{d}^{k,y} \right)$, $d=1,...,4$, четырехугольника $\bar{\Omega }_{p,q}^{k} $. Траектории $C_{A_{d}^{k} } \left(t,\tilde{x}\left(t;A_{d}^{k} \right),\tilde{y}\left(t;A_{d}^{k} \right)\right)$ движения точек $A_{d}^{k} $ обратно по переменной времени со слоя $t=t_{k} $ на слой $t=t_{k-s} $ образуют в плоскости $t=t_{k-s} $ точки $B_{d}^{k-s}=\left(t_{k-s} ,B_{d}^{k-s,x} ,B_{d}^{k-s,y} \right)$, $d=1,...,4$. Найдем приближенные координаты $\left(t_{k-s} ,B_{d}^{h,k-s,x} ,B_{d}^{h,k-s,y} \right)$ точек $B_{d}^{k-s} \left(t_{k-s} \right)$ и обозначим $B_{d}^{h,k-s} =\left(t_{k-s} ,B_{d}^{h,k-s,x} ,B_{d}^{h,k-s,y} \right)$, $d=1,...,4$. Для этого решим задачу Коши \eqref{GrindEQ__6_}-\eqref{GrindEQ__7_} методом Эйлера
% 
\begin{equation} \label{GrindEQ__10_}
\begin{gathered}
    B_{d}^{h,k-s,x} = A_{d}^{k,x} -s\tau u\left(t_{k} ,A_{d}^{k,x} ,A_{d}^{k,y} \right), \\
    B_{d}^{h,k-s,y} =A_{d}^{k,y}-s\tau v\left(t_{k} ,A_{d}^{k,x} ,A_{d}^{k,y} \right),    d=1,...,4.
\end{gathered}
\end{equation} 
%
Здесь $s=1$, если $\left(t_{k} ,x_{p} ,y_{q} \right)\in G_{\tau }^{D} $, и $s=2$, если $\left(t_{k} ,x_{p} ,y_{q} \right)\in G_{2\tau }^{D} $. Таким образом, каждую траекторию $C_{A_{d}^{k} } \left(t,\tilde{x}\left(t;A_{d}^{k} \right),\tilde{y}\left(t;A_{d}^{k} \right)\right)$ мы аппроксимируем прямым отрезком. Несложно показать второй порядок погрешность такой аппроксимации 
\begin{equation} \label{GrindEQ__11_} 
\left|B_{d}^{h,k-s,x} -B_{d}^{k-s,x} \right|=O\left(\tau ^{2} \right),   \left|B_{d}^{h,k-s,y} -B_{d}^{k-s,y} \right|=O\left(\tau ^{2} \right).  
\end{equation} 
%
Соединим точки $B_{d}^{h,k-s}$, $d=1,...,4$, прямыми отрезками и получим четырехугольник $P_{p,q}^{k-s} $, аппроксимирующий $Q_{p,q}^{k-s} $, рис.~\ref{Pic2}. 
%
\begin{figure}
    \centering
    \includegraphics[scale=0.9]{Pic2.png}
    \caption{Аппроксимация четырехугольника $Q_{p,q}^{k-s}$}
    \label{Pic2}
\end{figure}

Рассмотрим два множества $Q_{p,q}^{k-s} \backslash P_{p,q}^{k-s} $, $P_{p,q}^{k-s} \backslash Q_{p,q}^{k-s} $ и используем равенство 
%
\begin{equation}
    \label{GrindEQ__12_} 
    \begin{gathered}
     \int\limits_{Q_{p,q}^{k-s} }\rho  \left(t_{k-s} ,x,y\right)dxdy =  \int\limits_{P_{p,q}^{k-s} }\rho  \left(t_{k-s} ,x,y\right)dxdy + 
     \qquad \qquad \qquad \\
     \quad \qquad  + \int\limits_{Q_{p,q}^{k-s} \backslash P_{p,q}^{k-s} }\rho  \left(t_{k-s} ,x,y\right)dxdy - \int\limits _{P_{p,q}^{k-s} \backslash Q_{p,q}^{k-s} }\rho  \left(t_{k-s} ,x,y\right)dxdy.      
     \end{gathered}
\end{equation} 
%
Обозначим 
%
\begin{equation} 
\label{GrindEQ__13_} 
\delta Q_{p,q}^{h,k-s} =\int\limits_{Q_{p,q}^{k-s} \backslash P_{p,q}^{k-s} }\rho  \left(t_{k-s} ,x,y\right)dxdy-\int\limits_{P_{p,q}^{k-s} \backslash Q_{p,q}^{k-s} }\rho  \left(t_{k-s} ,x,y\right)dxdy. 
\end{equation} 
%
Известно \cite{24}, при выполнении равенства \eqref{GrindEQ__4_} оценку величины $\left|\delta Q_{p,q}^{h,k-s} \right|$ можно выписать в виде
%
\begin{equation} \label{GrindEQ__14_} 
\left|\delta Q_{p,q}^{h,k-s} \right|=c_{Q} \tau ^{2} h^{2} +o\left(\tau ^{2} h^{2} \right).  
\end{equation} 
%
Здесь $c_{Q} $ -- некоторая константа независящая от $\tau $ и $h$, а обозначение $o\left(\tau ^{2} h^{2} \right)$ означает величину асимптотически меньшую, чем $\tau^{2} h^{2}$, учитывая равенство \eqref{GrindEQ__4_}. 

Далее, рассмотрим интеграл от функции $f\left(t,x,y\right)$ в \eqref{GrindEQ__8_}. Пусть 
\begin{equation*}
    F\left(t\right)=\int\limits_{Q_{p,q} \left(t\right)}f\left(t,x,y\right)dxdy.  
\end{equation*}
\noindent 
%
Используя разложения функции $F\left(t\right)$ по формуле Тейлора с остаточным членом в форме Лагранжа, несложно показать, что 
%
\begin{equation} \label{GrindEQ__15_} 
\int\limits_{V_{p,q}^{k} }f\left(t,x\right)dtdxdy =\int\limits_{t_{k} -s\tau }^{t_{k} }F\left(t\right)dt =s\tau F\left(t_{k} \right)+\frac{\partial F}{\partial t} \left(t^{*} \right)s^{2} \tau ^{2} ,   
\end{equation} 
%
где $t^{*} $ -- некоторая величина, такая что $t^{*} \in \left[t_{k} -s\tau ,t_{k} \right]$. Определим множества $\Omega _{p,q} \backslash Q_{p,q} \left(t^{*} \right)$ и $Q_{p,q} \left(t^{*} \right)\backslash \Omega _{p,q} $ такими, что 
%
\begin{equation} \label{GrindEQ__16_} 
\begin{gathered}
    \int\limits_{\Omega _{p,q} }f\left(t,x,y\right)dxdy = \int \limits_{Q_{p,q} \left(t^{*} \right)}f\left(t,x,y\right)dxdy + 
    \qquad \qquad \qquad \qquad
    \\
    \qquad \qquad + \int \limits_{\Omega _{p,q} \backslash Q_{p,q} \left(t^{*} \right)}f\left(t,x,y\right)dxdy -  \int \limits_{Q_{p,q} \left(t^{*} \right)\backslash \Omega _{p,q} }f\left(t,x,y\right)dxdy .  
\end{gathered}
\end{equation} 
%
По аналогии с \eqref{GrindEQ__14_}, при выполнении равенства \eqref{GrindEQ__4_} справедливо
%
\begin{equation} \label{GrindEQ__17_} 
\begin{gathered}
    \left|\int \limits_{\Omega _{p,q} \backslash Q_{p,q} \left(t^{*} \right)}f\left(t_{k} ,x,y\right)dxdy -\int \limits_{Q_{p,q} \left(t^{*} \right)\backslash \Omega _{p,q} }f\left(t_{k} ,x,y\right)dxdy \right| = 
    \qquad
    \\ 
    \qquad
    \qquad \qquad \qquad \qquad
    = c_{Q}^{f} \left(t_{k} -t^{*} \right)^{2} {\rm meas}\left(\Omega _{p,q} \right)+o\left(\left(t_{k} -t^{*} \right)^{2} h^{2} \right). 
\end{gathered}
\end{equation} 
%
где $c_{Q}^{f} $ -- некоторая константа, независящая от величин $\left(t_{k} -t^{*} \right)$ и $h$. Поскольку функция $f(t,x,y)$ достаточно гладкая, то существует производная $\frac{\partial F}{\partial t} \left(t\right)$ в точке $t = t^{*}$. Оценим ее величину с учетом выражений \eqref{GrindEQ__16_}, \eqref{GrindEQ__17_} следующим образом:
%
\begin{equation*}
\begin{gathered} 
    \frac{\partial F}{\partial t} \left(t^{*} \right) 
    =
    \mathop{\lim }\limits_{t_{k} -t^{*} \to 0} \frac{F\left(t_{k} \right)-F\left(t^{*} \right)}{t_{k} -t^{*} } = \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad 
    \\ 
    = \mathop{\lim }\limits_{t_{k}-t^{*} \to 0} \frac{1}{t_{k} -t^{*} } \left(\int \limits_{\Omega _{p,q} }f\left(t_{k} ,x,y\right)dxdy -\int \limits_{Q_{p,q} \left(t^{*} \right)}f\left(t^{*} ,x,y\right)dxdy \right) \le 
    \\
    \qquad \qquad \qquad \qquad \qquad \qquad  \le
%     \mathop{\lim }\limits_{t_{k}-t^{*} \to 0} \frac{1}{t_{k} -t^{*} } \left(
%     c^f_Q  \left( t_k -t^* \right)^2 {\rm meas}\left(\Omega_{p,q} \right)   
%     +o\left(\left(\left( t_k -t^* \right)^2 h\right)^{2} \right)
%          \right)=
%          \\
     c^f_Q \left( t_k -t^* \right) {\rm meas}\left(\Omega_{p,q}^{\max} \right)   
     +o\left(\left( t_k -t^* \right) h^{2} \right).
\end{gathered}
\end{equation*}
%
Тогда,
%
\begin{equation} \label{GrindEQ__18_} 
\begin{gathered}
\int\limits _{V_{p,q}^{k} }f\left(t,x\right)dtdxdy 
=
\int \limits_{t_{k} -s\tau }^{t_{k} }F\left(t\right)dt 
=
s\tau F\left(t_{k} \right)+\delta F\left(t_{k} \right)
\end{gathered}
\end{equation} 
%
где $\left| \delta F\left(t_{k} \right)\right| \le  c^f_Q s^2 \tau^3 h^2 +o\left( \tau^3 h^{2} \right)$. По аналогии с \eqref{GrindEQ__9_} очевидно, что
%
\begin{equation} \label{GrindEQ__19_} 
\begin{array}{l} {\forall \left(x_{p} ,y_{q} \right)\notin \Gamma \; s\tau \int \limits_{\Omega _{p,q}^{k} }f\left(t_{k} ,x,y\right)dx =s\tau f^k_{p,q}{\rm meas}\left(\Omega _{p,q}^{k} \right)+O\left(\tau h^{4} \right),} 
\\
{\forall \left(x_{p} ,y_{q} \right)\in \Gamma \; s\tau \int \limits_{\Omega _{p,q}^{k} }f\left(t_{k} ,x,y\right)dx =s\tau f^k_{p,q}{\rm meas}\left(\Omega _{p,q}^{k} \right)+O\left(\tau h^{3} \right).} \end{array} 
\end{equation} 
%
Обозначение $O\left(\tau h^{3} \right)$ означает величину асимптотически равную $\tau h^{3}$, учитывая равенство \eqref{GrindEQ__4_}.

В итоге, подставим \eqref{GrindEQ__9_}, \eqref{GrindEQ__12_}, \eqref{GrindEQ__18_}, \eqref{GrindEQ__19_} в \eqref{GrindEQ__8_} и выпишем формулу для определения численного решения $\rho _{p,q}^{h,k} $ 
%
\begin{equation} \label{GrindEQ__20_} 
\rho _{p,q}^{h,k} =\frac{1}{{\rm meas}\left(\Omega _{p,q}^{k} \right)} \int \limits_{P_{p,q}^{k-s} }\rho ^{h} \left(t_{k-s} ,x,y\right)dxdy +s\tau f_{p,q}^{k} .  
\end{equation} 
%
Напомним, что функция $\rho ^{h} \left(t_{m-s} ,x,y\right)$ является кусочно-постоянной. Будем считать, что интеграл в \eqref{GrindEQ__20_} вычисляется точно. Этого можно достичь путем разложения четырёхугольника $P_{p,q}^{k-s} $ на составные части, в каждой из которых функция $\rho ^{h} \left(t_{m-s} ,x,y\right)$ равна константе.

\section{Построение численного решения в узлах, принадлежащих границам сеток на стыке сеток}

Опишем алгоритм построения решения в узлах сетки $\left(t_{k} ,x_{p} ,y_{q} \right)\in \Gamma _{h}^{D} $. Чтобы построить численное решение в узле $\left(t_{k} ,x_{p} ,y_{q} \right)\in \Gamma _{h}^{D} $, рассмотрим его окрестность $\Omega _{p,q}^{k} $ и точки $B_{d}^{h,k-1} =\left(t_{k-1} ,B_{d}^{h,k-1,x} ,B_{d}^{h,k-1,y} \right)$, $d=1,...,4$, вычисленные по формуле \eqref{GrindEQ__10_} при $s=1$. 

Пусть $k$ является нечетным числом. Поскольку численное решение $\rho ^{h} \left(t_{k-1} ,x,y\right)$ на четном слое $t=t_{k-1} $ считаем уже вычисленным на всем множестве $\left\{t=t_{k-1} \right\}\times \Omega $, то для вычисления $\rho _{p,q}^{h,k} $ используем формулу \eqref{GrindEQ__20_} при $s=1$.

Пусть $k$ является четным числом. Определим координаты вершин четырехугольника $P_{p,q}^{k-1} $ по формуле \eqref{GrindEQ__10_} при $s=1$. В зависимости от значений функций скорости $u\left(t,x,y\right)$, $v\left(t,x,y\right)$ в вершинах окрестности $\Omega _{p,q}^{k} $ возможны различные варианты расположения четырехугольника $P_{p,q}^{k-1} $ относительно множеств $\Omega _{\tau }^{k-1} =\left\{t=t_{k-1} \right\}\times \Omega _{\tau } $, $\Omega _{2\tau }^{k-1} =\left\{t=t_{k-1} \right\}\times \Omega _{2\tau } $. Если $P_{p,q}^{k-1} \subset \Omega _{\tau }^{k-1} $, как показано на рис.~\ref{Pic4}(а), то для вычисления $\rho _{p,q}^{h,k} $ используем формулу \eqref{GrindEQ__20_} при $s=1$. Если $P_{p,q}^{k-1} \not\subset \Omega _{\tau }^{k-1} $, то разделим множество $P_{p,q}^{k-1} $ на два множества следующим образом 
%
\begin{equation}\label{GrindEQ__21_}
    \begin{gathered}
    P_{p,q}^{k-1} =P_{\tau ,p,q}^{k-1} \cup P_{2\tau ,p,q}^{k-1},
    \\    
    \text{где} \; P_{\tau ,p,q}^{k-1} =P_{p,q}^{k-1} \cap \Omega _{\tau }^{k-1} \; \text{и} \; P_{2\tau ,p,q}^{k-1} =P_{p,q}^{k-1} \cap \Omega _{2\tau }^{k-1}. 
\end{gathered}
\end{equation}
%
В общем случае каждое из множеств $P_{\tau ,p,q}^{k-1} $, $P_{2\tau ,p,q}^{k-1} $ может быть многоугольным или пустым. В силу теоремы 1 и с учетом \eqref{GrindEQ__12_}, \eqref{GrindEQ__13_}, \eqref{GrindEQ__21_} получим 
% 
\begin{equation} \label{GrindEQ__22_} 
\begin{gathered}
    \int \limits_{\Omega _{p,q}^{k} }\rho \left(t_{k} ,x,y\right)dxdy  = \int \limits_{P_{\tau ,p,q}^{k-1} }\rho \left(t_{k-1} ,x,y\right)dxdy  + \qquad \qquad \qquad \qquad \\ \qquad \qquad  
    + \int \limits_{P_{2\tau ,p,q}^{k-1} }\rho \left(t_{k-1} ,x,y\right)dxdy  +  \int \limits_{V_{p,q}^{k} }f\left(t,x,y\right)dtdxdy +\delta Q_{p,q}^{h,k-1} ,
\end{gathered}
\end{equation} 
%
где $\left|\delta Q_{p,q}^{h,k-1} \right|=O\left(h^{4} \right)$ при условии выполнения равенства \eqref{GrindEQ__4_}. Найдем точки пересечения сторон многоугольника $P_{2\tau ,p,q}^{k-1} $ с поверхностью $S$, образованной на общей границе смежных множеств $D_{\tau } $ и $D_{2\tau } $.
%
\begin{figure}
    \centering
    a) \includegraphics[scale=0.95]{Pic4.png} \\
    б) \includegraphics[scale=0.95]{Pic5.png}
    \caption{ а) Четырехугольник $P_{p,N_1}^{k-1} \subset \Omega_{\tau}^{k-1}$ б) Разложение четырехугольника $P_{p,N_{1+1} }^{k-1}$ на $P_{\tau, p,N_{1+1} }^{k-1}$ и $P_{2\tau, p,N_{1+1} }^{k-1} $}
    \label{Pic4}
\end{figure}
%
Точки пересечения обозначим $B_{n,m}^{h,k-1} =\left(t_{k-1} ,B_{n,m}^{h,k-1,x} ,B_{n,m}^{h,k-1,y} \right)$, $n,m\in 1,...,4$, рис.~\ref{Pic4}~(б). Из множества вершин $B_{d}^{h,k-1} $, $d=1,...,4$, выберем точки, принадлежащие многоугольнику $P_{{\rm 2}\tau ,p,q}^{k-1} $, и обозначим их через $B_{2\tau ,d}^{h,k-1} $. Из множества вершин $A_{d}^{k-1} $, $d=1,...,4$, выберем точки, лежащие на границе многоугольника $P_{2\tau ,p,q}^{k-1} $, и обозначим их через $A_{2\tau ,d}^{k-1} $. На рис.~\ref{Pic4}~(б) такой точкой является вершина $A_{2\tau ,4}^{k-1} $.
Чтобы вычислить интеграл по множеству $P_{2\tau ,p,q}^{k-1} $, применим теорему 1, где в качестве $\Omega _{p,q}^{k} $ используем четырехугольник $P_{2\tau ,p,q}^{k-1} $. Тогда получим
%
\begin{equation} \label{GrindEQ__23_} 
\begin{gathered}
\int \limits_{P_{2\tau ,p,q}^{k-1} }\rho \left(t_{k-1} ,x,y\right)dxdy = 
\qquad \qquad \qquad \qquad \qquad \qquad\qquad \qquad \qquad  \\  \qquad \qquad \qquad
= \int \limits_{Q_{2\tau ,p,q}^{k-2} }\rho \left(t_{k-2} ,x,y\right)dxdy +\int \limits_{V_{2\tau ,p,q}^{k-1} }f\left(t,x,y\right) dtdxdy . 
\end{gathered}
\end{equation} 

Чтобы аппроксимировать многоугольник $Q_{2\tau ,p,q}^{k-2} $ с вершинами в точках $B_{2\tau ,d}^{k-2} =\left(t_{k-2} ,B_{2\tau ,d}^{k-2,x} ,B_{2\tau ,d}^{k-2,y} \right)$, $B_{2\tau ,n,m}^{k-2} =\left(t_{k-2} ,B_{2\tau ,n,m}^{k-2,x} ,B_{2\tau ,n,m}^{k-2,y} \right)$, опустим прямые отрезки, аппроксимирующие траектории движения вершин многоугольника $P_{2\tau ,p,q}^{k-1} $ до плоскости $t=t_{k-2} $, рис.~\ref{Pic5}. При этом, вместо точных значений функций скорости $u\left(t,x,y\right)$, $v\left(t,x,y\right)$ в вершинах четырехугольника $P_{2\tau ,p,q}^{k-1} $ используем их приближенные значения. Так, в узлах $B_{2\tau ,d}^{h,k-1} $, используем аппроксимацию
\begin{equation*}
\begin{split}
    & u\left(t_{k-1} ,B_{2\tau ,d}^{h,k-1,x} ,B_{2\tau ,d}^{h,k-1,y} \right)\approx u\left(t_{k} ,A_{d}^{k,x} ,A_{d}^{k,y} \right),   \\
    & v\left(t_{k-1} ,B_{2\tau ,d}^{h,k-1,x} ,B_{2\tau ,d}^{h,k-1,y} \right)\approx v\left(t_{k} ,A_{d}^{k,x} ,A_{d}^{k,y} \right).
\end{split}  
\end{equation*}
Соответственно координаты вершин $B_{2\tau ,d}^{k-2} $ вычислим приближенно в следующем виде:
\begin{equation} \label{GrindEQ__24_}
\begin{split}
    & B_{2\tau ,d}^{h,k-2,x} =B_{2\tau ,d}^{h,k-1,x} -\tau u\left(t_{k} ,A_{d}^{k,x} ,A_{d}^{k,y} \right),  \\
    & B_{2\tau ,d}^{h,k-2,y} =B_{2\tau ,d}^{h,k-1,y} -\tau v\left(t_{k} ,A_{d}^{k,x} ,A_{d}^{k,y} \right).
\end{split}
\end{equation} 
%
Поскольку $\left|B_{2\tau ,d}^{h,k-1,x} -A_{d}^{k,x} \right|<h/2$ и $\left|B_{2\tau ,d}^{h,k-1,y} -A_{d}^{k,y} \right|<h/2$, то по аналогии с \eqref{GrindEQ__11_} несложно показать, что
%
\begin{figure}
    \centering
    \includegraphics{Pic6-new.png}
    \caption{Формирование многоугольника $P_{2\tau ,p,N_1+1}^{k-2} $}
    \label{Pic5}
\end{figure}
%
%
\begin{equation} \label{GrindEQ__25_} 
\begin{split}
& \left|B_{2\tau ,d}^{k-2,x} -B_{2\tau ,d}^{h,k-2,x} \right|=O\left(\tau^2 + \tau h\right),\\
& \left|B_{2\tau ,d}^{k-2,y} -B_{2\tau ,d}^{h,k-2,y} \right|=O\left(\tau^2 + \tau h\right).  
\end{split}
\end{equation} 
%
Значения функций скорости $u\left(t,x,y\right)$ и $v\left(t,x,y\right)$ в узле $B_{n,m}^{h,k-1} $ аппроксимируем линейной комбинаций значений функций скорости в узлах $A_{l}^{k-1} $ и $A_{r}^{k-1} $, между которыми находится узел $B_{n,m}^{h,k-1} $ 
%
\begin{equation} \label{GrindEQ__26_} 
%\begin{array}{l} 
\begin{split}
	u\left(t_{k-1} ,B_{n,m}^{h,k-1,x} ,B_{n,m}^{h,k-1,y} \right)\approx \alpha _{n,m} u\left(t_{k-1} ,A_{l}^{k-1,x} ,A_{l}^{k-1,y} \right) + \\
	\qquad \qquad \qquad \qquad \qquad \qquad + \left(1-\alpha _{n,m} \right)u\left(t_{k-1} ,A_{r}^{k-1,x} ,A_{r}^{k-1,y} \right), \\ 
	v\left(t_{k-1} ,B_{n,m}^{h,k-1,x} ,B_{n,m}^{h,k-1,y} \right)\approx \alpha _{n,m} v\left(t_{k-1} ,A_{l}^{k-1,x} ,A_{l}^{k-1,y} \right) + \\
	\qquad \qquad \qquad \qquad \qquad \qquad + \left(1-\alpha _{n,m} \right)v\left(t_{k-1} ,A_{r}^{k-1,x} ,A_{r}^{k-1,y} \right), 
%\end{array} 
\end{split}
\end{equation} 
%
где $\alpha _{n,m} =\left(A_{r}^{k-1,x} -B_{n,m}^{h,k-1,x} \right)/h$. Таким образом, координаты вершины $B_{n,m}^{k-2} $ аппроксимируем в виде

\begin{equation*}
    \begin{gathered}      
    B_{n,m}^{h,k-2,x} = B_{n,m}^{h,x} -\tau \left[\alpha _{n,m} u\left(t_{k-1} ,A_{l}^{k-1,x} ,A_{l}^{k-1,y} \right) + \right.
     \qquad \qquad \qquad \qquad \\ \qquad \qquad \qquad  \qquad \qquad \qquad \qquad
    \left. + \left(1-\alpha _{n,m} \right)u\left(t_{k-1} ,A_{r}^{k-1,x} ,A_{r}^{k-1,y} \right)\right], \\ 
    B_{n,m}^{h,k-2,y} = B_{n,m}^{h,y} -\tau \left[\alpha _{n,m} v\left(t_{k-1} ,A_{l}^{k-1,x} ,A_{l}^{k-1,y} \right) + \right. 
    \qquad \qquad \qquad \qquad \\ \qquad \qquad \qquad  \qquad \qquad \qquad \qquad
    \left. + \left(1-\alpha _{n,m} \right)v\left(t_{k-1} ,A_{r}^{k-1,x} ,A_{r}^{k-1,y} \right)\right]. \end{gathered}
\end{equation*}
%
Поскольку аппроксимация \eqref{GrindEQ__26_} имеет погрешность порядка $O\left(h\right)$, то несложно показать, что
%
\begin{equation} 
\label{GrindEQ__27_} 
\begin{split}
\left|B_{n,m}^{h,k-2,x} -B_{n,m}^{k-2,x} \right|=O\left(\tau \left(\tau +h\right)\right), \\
\left|B_{n,m}^{h,k-2,y} -B_{n,m}^{k-2,y} \right|=O\left(\tau \left(\tau +h\right)\right).  
\end{split}
\end{equation} 
%
Отметим, что точка $B_{n,m}^{h,k-2} $ лежит на прямом отрезке, соединяющем узлы $B_{l}^{h,k-2} $ и $B_{r}^{h,k-2} $, где

\begin{equation} \label{GrindEQ__28_}
\begin{split}
    B_{d}^{h,k-2,x} =A_{d}^{k-1,x} -\tau u\left(t_{k-1} ,A_{d}^{k-1,x} ,A_{d}^{k-1,y} \right),  \\
    B_{d}^{h,k-2,y} =A_{d}^{k-1,y} -\tau v\left(t_{k-1} ,A_{d}^{k-1,x} ,A_{d}^{k-1,y} \right). 
\end{split}
\end{equation} 

Если множество вершин $A_{2\tau ,d}^{k-1} $ не пустое, то построим прямые отрезки, аппроксимирующие траектории движения вершин $A_{2\tau ,d}^{k-1} $ до плоскости $t=t_{k-2} $, рис 5. Координаты прямых отрезков в плоскости $t=t_{k-2} $ определим в следующем виде:

\begin{equation*}
\begin{split}
    B_{d}^{h,k-2,x} =A_{2\tau ,d}^{k-1,x} -\tau u\left(t_{k-1} ,A_{2\tau ,d}^{k-1,x} ,A_{2\tau ,d}^{k-1,y} \right),   \\
    B_{d}^{h,k-2,y} =A_{2\tau ,d}^{k-1,y} -\tau v\left(t_{k-1} ,A_{2\tau ,d}^{k-1,x} ,A_{2\tau ,d}^{k-1,y} \right). 
\end{split}
\end{equation*}
%
Так, на рис.~\ref{Pic5}, точка $A_{2\tau ,4}^{k-1} $ лежит между точками $B_{1,2}^{h,k-1} $ и $B_{3,4}^{h,k-1} $, поэтому вычислим координаты точки $B_{4}^{h,k-2} $ по формуле \eqref{GrindEQ__28_}. В итоге в общем случае мы получим точки $B_{2\tau ,d}^{h,k-2} $, $B_{n,m}^{h,k-2} $, $B_{d}^{h,k-2} $, являющиеся вершинами многоугольника $P_{2\tau ,p,q}^{k-2} $. Поскольку расстояние между вершинами $B_{2\tau ,d}^{h,k-2} $, $B_{n,m}^{h,k-2} $, $B_{d}^{h,k-2} $ есть $O\left(h\right)$ или асимптотически меньше, то с учетом \eqref{GrindEQ__25_}, \eqref{GrindEQ__27_} несложно показать, что 
%
\begin{equation} \label{GrindEQ__29_} 
\int \limits_{Q_{2\tau ,p,q}^{k-2} }\rho  \left(t_{k-2} ,x,y\right)dxdy=\int \limits_{P_{2\tau ,p,q}^{k-2} }\rho  \left(t_{k-2} ,x,y\right)dxdy+\delta Q_{2\tau ,p,q}^{h,k-2} ,  
\end{equation} 
%
где 
%
\begin{equation} \label{GrindEQ__30_} 
\left|\delta Q_{2\tau ,p,q}^{h,k-2} \right|=O\left(h\tau \left(\tau +h\right)\right).  
\end{equation} 

Подставим \eqref{GrindEQ__23_} в \eqref{GrindEQ__22_} и рассмотрим сумму интегралов от функции $f\left(t,x,y\right)$. Несложно понять, что

\begin{equation} \label{GrindEQ__31_} 
\int \limits_{V_{p,q}^{k} }f\left(t,x,y\right)dtdxdy +\int \limits_{V_{2\tau ,p,q}^{k-1} }f\left(t,x,y\right)dtdxdy \approx 2\tau f_{p,q}^{k} {\rm meas}\left(\Omega _{p,q}^{k} \right).  
\end{equation} 
%
По аналогии с \eqref{GrindEQ__18_}, \eqref{GrindEQ__19_} можно показать, что такая аппроксимация имеет порядок точности $O\left(h^{4} \right)$ при использовании соотношения \eqref{GrindEQ__4_}. Таким образом, подставляя \eqref{GrindEQ__9_}, \eqref{GrindEQ__23_}, \eqref{GrindEQ__29_} и \eqref{GrindEQ__31_} в \eqref{GrindEQ__22_}, определим следующую формулу:
%
\begin{equation} \label{GrindEQ__32_} 
\begin{split}
\rho _{p,q}^{h,k} =\frac{1}{h^{2} } \int \limits_{P_{\tau ,p,q}^{k-1} }\rho ^{h} \left(t_{k-1} ,x,y\right)dxdy + \frac{1}{h^{2} } \int \limits_{P_{2\tau ,p,q}^{k-2} }\rho ^{h}  \left(t_{k-2} ,x,y\right)&dxdy+ \\
& + s\tau f_{p,q}^{k} .  
\end{split}
\end{equation} 
%
Здесь $s=1$, если $\left(t_{k} ,x_{p} ,y_{q} \right)\in G_{\tau }^{D} $, и $s=2$, если $\left(t_{k} ,x_{p} ,y_{q} \right)\in G_{2\tau }^{D} $. 

Отметим важные свойства описанного алгоритма, которое позволяется добиться выполнения закона сохранения для решения $\rho ^{h} $ в норме пространства $L_{1} $. Первое свойство, попарное пересечение всех областей интегрирования $P_{p,q}^{k-2} $, $P_{2\tau ,p,q}^{k-2} $ является пустым. Так, $\forall \ p,q,i,j = 0,\dots,N : (p,q) \ne (i,j)$ справедливо
%
\begin{equation} \label{GrindEQ__33_}
P_{p,q}^{k-2} \cap P_{i,j}^{k-2} =\emptyset,\quad
        P_{2\tau ,p,q}^{k-2} \cap P_{2\tau ,i,j}^{k-2}  = \emptyset, \quad P_{2\tau ,p,q}^{k-2} \cap P_{i,j}^{k-2} =\emptyset .
\end{equation} 
%
Второе свойство, объединение всех областей интегрирования на слое $t=t_{k-2} $ равно $\Omega $

\begin{equation} \label{GrindEQ__34_} 
\bigcup _{p,q=0}^{N}\left(P_{2\tau ,p,q}^{k-2} \cup P_{p,q}^{k-2} \right) =\Omega .  
\end{equation} 

\section{Закон сохранения и сходимость}

Проверку выполнения закона сохранения выполним при переходе с одного четного слоя по времени на следующий четный слой, поскольку только на этих временных слоях решение $\rho ^{h} $ определено на всем множестве $\Omega =\left[0,1\right]\times \left[0,1\right]$. В этом случае закон сохранения массы с учетом функции правой части имеет следующий вид [1]
%
\begin{equation} \label{GrindEQ__35_} 
\int \limits_{\Omega }\rho \left(t_{2k} ,x,y\right)dxdy =\int \limits_{\Omega }\rho \left(t_{2k-2} ,x,y\right)dxdy +\int \limits_{t_{2k-2} }^{t_{2k} } \int\limits_{\Omega }f\left(t,x,y\right)dxdy dt.  
\end{equation} 
%
Покажем, что \eqref{GrindEQ__35_} справедливо для функции численного решения $\rho ^{h}$ и функции $f^{h}$, где $f^{h} $ -- кусочно-постоянная функция, определенная в следующем виде:
%
\begin{equation}
\begin{gathered}
    \forall \left(t,x,y\right) \in 
    \left\{
        \begin{array}{c} 
        {
            \left[t_{2k-1}, t_{2k} \right] \times \Omega _{p,q}, \; \text{при}  \; \left(t_{2k}, x_{p}, y_{q} \right) \in G_{\tau }^{D},
        } 
        \\
        {
            \left[t_{2k-2}, t_{2k} \right] \times \Omega _{p,q}, \; \text{при} \;  \left(t_{2k}, x_{p}, y_{q} \right) \in G_{2\tau }^{D}
        }
        \end{array}
    \right.   
\\
    f^{h} \left(t,x,y\right) = f^{h} \left(t_{2k}, x_{p}, y_{q} \right).
\end{gathered}
\end{equation}

Рассмотрим формулу \eqref{GrindEQ__20_} для вычисления решения $\rho ^{h} $ на слое $t=t_{2k} $ при $s=1$ и просуммируем её по всем $p,q:\left(x_{p} ,y_{q} \right)\in G_{\tau }^{\Omega } \backslash \Gamma _{h}^{\Omega } $ 
%
\begin{equation} \label{GrindEQ__36_} 
\begin{split}
    \int \limits_{\Omega _{\tau } \backslash \Omega _{\Gamma } }\rho ^{h} \left(t_{2k} ,x,y\right)dxdy =  \sum _{(p,q):G_{\tau }^{\Omega } \backslash \Gamma _{h}^{\Omega } }\int \limits_{P_{p,q}^{2k-1} }&\rho ^{h} \left(t_{2k-1} ,x,y\right)dxdy  +\\
    & + \int \limits_{2k-1}^{2k} \int \limits_{\Omega _{\tau } \backslash \Omega _{\Gamma } }f^{h} dxdydt.    
\end{split}
\end{equation} 
%
Рассмотрим равенство \eqref{GrindEQ__32_} для вычисления $\rho ^{h} $ на слое $t=t_{2k} $ и просуммируем его по всем $p,q:\left(x_{p} ,y_{q} \right)\in \Gamma _{\tau }^{\Omega } $
%
\begin{equation*}
    \begin{gathered}
        \int \limits_{\Omega _{\Gamma } \cap \Omega _{\tau } }\rho ^{h} \left(t_{2k} ,x,y\right)dxdy = \sum _{(p,q):\Gamma _{\tau }^{\Omega } }\int \limits_{P_{\tau ,p,q}^{k-1} }\rho ^{h} \left(t_{2k-1} ,x,y\right)dxdy  + 
        \qquad \qquad \\ \qquad \qquad 
        + \sum _{(p,q):\Gamma _{\tau }^{\Omega } }\int \limits_{P_{2\tau ,p,q}^{k-2} }\rho ^{h}  \left(t_{2k-2} ,x,y\right)dxdy +\int \limits_{2k-1}^{2k}\int \limits_{\Omega _{\Gamma } \cap \Omega _{\tau } }f^{h} dxdy dt,
    \end{gathered}
\end{equation*}
\noindent 
%
а также по всем $p,q:\left(x_{p} ,y_{q} \right)\in \Gamma _{2\tau }^{\Omega } $
%
\begin{equation*}
    \begin{gathered}
        \int \limits_{\Omega _{\Gamma } \cap \Omega _{2\tau } }\rho ^{h} \left(t_{2k} ,x,y\right)dxdy = \sum _{(p,q):\Gamma _{2\tau }^{\Omega } }\int \limits_{P_{\tau ,p,q}^{k-1} }\rho ^{h} \left(t_{2k-1} ,x,y\right)dxdy + 
        \qquad \qquad \\  \qquad
        + \sum _{(p,q):\Gamma _{2\tau }^{\Omega } }\int \limits_{P_{2\tau ,p,q}^{k-2} }\rho ^{h}  \left(t_{2k-2} ,x,y\right)dxdy + 2\int \limits_{2k-1}^{2k}\int \limits_{\Omega _{\Gamma } \cap \Omega _{2\tau } }f^{h} dx dy dt.    
    \end{gathered}
\end{equation*}
%
Сложим два последних равенства и получим 
%
\begin{equation}\label{GrindEQ__37_}
    \begin{gathered}
        \int \limits_{\Omega _{\Gamma } }\rho ^{h} \left(t_{2k} ,x,y\right)dxdy = \sum _{(p,q):\Gamma _{h}^{\Omega } }\int \limits_{P_{\tau ,p,q}^{k-1} }\rho ^{h} \left(t_{2k-1} ,x,y\right)dxdy  + 
        \qquad \qquad \\
        + \sum _{(p,q):\Gamma _{h}^{\Omega } }\int \limits_{P_{2\tau ,p,q}^{k-2} }\rho ^{h}  \left(t_{2k-2} ,x,y\right)dxdy + \int \limits_{2k-1}^{2k} \int \limits_{\Omega _{\Gamma } \cap \Omega _{\tau } }f^{h} dx dy dt + 
        \\ \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad
        + 2\int \limits_{2k-1}^{2k} \int \limits_{\Omega _{\Gamma } \cap \Omega _{2\tau } }f^{h} dx dy dt.
    \end{gathered}
\end{equation}
%
Сложим \eqref{GrindEQ__36_}, \eqref{GrindEQ__37_} и получим
%
\begin{equation}\label{GrindEQ__38_}
    \begin{gathered}
        \int\limits_{\Omega _{\tau } \cup \Omega _{\Gamma } }\rho ^{h} \left(t_{2k} ,x,y\right)dxdy = \sum _{(p,q):G_{\tau }^{\Omega } \backslash \Gamma _{h}^{\Omega } }\int \limits_{P_{p,q}^{2k-1} }\rho ^{h} \left(t_{2k-1} ,x,y\right)dxdy  + 
        \qquad \\
        \sum _{(p,q):\Gamma _{h}^{\Omega } }\int \limits_{P_{\tau ,p,q}^{2k-1} }\rho ^{h} \left(t_{2k-1} ,x,y\right)dxdy  +  \sum _{(p,q):\Gamma _{h}^{\Omega } }\int \limits_{P_{2\tau ,p,q}^{2k-2} }\rho ^{h}  \left(t_{2k-2} ,x,y\right)dxdy + 
        \\ \qquad \qquad \qquad \qquad \qquad
        + \int \limits_{2k-1}^{2k} \int \limits_{\Omega _{\tau } }f^{h} dx dy dt + 2\int \limits_{2k-1}^{2k} \int \limits_{\Omega _{\Gamma } \cap \Omega _{2\tau } }f^{h} dx dy dt.
    \end{gathered}
\end{equation}
%
Отметим, что 
%
\begin{equation} \label{GrindEQ__39_} 
\begin{gathered}
    \sum_{(p,q):G_{\tau }^{\Omega } \backslash \Gamma _{h}^{\Omega } }  \int \limits_{P_{p,q}^{2k-1} }\rho ^{h} \left(t_{2k-1} ,x,y\right)dxdy  + 
    \qquad \qquad \qquad \qquad \qquad \qquad \qquad \\ \qquad \qquad
    +\sum_{(p,q):\Gamma _{h}^{\Omega } } \int \limits_{P_{\tau ,p,q}^{2k-1} }\rho ^{h} \left(t_{2k-1} ,x,y\right)dxdy  =
    \int \limits_{\Omega _{\tau } }\rho ^{h} \left(t_{2k-1} ,x,y\right)dxdy .    
\end{gathered}
\end{equation} 

Рассмотрим \eqref{GrindEQ__20_} для вычисления решения $\rho ^{h} $ на слое по времени $t=t_{2k-1} $ и просуммируем её по всем $p,q:\left(x_{p} ,y_{q} \right)\in G_{\tau }^{\Omega }$, получим 
%
\begin{equation} \label{GrindEQ__40_} 
\begin{gathered}
    \int \limits_{\Omega _{\tau } }\rho ^{h} \left(t_{2k-1} ,x,y\right)dxdy =\sum _{(p,q):G_{\tau }^{\Omega } }\int \limits_{P_{p,q}^{2k-2} }\rho ^{h} \left(t_{2k-2} ,x,y\right)dxdy  +
    \qquad \qquad  \\ \qquad  \qquad  \qquad  \qquad  \qquad \qquad  \qquad  \qquad  \qquad \qquad 
    + \int\limits_{2k-2}^{2k-1} \int \limits_{\Omega _{\tau } }f^{h} dx dy dt.      
\end{gathered}
\end{equation} 
%
Подставим \eqref{GrindEQ__39_}, \eqref{GrindEQ__40_} в \eqref{GrindEQ__38_} и получим
%
\begin{equation} \label{GrindEQ__41_}
    \begin{gathered}
     \int \limits_{\Omega _{\tau } \cup \Omega _{\Gamma } }\rho ^{h} \left(t_{2k} ,x,y\right)dxdy =\sum _{(p,q):G_{\tau }^{\Omega } }\int \limits_{P_{p,q}^{2k-2} }\rho ^{h} \left(t_{2k-2} ,x,y\right)dxdy  + \qquad  \\
     + \int \limits_{2k-2}^{2k-1} \int \limits_{\Omega _{\tau } }f^{h} dx dy dt + \sum _{(p,q):\Gamma _{h}^{\Omega } }\int \limits_{P_{2\tau ,p,q}^{2k-2} }\rho ^{h}  \left(t_{2k-2} ,x,y\right)dxdy + 
     \\ \qquad  \qquad  \qquad \qquad \qquad 
    + \int \limits_{2k-1}^{2k} \int \limits_{\Omega _{\tau } }f^{h} dx dy dt + 2 \int \limits_{2k-1}^{2k} \int \limits_{\Omega _{\Gamma } \cap \Omega _{2\tau } }f^{h} dx dy dt.
    \end{gathered}
\end{equation}

Далее, рассмотрим формулу \eqref{GrindEQ__20_} для вычисления $\rho ^{h} $ на слое по времени $t=t_{2k} $ и просуммируем её по всем $p,q:\left(x_{p} ,y_{q} \right)\in G_{2\tau }^{\Omega } \backslash \Gamma _{h}^{\Omega } $ 
\begin{equation} \label{GrindEQ__42_} 
\begin{split}
    \int \limits_{\Omega _{2\tau } \backslash \Omega _{\Gamma } }\rho ^{h} \left(t_{2k} ,x,y\right)dxdy =\sum _{(p,q):G_{2\tau }^{\Omega } \backslash \Gamma _{h}^{\Omega } }\int \limits_{P_{p,q}^{2k-2} }\rho ^{h} \left(t_{2k-2} ,x,y\right)dxdy  + \\
    + \int \limits_{2k-2}^{2k} \int \limits_{\Omega _{2\tau } \backslash \Omega _{\Gamma } }f^{h} dx dy dt.    
\end{split}
\end{equation} 


\noindent Отметим, рис.~\ref{Pic6}, что справедливо равенство
% 
\begin{equation} \label{GrindEQ__43_} 
\bigcup _{(p,q):G_{\tau }^{\Omega } }P_{p,q}^{k-2}  \cup \bigcup _{(p,q):\Gamma _{h}^{\Omega } }P_{2\tau ,p,q}^{k-2}  \cup \bigcup _{(p,q):G_{2\tau }^{\Omega } \backslash \Gamma _{h}^{\Omega } }P_{p,q}^{k-2}  =\Omega .  
\end{equation} 

\begin{figure}
    \centering
    \includegraphics[scale=0.9]{Pic6a.png}
    \caption{Разложение вычислительной области $\Omega $ на три объединения множеств}
    \label{Pic6}
\end{figure}

Сложим \eqref{GrindEQ__41_} с  \eqref{GrindEQ__42_} и с учетом \eqref{GrindEQ__43_} получим 
%
\begin{equation} \label{GrindEQ__44_} 
\int \limits_{\Omega }\rho ^{h} \left(t_{2k} ,x,y\right)dxdy =\int \limits_{\Omega }\rho ^{h} \left(t_{2k-2} ,x,y\right)dxdy +\int \limits_{2k-2}^{2k} \int \limits_{\Omega }f^{h} dx dy dt.  
\end{equation} 
%
Тождество \eqref{GrindEQ__44_} представляет собой закон сохранения для функции $\rho ^{h} $ при переходе с временного слоя $t=t_{2k-2} $ на слой $t=t_{2k} $ с учетом внутренних источников, описываемых функций $f^h\left(t,x,y\right)$. 

Для проверки сходимости численного решения $\rho ^{h} \left(t,x,y\right)$ к точному решению $\rho \left(t,x,y\right)$ задачи \eqref{GrindEQ__1_}--\eqref{GrindEQ__3_} используем дискретный аналог нормы пространства \textit{L}${}_{1}$
%
\begin{equation} \label{GrindEQ__45_} 
\left\| \rho ^{h} \left(t_{k} ,x,y\right)-\rho \left(t_{k} ,x,y\right)\right\| _{L_{1}^{h} } =\sum _{p,q=0}^{N}{\rm meas}\left(\Omega _{p,q} \right)\left|\rho _{p,q}^{h,k} -\rho _{p,q}^{k} \right| .  
\end{equation} 
%

%Рассмотрим теорему.
\begin{theorem}
    \textit{Для численного решения }$\rho ^{h} $\textit{ задачи }(\textit{1})--(\textit{3})\textit{, с учетом ограничения }\eqref{GrindEQ__4_},\textit{ справедливо}
    %
    \begin{equation} \label{GrindEQ__46_} 
    \left\| \rho ^{h} \left(t_{2k} ,x,y\right)-\rho \left(t_{2k} ,x,y\right)\right\| _{L_{1}^{h} } \le 2k\left(c_{1} h^{2} +c_{2} h^{3} \right),  
    \end{equation} 
    %
    где константы $c_{1} $, $c_{2} $ не зависят от величин $k$, $\tau $ и $h$.
\end{theorem}

\begin{proof}
Докажем теорему методом индукции по \textit{$k$}. В силу начального условия \eqref{GrindEQ__3_} и тождества \eqref{GrindEQ__45_} очевидно, что \eqref{GrindEQ__46_} справедливо для $k=0$. Предположим, что \eqref{GrindEQ__46_} справедливо для $k-1$ и покажем, что \eqref{GrindEQ__46_} верно для \textit{$k$}. 
Для этого рассмотрим равенство из теоремы 1 для окрестности узла $\left(t_{2k-1} ,x_{p} ,y_{q} \right)\in G_{\tau }^{D} $. Подставим \eqref{GrindEQ__9_}, \eqref{GrindEQ__12_}, \eqref{GrindEQ__18_}, \eqref{GrindEQ__19_} в \eqref{GrindEQ__8_} и получим
$\forall \left(t_{2m-1} ,x_{p} ,y_{q} \right)\in G_{\tau }^{D}$%
\begin{equation} \label{GrindEQ__47_} 
\begin{gathered}
    {\rm meas}\left(\Omega _{p,q} \right)\rho _{p,q}^{2k-1} = \int \limits_{P_{p,q}^{2k-2} }\rho \left(t_{2k-2} ,x,y\right)dxdy + 
    \qquad \qquad \qquad \qquad \qquad \\ \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \
    + \tau {\rm meas}\left(\Omega _{p,q} \right)f_{p,q}^{2k-1} +\delta _{p,q}^{2k-1}.  
\end{gathered}
\end{equation} 
%
При использовании равенства \eqref{GrindEQ__4_} и с учетом оценок погрешностей аппроксимации в \eqref{GrindEQ__9_}, \eqref{GrindEQ__14_}, \eqref{GrindEQ__18_} и \eqref{GrindEQ__19_} справедливо
%
\begin{equation} \label{GrindEQ__48_} 
\delta _{p,q}^{2k-1} =\left\{\begin{array}{c} {\delta _{{\rm inner},p,q}^{2m-1} =\tilde{c}_{1} h^{4} +\tilde{c}_{2} h^{5} +O\left(h^{6} \right),{\rm \; \; \; }\forall \left(x_{p} ,y_{q} \right)\in G_{\tau }^{\Omega } \backslash \Gamma ,} \\ {\delta _{{\rm bound},p,q}^{2m-1} =\hat{c}_{1} h^{3} +\hat{c}_{2} h^{4} +O\left(h^{5} \right),{\rm \; \; \; }\forall \left(x_{p} ,y_{q} \right)\in \Gamma .{\rm \; \; \; \; \; \; \; \; }} \end{array}\right.  
\end{equation} 
%
Численное решение $\rho ^{h} $ в узле $\left(t_{2k-1} ,x_{p} ,y_{q} \right)\in G_{\tau }^{D} $ определяется формулой
%
\begin{equation} \label{GrindEQ__49_}
\begin{split}
    \rho _{p,q}^{h,2k-1}  = \frac{1}{{\rm meas}\left(\Omega _{p,q}^{2k-1} \right)} \int \limits_{P_{p,q}^{2k-2} }\rho ^{h} \left(t_{2k-2} ,x,y\right)dxdy 
    + \tau f_{p,q}^{2k-1} .
\end{split}  
\end{equation} 
%
Рассмотрим модуль разности между равенством \eqref{GrindEQ__47_} и формулой \eqref{GrindEQ__49_}, умноженной на ${\rm meas}\left(\Omega _{p,q}^{2k-1} \right)$
%
\begin{equation} \label{GrindEQ__50_} 
\begin{gathered}
    {\rm meas}\left(\Omega _{p,q} \right)\left|\rho _{p,q}^{2k-1} -\rho _{p,q}^{h,2k-1} \right| = 
    \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \\
    = \left|\int \limits_{P_{p,q}^{2k-2} }\rho \left(t_{2k-2} ,x,y\right)dxdy -\int \limits_{P_{p,q}^{2k-2} }\rho ^{h} \left(t_{2k-2} ,x,y\right)dxdy +\delta _{p,q}^{2k-1} \right|\le 
    \\ \qquad \qquad \qquad
    \le \int \limits_{P_{p,q}^{2k-2} }\left|\rho \left(t_{2k-2} ,x,y\right)-\rho ^{h} \left(t_{2k-2} ,x,y\right)\right|dxdy +\left|\delta _{p,q}^{2k-1} \right|.
\end{gathered}
\end{equation} 
%
Просуммируем \eqref{GrindEQ__50_} по всем узлам $\left(t_{2k-1} ,x_{p} ,y_{q} \right)\in G_{\tau }^{D} $ на слое $t=t_{2k-1} $ и получим 
%
\begin{equation} \label{GrindEQ__51_} 
\begin{gathered}
   \sum _{\left(t_{2k-1} ,x_{p} ,y_{q} \right)\in G_{\tau }^{D} }{\rm meas}\left(\Omega _{p,q} \right)\left|\rho _{p,q}^{2k-1} -\rho _{p,q}^{h,2k-1} \right| \le 
   \qquad \qquad \qquad \qquad \qquad \\
   \le \sum _{\left(x_{p} ,y_{q} \right)\in G_{\tau }^{\Omega } }\int \limits_{P_{p,q}^{2k-2} }\left|\rho \left(t_{2k-2} ,x,y\right)-\rho ^{h} \left(t_{2k-2} ,x,y\right)\right|dxdy  + 
   \\ \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad
   + N^{2} \left|\delta _{{\rm inner},p,q}^{2k-1} \right|+4N\left|\delta _{{\rm bound},p,q}^{2k-1} \right|.
\end{gathered}
\end{equation} 

Далее, используем теорему 1 для окрестности узла $\left(t_{2k} ,x_{p} ,y_{q} \right)\in G_{\tau }^{D} \backslash \Gamma _{h}^{D} $. В силу \eqref{GrindEQ__9_}, \eqref{GrindEQ__12_}, \eqref{GrindEQ__14_}, \eqref{GrindEQ__18_} и \eqref{GrindEQ__19_} получим: $\forall\ \left(t_{2k} ,x_{p} ,y_{q} \right)\in G_{\tau }^{D} \backslash \Gamma _{h}^{D}$
%
\begin{equation} \label{GrindEQ__52_} 
    {\rm meas}\left(\Omega _{p,q} \right)\rho _{p,q}^{2k} =\int \limits_{P_{p,q}^{2k-1} }\rho \left(t_{2k-1} ,x,y\right)dxdy 
    + \tau {\rm meas}\left(\Omega _{p,q} \right)f_{p,q}^{2k} +\delta _{p,q}^{\tau ,2k} ,     
\end{equation} 
%
где по аналогии с \eqref{GrindEQ__48_} справедливо 
%
\begin{equation} \label{GrindEQ__53_} 
\delta _{p,q}^{\tau ,2k} 
=
\left\{\begin{array}{c} 
{\delta _{{\rm inner},p,q}^{\tau ,2k} =\tilde{c}_{3} h^{4} +O\left(h^{5} \right),{\rm \; \; \; }\forall \left(x_{p} ,y_{q} \right)\in \left(G_{\tau }^{\Omega } \backslash \Gamma _{h}^{\Omega } \right)\backslash \Gamma ,} 
\\ 
{\delta _{{\rm bound},p,q}^{\tau ,2k} =\hat{c}_{3} h^{3} +O\left(h^{4} \right),{\rm \; \; \; }\forall \left(x_{p} ,y_{q} \right)\in \Gamma .{\rm \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; }} \end{array}\right.  
\end{equation} 
%
Из \eqref{GrindEQ__20_} следует, что $\forall \left(t_{2k} ,x_{p} ,y_{q} \right)\in G_{\tau }^{D} \backslash \Gamma _{h}^{D}$
%
\begin{equation} \label{GrindEQ__54_} 
\begin{split}
    {\rm meas}\left(\Omega _{p,q} \right)\rho _{p,q}^{h,2k} =\int \limits_{P_{p,q}^{2k-1} }\rho ^{h} \left(t_{2k-1} ,x,y\right)dxdy
    + \tau {\rm meas}\left(\Omega _{p,q} \right)f_{p,q}^{2k} ,    
\end{split}
\end{equation} 
%
Рассмотрим модуль разности между \eqref{GrindEQ__53_} и \eqref{GrindEQ__54_}
%
\begin{equation} \label{GrindEQ__55_} 
\begin{gathered}
    {\rm meas} \left(\Omega _{p,q} \right)\left|\rho _{p,q}^{2k} -\rho _{p,q}^{h,2k} \right| =
    \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \\
    = \left|\int \limits_{P_{p,q}^{2k-1} }\rho \left(t_{2k-1} ,x,y\right)-\rho ^{h} \left(t_{2k-1} ,x,y\right)dxdy +\delta _{p,q}^{\tau ,2k} \right|\le 
    \\ \qquad \qquad \qquad
    \le \int \limits_{P_{p,q}^{2k-1} }\left|\rho \left(t_{2k-1} ,x,y\right)-\rho ^{h} \left(t_{2k-1} ,x,y\right)\right|dxdy +\left|\delta _{p,q}^{\tau ,2k} \right|.
\end{gathered}
\end{equation} 
%
Просуммируем \eqref{GrindEQ__55_} по всем по всем узлам $\left(t_{2k} ,x_{p} ,y_{q} \right)\in G_{\tau }^{D} \backslash \Gamma _{h}^{D} $ на слое $t=t_{2k} $ и получим
%
\begin{equation} \label{GrindEQ__56_} 
\begin{gathered}
    \sum _{\left(x_{p} ,y_{q} \right)\in G_{\tau }^{\Omega } \backslash \Gamma _{h}^{\Omega } }{\rm meas}\left(\Omega _{p,q} \right)\left|\rho _{p,q}^{2k} -\rho _{p,q}^{h,2k} \right| \le 
    \qquad \qquad \qquad \qquad \qquad \qquad \\
    \le \sum _{\left(x_{p} ,y_{q} \right)\in G_{\tau }^{\Omega } \backslash \Gamma _{h}^{\Omega } }\int \limits_{P_{p,q}^{2k-1} }\left|\rho \left(t_{2k-1} ,x,y\right)-\rho ^{h} \left(t_{2k-1} ,x,y\right)\right|dxdy  +
    \\ \qquad \qquad \qquad \qquad \qquad \qquad \qquad
    + N^{2} \left|\delta _{{\rm inner},p,q}^{\tau ,2k} \right|+4N\left|\delta _{{\rm bound},p,q}^{\tau ,2k} \right|.
\end{gathered}
\end{equation} 

Пусть $\left(t_{2k} ,x_{p} ,y_{q} \right)\in G_{\tau }^{D} \cap \Gamma _{h}^{D} $. Подставим \eqref{GrindEQ__9_}, \eqref{GrindEQ__23_}, \eqref{GrindEQ__29_}, \eqref{GrindEQ__31_} в \eqref{GrindEQ__22_} и получим
%
\begin{equation} \label{GrindEQ__57_} 
\begin{split}
    h^{2} \rho _{p,q}^{2k} = \left(\int \limits_{P_{\tau ,p,q}^{2k-1} }\rho \left(t_{2k-1} ,x,y\right)dxdy +\int \limits_{P_{2\tau ,p,q}^{2k-2} }\rho  \left(t_{2k-2} ,x,y\right)dxdy\right) + 
    \\
    + s\tau h^{2} f_{p,q}^{2k} +\delta _{p,q}^{\Gamma _{\tau } ,2k} .   
\end{split}
\end{equation} 
%
Здесь $\delta _{p,q}^{\Gamma _{\tau } ,2k} =\tilde{c}_{4} h^{4} +\tilde{c}_{5} h^{5} +O\left(h^{6} \right)$. Численное решение $\rho ^{h} $ в узле $\left(t_{2k} ,x_{p} ,y_{q} \right)\in G_{\tau }^{D} \cap \Gamma _{h}^{D} $ определяется формулой \eqref{GrindEQ__32_}. Рассмотрим модуль разности между \eqref{GrindEQ__57_} и формулой \eqref{GrindEQ__32_}, умноженной на $h^{2} $
%
\begin{equation} \label{GrindEQ__58_} 
\begin{gathered}
    h^{2} \left|\rho _{p,q}^{2k} -\rho _{p,q}^{h,2k} \right|\le \int \limits_{P_{\tau ,p,q}^{2k-1} }\left|\rho \left(t_{2k-1} ,x,y\right)-\rho ^{h} \left(t_{2k-1} ,x,y\right)\right|dxdy + 
    \qquad \\ \qquad \qquad
    +\int \limits_{P_{2\tau ,p,q}^{2k-2} }\left|\rho \left(t_{2k-2} ,x,y\right)-\rho ^{h} \left(t_{2k-2} ,x,y\right)\right|dxdy +\left|\delta _{p,q}^{\Gamma _{\tau } ,2k} \right|.
\end{gathered}
\end{equation} 
%
Просуммируем \eqref{GrindEQ__58_} по всем узлам $\left(t_{2k} ,x_{p} ,y_{q} \right)\in G_{\tau }^{D} \cap \Gamma _{h}^{D} $ и, учитывая, что $h^{2} ={\rm meas}\left(\Omega _{p,q} \right)$, получим 
%
\begin{equation} \label{GrindEQ__59_} 
\begin{gathered}
    \sum _{\left(x_{p} ,y_{q} \right)\in G_{\tau }^{\Omega } \cap \Gamma _{h}^{\Omega } }{\rm meas}\left(\Omega _{p,q} \right) \left|\rho _{p,q}^{2k} -\rho _{p,q}^{h,2k} \right| \le \qquad 
    \qquad \qquad \qquad \qquad \qquad\\
    \le \sum _{\left(x_{p} ,y_{q} \right)\in G_{\tau }^{\Omega } \cap \Gamma _{h}^{\Omega } } \int \limits_{P_{\tau ,p,q}^{2k-1} }\left|\rho \left(t_{2k-1} ,x,y\right)-\rho ^{h} \left(t_{2k-1} ,x,y\right)\right|dxdy  + 
    \\
    +\sum _{\left(x_{p} ,y_{q} \right)\in G_{\tau }^{\Omega } \cap \Gamma _{h}^{\Omega } } \int \limits_{P_{2\tau ,p,q}^{2k-2} }\left|\rho \left(t_{2k-2} ,x,y\right)-\rho ^{h} \left(t_{2k-2} ,x,y\right)\right|dxdy + 
    \\ \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad +N\left|\delta _{p,q}^{\Gamma _{\tau } ,2k} \right|.
\end{gathered}
\end{equation} 

Пусть $\left(t_{2k} ,x_{p} ,y_{q} \right)\in G_{2\tau }^{D} \cap \Gamma _{h}^{D} $. Рассуждая аналогично получим
%
\begin{equation} \label{GrindEQ__60_}
\begin{gathered}
    \sum _{\left(x_{p} ,y_{q} \right)\in G_{2\tau }^{\Omega } \cap \Gamma _{h}^{\Omega } }{\rm meas}\left(\Omega _{p,q} \right)\left|\rho _{p,q}^{2k} -\rho _{p,q}^{h,2k} \right| \le 
    \qquad \qquad \qquad \qquad \qquad  \qquad \\
    \le \sum _{\left(x_{p} ,y_{q} \right)\in G_{2\tau }^{\Omega } \cap \Gamma _{h}^{\Omega } }\int \limits_{P_{\tau ,p,q}^{2k-1} }\left|\rho \left(t_{2k-1} ,x,y\right)-\rho ^{h} \left(t_{2k-1} ,x,y\right)\right|dxdy  + 
    \\
    +\sum _{\left(x_{p} ,y_{q} \right)\in G_{2\tau }^{\Omega } \cap \Gamma _{h}^{\Omega } }\int \limits_{P_{2\tau ,p,q}^{2k-2} }\left|\rho \left(t_{2k-2} ,x,y\right)-\rho ^{h} \left(t_{2k-2} ,x,y\right)\right|dxdy + 
    \\ 
    \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad + N\left|\delta _{p,q}^{\Gamma _{2\tau } ,2k} \right|.
\end{gathered}
\end{equation} 
%
где $\delta _{p,q}^{\Gamma _{2\tau } ,2k} =\tilde{c}_{6} h^{4} +\tilde{c}_{7} h^{5} +O\left(h^{6} \right)$. Пусть $\left(t_{2k} ,x_{p} ,y_{q} \right)\in G_{2\tau }^{D} \backslash \Gamma _{h}^{D} $. Повторяя рассуждения, получим
%
\begin{equation} \label{GrindEQ__61_} 
\begin{gathered}
    \sum _{\left(x_{p} ,y_{q} \right)\in G_{2\tau }^{\Omega } \backslash \Gamma _{h}^{\Omega } }{\rm meas}\left(\Omega _{p,q} \right)\left|\rho _{p,q}^{2k} -\rho _{p,q}^{h,2k} \right| \le 
     \qquad \qquad \qquad \qquad \qquad \qquad\\
    \le \sum _{\left(x_{p} ,y_{q} \right)\in G_{2\tau }^{\Omega } \backslash \Gamma _{h}^{\Omega } }\int \limits_{P_{p,q}^{2k-2} }\left|\rho \left(t_{2k-2} ,x,y\right)-\rho ^{h} \left(t_{2k-2} ,x,y\right)\right|dxdy  + 
    \\ 
    \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad + N^2 \left|\delta _{p,q}^{2\tau ,2k} \right|.
\end{gathered}
\end{equation} 
%
где $\delta _{p,q}^{2\tau ,2k} =\tilde{c}_{8} h^{4} +\tilde{c}_{9} h^{5} +O\left(h^{6} \right)$. Сложим \eqref{GrindEQ__56_}, \eqref{GrindEQ__59_}, \eqref{GrindEQ__60_} и получим
%
\begin{equation} \label{GrindEQ__62_} 
\begin{gathered}
    \sum _{\left(x_{p} ,y_{q} \right)\in G_{\tau }^{\Omega } \cup \Gamma _{h}^{\Omega } }{\rm meas}\left(\Omega _{p,q} \right)\left|\rho _{p,q}^{2k} -\rho _{p,q}^{h,2k} \right| \le 
     \qquad  \qquad  \qquad  \qquad  \qquad  \qquad \\
    \le \sum _{\left(x_{p} ,y_{q} \right)\in G_{\tau }^{\Omega } \backslash \Gamma _{h}^{\Omega } }\int \limits_{P_{p,q}^{2k-1} }\left|\rho \left(t_{2k-1} ,x,y\right)-\rho ^{h} \left(t_{2k-1} ,x,y\right)\right|dxdy  + 
    \\
    +\sum _{\left(x_{p} ,y_{q} \right)\in \Gamma _{h}^{\Omega } }\int \limits_{P_{\tau ,p,q}^{2k-1} }\left|\rho \left(t_{2k-1} ,x,y\right)-\rho ^{h} \left(t_{2k-1} ,x,y\right)\right|dxdy  +
    \\
    + \sum _{\left(x_{p} ,y_{q} \right)\in \Gamma _{h}^{\Omega } }\int \limits_{P_{2\tau ,p,q}^{2k-2} }\left|\rho \left(t_{2k-2} ,x,y\right)-\rho ^{h} \left(t_{2k-2} ,x,y\right)\right|dxdy  + 
    \\  \qquad  \qquad 
    + N^{2} \left|\delta _{{\rm inner},p,q}^{\tau ,2k} \right|+4N\left|\delta _{{\rm bound},p,q}^{\tau ,2k} \right|+N\left|\delta _{p,q}^{\Gamma _{\tau } ,2k} \right|+N\left|\delta _{p,q}^{\Gamma _{2\tau } ,2k} \right|.
\end{gathered}
\end{equation} 
%
Отметим, что 
%
\begin{equation} \label{GrindEQ__63_} 
\bigcup _{\left(x_{p} ,y_{q} \right)\in G_{\tau }^{\Omega } \backslash \Gamma _{h}^{\Omega } }P_{p,q}^{2k-1}  \cup \bigcup _{\left(x_{p} ,y_{q} \right)\in \Gamma _{h}^{\Omega } }P_{\tau ,p,q}^{2m-1}  =\Omega _{\tau } .  
\end{equation} 
%
Тогда \eqref{GrindEQ__62_} сводится к виду
%
\begin{equation} \label{GrindEQ__64_} 
\begin{gathered}
    \sum _{\left(x_{p} ,y_{q} \right)\in G_{\tau }^{\Omega } \cup \Gamma _{h}^{\Omega } }{\rm meas}\left(\Omega _{p,q} \right)\left|\rho _{p,q}^{2k} -\rho _{p,q}^{h,2k} \right| \le 
     \qquad  \qquad \qquad \qquad \qquad \qquad \\
    \le \int \limits_{\Omega _{\tau } }\left|\rho \left(t_{2k-1} ,x,y\right)-\rho ^{h} \left(t_{2k-1} ,x,y\right)\right|dxdy + 
    \\
    +\sum _{\left(x_{p} ,y_{q} \right)\in \Gamma _{h}^{\Omega } }\int \limits_{P_{2\tau ,p,q}^{2k-2} }\left|\rho \left(t_{2k-2} ,x,y\right)-\rho ^{h} \left(t_{2k-2} ,x,y\right)\right|dxdy  + 
    \\  \qquad  \qquad 
    + N^{2} \left|\delta _{{\rm inner},p,q}^{\tau ,2k} \right|+4N\left|\delta _{{\rm bound},p,q}^{\tau ,2k} \right|+N\left|\delta _{p,q}^{\Gamma _{\tau } ,2k} \right|+N\left|\delta _{p,q}^{\Gamma _{2\tau } ,2k} \right|.
\end{gathered}
\end{equation} 
%
Используя \eqref{GrindEQ__9_}, несложно увидеть, что 
% 
\begin{equation} \label{GrindEQ__65_}
\begin{gathered}
    \int \limits_{\Omega _{\tau } }\left|\rho \left(t_{2k-1} ,x,y\right)-\rho ^{h} \left(t_{2k-1} ,x,y\right)\right|dxdy = 
     \qquad \qquad \qquad \qquad \qquad\\
    \sum _{\left(x_{p} ,y_{q} \right)\in G_{\tau }^{\Omega } }\int \limits_{\Omega _{p,q} }\left|\rho \left(t_{2k-1} ,x,y\right)-\rho ^{h} \left(t_{2k-1} ,x_{p} ,y_{q} \right)\right|dxdy  \le 
    \\
    \le \sum _{\left(x_{p} ,y_{q} \right)\in G_{\tau }^{\Omega } }\int \limits_{\Omega _{p,q} }\left|\rho \left(t_{2k-1} ,x_{p} ,y_{q} \right)-\rho ^{h} \left(t_{2k-1} ,x_{p} ,y_{q} \right)\right|dxdy + 
    \\ \qquad \qquad \qquad \qquad \qquad \qquad \qquad
    + N^{2} \left|\varepsilon _{{\rm inner},p,q}^{2k-1} \right| + 4N\left|\varepsilon _{{\rm bound},p,q}^{2k-1} \right| ,
\end{gathered}
\end{equation} 
%
где 
\begin{equation*}
    \begin{split}
        \varepsilon _{{\rm inner},p,q}^{2k-1} =\tilde{c}_{10} h^{4} +\tilde{c}_{11} h^{5} +O\left(h^{6} \right),{\rm \; \; \; }\forall \left(x_{p} ,y_{q} \right)\in G_{\tau }^{\Omega } \backslash \Gamma, \\
        \varepsilon _{{\rm bound},p,q}^{2k-1} =\hat{c}_{10} h^{3} +\hat{c}_{11} h^{4} +O\left(h^{5} \right),{\rm \; \; \; }\forall \left(x_{p} ,y_{q} \right)\in \Gamma .
    \end{split}
\end{equation*}
%
Подставим \eqref{GrindEQ__51_}, \eqref{GrindEQ__65_} в \eqref{GrindEQ__64_} и получим
%
\begin{equation} \label{GrindEQ__66_} 
\begin{gathered}
    \sum _{\left(x_{p} ,y_{q} \right)\in G_{\tau }^{\Omega } \cup \Gamma _{h}^{\Omega } }{\rm meas}\left(\Omega _{p,q} \right)\left|\rho _{p,q}^{2m} -\rho _{p,q}^{h,2m} \right| \le 
     \qquad  \qquad  \qquad  \qquad  \qquad  \qquad \\
    \le \sum _{\left(t_{k-1,} x_{p} ,y_{q} \right)\in G_{\tau }^{D} }\int \limits_{P_{p,q}^{2k-2} }\left|\rho \left(t_{2k-2} ,x,y\right)-\rho ^{h} \left(t_{2k-2} ,x,y\right)\right|dxdy  + 
    \\
    +\sum _{\left(x_{p} ,y_{q} \right)\in \Gamma _{h}^{\Omega } }\int \limits_{P_{2\tau ,p,q}^{2k-2} }\left|\rho \left(t_{2k-2} ,x,y\right)-\rho ^{h} \left(t_{2k-2} ,x,y\right)\right|dxdy  + 
    \\
    + N^{2} \left|\delta _{{\rm inner},p,q}^{\tau ,2k} \right| + 4N\left|\delta _{{\rm bound},p,q}^{\tau ,2k} \right|+ N\left|\delta _{p,q}^{\Gamma _{\tau } ,2k} \right| + N\left|\delta _{p,q}^{\Gamma _{2\tau } ,2k} \right| + 
    \\  \quad
    + N^{2} \left|\varepsilon _{{\rm inner},p,q}^{2k-1} \right|+4N\left|\varepsilon _{{\rm bound},p,q}^{2k-1} \right|+N^{2} \left|\delta _{{\rm inner},p,q}^{2k-1} \right|+4N\left|\delta _{{\rm bound},p,q}^{2k-1} \right|.
\end{gathered}
\end{equation} 
%
Сложим \eqref{GrindEQ__66_} c \eqref{GrindEQ__61_} и получим
\begin{equation} \label{GrindEQ__67_} 
\begin{gathered}
    \sum _{\left(x_{p} ,y_{q} \right)\in G_{h}^{\Omega } }{\rm meas}\left(\Omega _{p,q} \right)\left|\rho _{p,q}^{2k} -\rho _{p,q}^{h,2k} \right| \le 
     \qquad \qquad \qquad \qquad \qquad \qquad \qquad \\
    \le \sum _{\left(t_{2k-1,} x_{p} ,y_{q} \right)\in G_{\tau }^{D} }\int \limits_{P_{p,q}^{2k-2} }\left|\rho \left(t_{2k-2} ,x,y\right)-\rho ^{h} \left(t_{2k-2} ,x,y\right)\right|dxdy  + \\
    + \sum _{\left(t_{2m,} x_{p} ,y_{q} \right)\in \Gamma _{h}^{D} }\int \limits_{P_{2\tau ,p,q}^{2k-2} }\left|\rho \left(t_{2k-2} ,x,y\right)-\rho ^{h} \left(t_{2k-2} ,x,y\right)\right|dxdy  + \\
    + \sum _{\left(t_{2m,} x_{p} ,y_{q} \right)\in G_{2\tau }^{D} \backslash \Gamma _{h}^{D} }\int \limits_{P_{p,q}^{2k-2} }\left|\rho \left(t_{2k-2} ,x,y\right)-\rho ^{h} \left(t_{2k-2} ,x,y\right)\right|dxdy  + \\
    + N^2\left|\delta _{p,q}^{2\tau ,2k} \right|+N^{2} \left|\delta _{{\rm inner},p,q}^{\tau ,2k} \right|+4N\left|\delta _{{\rm bound},p,q}^{\tau ,2k} \right|+N\left|\delta _{p,q}^{\Gamma _{\tau } ,2k} \right|+ 
    \\
    + N\left|\delta _{p,q}^{\Gamma _{2\tau } ,2k} \right|+ 
    N^{2} \left|\varepsilon _{{\rm inner},p,q}^{2k-1} \right|+4N\left|\varepsilon _{{\rm bound},p,q}^{2k-1}\right|+
    N^{2} \left|\delta _{{\rm inner},p,q}^{2k-1} \right| +
    \\
    %&
    \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad +4N\left|\delta _{{\rm bound},p,q}^{2k-1} \right|.
\end{gathered}
\end{equation} 
%
Отметим, что 
%
\begin{equation} \label{GrindEQ__68_} 
\begin{gathered}
\bigcup _{\left(t_{2k-1,} x_{p} ,y_{q} \right)\in G_{\tau }^{D} }P_{p,q}^{2k-2} 
 \quad \cup \quad
\bigcup _{\left(t_{2k,} x_{p} ,y_{q} \right)\in \Gamma _{h}^{D} }P_{2\tau ,p,q}^{2k-2}  
 \quad \cup 
   \qquad \qquad \qquad\\ \qquad \qquad \qquad \qquad \qquad \qquad \qquad
 \cup \quad 
\bigcup _{\left(t_{2k,} x_{p} ,y_{q} \right)\in G_{2\tau }^{D} \backslash \Gamma _{h}^{D} }P_{p,q}^{2k-2}  =\Omega .  
\end{gathered}
\end{equation} 
%
Таким образом, \eqref{GrindEQ__67_} сводится к виду
%
\begin{equation*}
    \begin{gathered}
                \sum _{\left(x_{p} ,y_{q} \right)\in G_{h}^{\Omega } }{\rm meas}\left(\Omega _{p,q} \right)\left|\rho _{p,q}^{2k} -\rho _{p,q}^{h,2k} \right| \le 
        \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \\
        \le \int \limits_{\Omega }\left|\rho \left(t_{2k-2} ,x,y\right)-\rho ^{h} \left(t_{2k-2} ,x,y\right)\right|dxdy + 
        \\
        + N \left|\delta _{p,q}^{2\tau ,2k} \right| + N^{2} \left|\delta _{{\rm inner},p,q}^{\tau ,2k} \right| + 4N\left|\delta _{{\rm bound},p,q}^{\tau ,2k} \right| + N\left|\delta _{p,q}^{\Gamma _{\tau } ,2k} \right| + N\left|\delta _{p,q}^{\Gamma _{2\tau } ,2k} \right| + 
        \\ \qquad \qquad 
        + N^{2} \left|\varepsilon _{{\rm inner},p,q}^{2k-1} \right|+4N\left|\varepsilon _{{\rm bound},p,q}^{2k-1} \right|+N^{2} \left|\delta _{{\rm inner},p,q}^{2k-1} \right|+4N\left|\delta _{{\rm bound},p,q}^{2k-1} \right|.
    \end{gathered}
\end{equation*}
%
Используем неравенство аналогичное \eqref{GrindEQ__65_}, но для слоя $t=t_{2k-2}$, и получим
%
\begin{equation} \label{GrindEQ__69_} 
\begin{gathered}
    \sum _{\left(x_{p} ,y_{q} \right)\in G_{h}^{\Omega } }{\rm meas}\left(\Omega _{p,q} \right)\left|\rho _{p,q}^{2k} -\rho _{p,q}^{h,2k} \right| \le 
    \qquad \qquad \qquad \qquad \qquad \qquad \qquad \\
    \le \sum _{\left(x_{p} ,y_{q} \right)\in G_{h}^{\Omega } }\int \limits_{\Omega _{p,q} }\left|\rho \left(t_{2k-2} ,x_{p} ,y_{q} \right)-\rho ^{h} \left(t_{2k-2} ,x_{p} ,y_{q} \right)\right|dxdy  + 
    \\
    + N^{2} \left|\varepsilon _{{\rm inner},p,q}^{2k-2} \right| + 4N\left|\varepsilon _{{\rm bound},p,q}^{2k-2} \right|+N^{2} \left|\delta _{{\rm inner},p,q}^{\tau ,2k} \right|+4N\left|\delta _{{\rm bound},p,q}^{\tau ,2k} \right|+ 
    \\
    + N\left|\delta _{p,q}^{\Gamma _{\tau } ,2k} \right|+N\left|\delta _{p,q}^{\Gamma _{2\tau } ,2k} \right|+ N^{2} \left|\varepsilon _{{\rm inner},p,q}^{2k-1} \right|+4N\left|\varepsilon _{{\rm bound},p,q}^{2k-1} \right|+
    \\ \qquad \qquad \qquad \qquad \qquad
    + N^{2} \left|\delta _{{\rm inner},p,q}^{2k-1} \right|+4N\left|\delta _{{\rm bound},p,q}^{2k-1} \right|+N\left|\delta _{p,q}^{\Gamma _{2\tau } ,2k} \right|.
\end{gathered}
\end{equation} 
%
где 
\begin{equation*}
  \begin{gathered}
      \varepsilon _{{\rm inner},p,q}^{2k-2} =\tilde{c}_{12} h^{4} +\tilde{c}_{13} h^{5} +O\left(h^{6} \right),{\rm \; \; \; }\forall \left(x_{p} ,y_{q} \right)\in G_{h}^{\Omega } \backslash \Gamma , 
      \\
      \varepsilon _{{\rm bound},p,q}^{2k-2} =\hat{c}_{12} h^{3} +\hat{c}_{13} h^{4} +O\left(h^{5} \right),{\rm \; \; \; }\forall \left(x_{p} ,y_{q} \right)\in \Gamma .
  \end{gathered}  
\end{equation*}
%
Несложно понять, что существуют константы $\tilde{c}_{1}^{2k-2} $ и $\tilde{c}_{2}^{2k-2} $ такие, что 
%
\begin{equation*}
    \begin{gathered}
        2\tilde{c}_{1}^{2k-2} h^{2} +2\tilde{c}_{2}^{2k-2} h^{3} = N^{2} \left|\varepsilon _{{\rm inner},p,q}^{2k-2} \right| + 4N\left|\varepsilon _{{\rm bound},p,q}^{2k-2} \right| + N^{2} \left|\delta _{{\rm inner},p,q}^{\tau ,2k} \right| + 
        \qquad \\
        + 4N\left|\delta _{{\rm bound},p,q}^{\tau ,2k} \right| + N\left|\delta _{p,q}^{\Gamma _{\tau } ,2k} \right| + N\left|\delta _{p,q}^{\Gamma _{2\tau } ,2k} \right| + N^{2} \left|\varepsilon _{{\rm inner},p,q}^{2k-1} \right| + 
        \\ \qquad \qquad
        +4N\left|\varepsilon _{{\rm bound},p,q}^{2k-1} \right|+N^{2} \left|\delta _{{\rm inner},p,q}^{2k-1} \right|+4N\left|\delta _{{\rm bound},p,q}^{2k-1} \right|+N\left|\delta _{p,q}^{\Gamma _{2\tau } ,2k} \right|.
    \end{gathered}
\end{equation*}
%
Пусть $c_{1}^{} =\mathop{\max }\limits_{k=1,...,M/2} \tilde{c}_{1}^{2k-2} $, $c_{2}^{} =\mathop{\max }\limits_{k=1,...,M/2} \tilde{c}_{2}^{2k-2} $. Тогда из \eqref{GrindEQ__69_} следует
%
\begin{equation} \label{GrindEQ__70_} 
\begin{gathered}
    \left\| \rho ^{h} \left(t_{2k} ,x,y\right)-\rho \left(t_{2k} ,x,y\right)\right\| _{L_{1}^{h} } \le 
    \qquad \qquad \qquad \qquad \qquad  \qquad \qquad \\ \qquad \qquad 
     \le \left\| \rho ^{h} \left(t_{2k-2} ,x,y\right)-\rho \left(t_{2k-2} ,x,y\right)\right\| _{L_{1}^{h} } + 
    2\left(c_{1} h^{2} +c_{2} h^{3} \right).
\end{gathered}
\end{equation} 

\noindent Из гипотезы индукции следует, что
%
\begin{equation*}
    \left\| \rho ^{h} \left(t_{2k-2} ,x,y\right)-\rho \left(t_{2k-2} ,x,y\right)\right\| _{L_{1}^{h} } \le \left(2k-2\right)\left(c_{1} h^{2} +c_{2} h^{3} \right).   
\end{equation*}
%
Тогда из \eqref{GrindEQ__70_} получаем \eqref{GrindEQ__46_}. 
\end{proof}

\noindent 

\noindent \textbf{Следствие. }\textit{Для численного решения }$\rho ^{h} $\textit{ задачи }(\textit{1})--(\textit{3})\textit{, с учетом ограничения }(\textit{4}),\textit{ справедливо}

\noindent 
\begin{equation} \label{GrindEQ__71_} 
\left\| \rho ^{h} \left(T,x,y\right)-\rho \left(T,x,y\right)\right\| _{L_{1}^{h} } \le T\left(c_{1}^{T} h+c_{2}^{T} h^{2} \right),  
\end{equation} 
\textbf{}

\noindent \textit{где константы $c_{1}^{T} $, $c_{2}^{T} $ не зависят от $\tau $ и $h$.}

\noindent \textbf{}

\noindent Неравенство \eqref{GrindEQ__71_} доказывает сходимость численного решения $\rho ^{h} \left(t,x,y\right)$ к точному решению $\rho \left(t,x,y\right)$ с первым порядком точности в дискретном аналоге нормы пространства $L_{1} $, определенном по формуле \eqref{GrindEQ__45_}.

\section{Результаты вычислительного эксперимента}

Для проведения первого теста рассмотрим функции
%
\begin{equation} \label{GrindEQ__72_}
\begin{gathered}
     \rho \left(t,x,y\right) = 2 + t \sin \left( 2\,\pi x y \right), \\
    u\left(t,x,y\right) =  t^2 x^2 y^2 \left( 1-x \right)  \left( 1-y \right), \quad
    v\left(t,x,y\right) =  t^2 x y \left( 1-x \right)  \left( 1-y \right).
\end{gathered}
\end{equation} 
%
Используя функции из \eqref{GrindEQ__72_}, определим функции $\rho _{{\rm init}} (x,y)$, $f\left(t,x,y\right)$. Таким образом, получим задачу \eqref{GrindEQ__1_}--\eqref{GrindEQ__3_}. Численное решение $\rho ^{h} \left(t,x,y\right)$ найдем на серии сеток при $N=10\cdot 2^{n} $, $n=0,...,5$. При этом используем тождество $\tau = h/2$. Для оценки погрешности $err_{n} \left(t_{k} \right)$ численного решения используем норму \eqref{GrindEQ__45_}. Порядок сходимости численного решения к точному решению оценивается величиной 
$$
\nonumber
{\rm conv}_{n} =\log _{2} \left(err_{n-1} \left(T\right)/err_{n} \left(T\right)\right).
$$
Результаты расчетов показаны в таблице~\ref{tab:1}. Из таблицы видно, что алгоритм имеет первый порядок сходимости.
\begin{table}[h]
    \centering
    \begin{tabular}{|c|c|c|c|c|} \hline 
    $n$ & $N$ &  $K$ & $err_{n} \left(T\right)$ & ${\rm conv}_{n} $ \\ \hline 
    0   & 10  &   20 & 0,1234 &  \\ \hline 
    1   & 20  &   40 & 0,0603 & 1,0316 \\ \hline 
    2   & 40  &  80  & 0,0295 & 1,0318 \\ \hline 
    3   & 80  &  160 & 0,0146 & 1,0166 \\ \hline 
    4   & 160 &  320 & 0,0073 & 1,0078 \\ \hline 
    5   & 320 &  640 & 0,0036 & 1,0030 \\ \hline 
    \end{tabular}
    \vspace{2mm}
    \caption{О порядке сходимости метода}
    \label{tab:1}    
    \vspace{-5mm}
\end{table}
Чтобы проверить выполнение закона сохранения, используем следующие обозначения:
%
\begin{enumerate}
\item  интеграл на слое по времени $t=t_{2k} $
%
\begin{equation*}
    {\rm mass}_{2k} =\int \limits_{\Omega }\rho ^{h} \left(t_{2k} ,x,y\right)dxdy =\sum _{p,q=0}^{N}{\rm meas}\left(\Omega _{p,q} \right)\rho _{p,q}^{h,2k} ;
\end{equation*}
%
\item  интеграл от функции $f^{h} \left(t,x,y\right)$
%
\begin{equation*}
\begin{gathered}
        {\rm source}_{2k} =\int \limits_{2k-2}^{2k} dt \int \limits_{\Omega }f^{h} dx dy = 
    \sum_{\left(x_p,y_q\right) \in G_{\tau}^{\Omega}} \tau\, {\rm meas}\left(\Omega _{p,q} \right)f_{p,q}^{2k} + 
    \\ \qquad \qquad
    + \sum_{\left(x_p,y_q\right) \in G_{\tau}^{\Omega}} \tau\, {\rm meas}\left(\Omega _{p,q} \right)f_{p,q}^{2k-1} 
    + \sum_{\left(x_p,y_q\right) \in G_{2\tau}^{\Omega}} 2\tau \, {\rm meas}\left(\Omega _{p,q} \right)f_{p,q}^{2k};
\end{gathered}
\end{equation*}
%
\item  погрешность выполнения закона сохранения
%
\begin{equation*}
    \varepsilon _{2k} ={\rm mass}_{2k} -\left({\rm mass}_{2k-2} {\rm +source}_{2k} \right).
\end{equation*}
\end{enumerate}

\noindent В таблице~\ref{tab:2} показаны результаты расчетов для $N=320$. Видно, что закон сохранения массы справедлив для численного решения $\rho ^{h} \left(t,x,y\right)$ при переходе с четного слоя по времени на следующий четный слой с точностью до погрешности вычислений.

\begin{table}[h]
    \centering
    \begin{tabular}{|c|c|c|c|c|} \hline 
$k$ & ${\rm mass}_{2k-2} $ & ${\rm mass}_{k} $ & ${\rm source}_{2k} $ & $\varepsilon _{2m} $ \\ \hline 
0         & -        & 2,0      & -         & -         \\ \hline 
2         & 2,0      & 2,00061   & 0.608984E-3   & 1,585E-12 \\ \hline 
4         & 2,00061  & 2,00122   & 0.608984E-3   & 5,684E-14 \\ \hline 
6         & 2,00122  & 2,00183   & 0.608984E-3   & 9,770E-14 \\ \hline 
8         & 2,00183  & 2,00244   & 0.608984E-3   & 2,309E-14 \\ \hline 
10        & 2,00244  & 2,00304   & 0.608984E-3   & 9,992E-14 \\ \hline 
$\cdots$  & $\cdots$ & $\cdots$  & $\cdots$      & $\cdots$  \\ \hline 
636       & 2,38789  & 2,38850   & 0.608801E-3   & 7,150E-14 \\ \hline 
638       & 2,38850  & 2,38911   & 0.608800E-3   & 9.992E-14 \\ \hline 
640       & 2,38911  & 2,38971   & 0.608799E-3   & 1.297E-13 \\ \hline 
\end{tabular}
    \vspace{2mm}    
    \caption{О выполнении закона сохранения}
    \label{tab:2}
    \vspace{-7mm}    
\end{table}
Для проведения второго теста рассмотрим твердотельное вращение круга с постоянной угловой скоростью вокруг своего центра. Для этого в качестве вычислительной области $\Omega$ рассмотрим квадрат $\Omega = [-2,2] \times [-2,2]$. Положим $T=2$, а числа $N_1$, $N_2$ определим такими, что $N_1 \cdot h =-1,2$; $N_2 \cdot h = 1,2$. Функцию плотности в начальный момент времени зададим в виде
\begin{equation*}
\rho_{{\rm init}}(x,y) = 
\left\{
	\begin{gathered}
		1, \text{\ \rm если } x^2+y^2 \leq 1, \\
		0, \text{\ \rm если } x^2+y^2 > 1. \\
	\end{gathered}
\right.
\end{equation*}
Функции скорости зададим в виде 
\begin{equation}
\label{vel_init_test_2}	
\begin{gathered}
	u(t,x,y) = r\, \sin(\phi),\ v(t,x,y) = r\, \cos(\phi), \text{\ \rm если } x^2+y^2 \leq 1, \\
	u(t,x,y) = v(t,x,y) = 0, \text{\ \rm если } x^2+y^2 > 1, \\
\end{gathered}
\end{equation}
где $r = \sqrt{x^2+y^2}$, а угол $\phi \in [0, 2\pi)$ такой, что $\cos(\phi) = x / r,  \sin(\phi) = y / r$. Положим $f(t,x,y) = 0$. Расчёты проведём на серии сеток с шагами $h = 1/(10 \cdot 2^n)$, $n=0,...,5$, и $\tau = h / 4$.
Решение, вычисленное на сетке при $h = 1/40$ в моменты $t=0,5$ и $t=2$, показано на рис. 7.
%
\begin{figure}
	\centering
	а)
	\includegraphics[scale=0.55]{SecondTest_1.png} 
	\linebreak
	б)
	\includegraphics[scale=0.55]{SecondTest_2.png} \\
	\caption{\centering{Решение $\rho^h$ при моделировании твердотельного вращения с разрывными функциями скорости при а) $t=0,5$ б) $t=2$ }}
	\label{Pic_7}
\end{figure}
%
Из рисунка видно, что погрешность накапливается на границе круга и растёт с течением времени. 
При этом, как видно из таблицы~\ref{tab:3} с данными расчётов при $h=1/160$ и $\tau=h/4$, закон сохранения массы для численного решения $\rho^h$ выполняется с точностью до вычислительной погрешности. 
Результаты расчётов на серии сеток показаны в таблице~\ref{tab:4}, где $N= 1/h$, $K= 1/\tau$.
Из таблицы видно, что сходимость численного решения наблюдается, но сходимость не монотонная. Возникающие пики в численном решении на границе круга появляются в виду разрывности функций скорости. 
%
%
\begin{table}[h]
	\centering
	\begin{tabular}{|c|c|c|} \hline 
		$2k$ & ${\rm mass}_{2k}$ & ${\rm mass}_{2k} - {\rm mass}_{0} $   \\ \hline 
		0         & 3,13703125  & -  \\ \hline 
		2         & 3,13703125  & 1.73461E-12   \\ \hline 
		4         & 3,13703125  & 3,39284E-13   \\ \hline 
		$\cdots$  & $\cdots$ & $\cdots$        \\ \hline 
		1276      & 3,13703125  & 1,77636E-14   \\ \hline 
		1278      & 3,13703125  & 3,01981E-14   \\ \hline 
		1280      & 3,13703125  & 4,57412E-14   \\ \hline 
	\end{tabular}
	\vspace{2mm}    
	\caption{\centering{О выполнении закона сохранения при вращении с разрывными функциями скорости}}
	\label{tab:3}
	\vspace{-5mm}    
\end{table}
%
\begin{table}[h]
	\centering
	\begin{tabular}{|c|c|c|c|c|} \hline 
		$n$ & $N$ &  $K$ & $err_{n} \left(T\right)$ & ${\rm conv}_{n} $ \\ \hline 
		0 & 10 & 40    & 1,1173 &  \\ \hline
		1 & 20 & 80    & 1,0461 & 0,095 \\ \hline
		2 & 40 & 160   & 0,6206 & 0,753 \\ \hline
		3 & 80 & 320   & 0,3865 & 0,683 \\ \hline
		4 & 160 & 640  & 0,1968 & 0,974 \\ \hline
		5 & 320 & 1280 & 0,1125 & 0,807 \\ \hline
	\end{tabular}
	\vspace{2mm}
	\caption{\centering{О порядке сходимости в тесте с разрывными функциями скорости}}
	\label{tab:4}    
	\vspace{-5mm}
\end{table}
%
\begin{figure}
	\centering
	а)
	\includegraphics[scale=0.5]{ThirdTest_1.png} 
	\linebreak
	б)
	\includegraphics[scale=0.5]{ThirdTest_2.png} 
	\\
	\caption{\centering{Решение $\rho^h$ при моделировании твердотельного вращения без разрыва в функциях скорости при а) $t=0,5$ б) $t=2$}}
	\label{Pic_7}
	\vspace{-3mm}
\end{figure}
%

Заменим функции скорости (\ref{vel_init_test_2}) на функции вида
%
\begin{equation}
\label{vel_new_test_2}	
	u(t,x,y) = r\, \sin(\phi),\ v(t,x,y) = r\, \cos(\phi)\ \forall \ (x,y) \in \Omega.
\end{equation}
%
Указанные функции скорости не равны нулю на границы вычислительной области $\Omega$. 
%
%
Тем не менее, мы можем провести вычислительные расчёты до некоторого расчётного времени, поскольку в начальный момент времени искомая функция равна нулю на границе области $\Omega$ и в её близи. В результате вычислений на сетке при $h=1/40$ и $\tau=h/4$ получим решение, показанное на рис. 8 при $t=0,5$ и $t=2$.
%
Из рисунков видно, что численное решение $\rho^h$ размывается ввиду численной вязкости алгоритма. При этом, вычисления на серии сеток показывают сходимость описанного алгоритма с первым порядком точности и выполнение закон сохранения для численного решения. 
%
\section{Дискуссия}
К настоящему моменту алгоритмы из семейства полулагранжевых методов уже широко известны и весьма хорошо изучены. Они обладают четырьмя важными свойствами: сходимость, консервативность, монотонность \cite{16} и возможность использовать шаги по времени большие, чем шаги соответствующие алгебраическому ограничению, вытекающему из условия Куранта-Фрид\-рихса-Леви \cite{25}. Кроме этого, несомненным плюсом этого подхода является наличие теоретического обоснования сходимости. Одним из недостатков этого метода являются высокие вычислительные затраты для двумерных (по пространству) и особенно для трехмерных задач. Цель данной работы состоит в том, чтобы для двумерной задачи показать возможность построения алгоритма, использующего разные шаги по времени в разных частях вычислительной области. При решении некоторых задач такой подход, может позволить сократить вычислительные затраты и, соответственно, уменьшить время требуемое для проведения вычислений. Для одномерных задач такое исследование уже проведено \cite{26}. Основная вычислительная сложность этого метода состоит в аппроксимации интегралов на нижних слоях по времени $t_{k-1}$ и $t_{k-2}$ в узлах сетки близких к стыку двух смежных вычислительных областей. В двумерной задаче, в общем случае, область интегрирования может представлять собой пустое множество, треугольник, четырехугольник, пятиугольник и шестиугольник. С одной стороны, такой большой набор возможных областей интегрирования требует построения еще более вычислительно емких алгоритмов интегрирования. С другой стороны, число узлов для которых необходимо применить описанный усложненный алгоритм, требующий еще больших вычислительных затрат, есть $O(N)$. Поэтому на практике выигрыш от применения данного подхода следует ожидать при решении задач на сетке с числом узлов, превышающим некоторое пороговое значение. 
Нужно отметить, что обобщение этого подхода на трехмерный случай является весьма затруднительном, поскольку области интегрирования на нижних слоях по времени в такой задаче будут иметь сложные формы.

Один из рассмотренных вычислительных экспериментов показывает другой важный недостаток полулагранжевых методов --- наличие искусственной вязкости. Данный недостаток характерен для многих численных методов. Вместе с тем, разрабатываются алгоритмы из семейства полулагранжевых методов, в которых искусственная вязкость существенно ниже \cite{16}. Однако, недостатком такого подхода является подвижность вычислительной сетки в пространстве.

Возможность использовать шаги по времени большие, чем шаги соответствующие алгебраическому ограничению, вытекающему из условия Куранта-Фридрихса-Леви, безусловно являются существенным преимуществом. Однако, ограничение на шаг по времени все же остаётся. Оно вытекает из необходимости корректной аппроксимации области интегрирования на нижнем слое по времени. В одномерной задаче, это ограничение сводится к непересечению линий, аппроксимирующих траектории движения границ окрестности любого узла сетки с текущего (верхнего) слоя по времени на предыдущий (нижний) слой. Итоговое ограничение зависит от способа аппроксимации траектории движения. Часто, каждая траектория движения аппроксимируется прямым отрезком. В этом случае, для одномерной задачи несложно получить такое ограничение на шаг по времени. Однако, для двумерных и, тем более, трехмерных задач строгое теоретическое обоснование ограничения на шаг по времени в зависимости от способа аппроксимации траектории является отдельным сложным исследованием. Интуитивно понятно, что при моделировании, например, ламинарного течения в середине вычислительной области можно использовать шаги по времени большие, чем шаги соответствующие традиционному ограничению, вытекающему из условия Куранта-Фридрихса-Леви. Однако, если мы моделируем течение с вихрями, то использование больших шагов по времени в области наличия вихря не возможно. В данной работе, мы не касались вопроса построения и использования такого ограничения, просто полагая, что оно выполнено. В вычислительных экспериментах, для шага по времени мы использовали ограничение, вытекающее из условия Куранта-Фридрихса-Леви.

\section{Заключение}

В работе показана возможность использовать разные шаги по времени в разных частях вычислительной области при решении двумерного уравнения неразрывности полулагранжевым методом. При этом подтверждена сходимость численного решения к точному решению с первым порядком точности в дискретном аналоге нормы пространства $L_{1} $. Кроме этого доказано выполнение закона сохранения для численного решения. Использование большого шага по времени в части вычислительно области, в которой присутствует ламинарное течение, позволит сократить время расчетов без существенной потери точности. В дальнейшем планируется развить описанный полулагранжевый метод для решения двумерного уравнения конвекции-диффузии.

\bigskip
\begin{thebibliography}{1} % Ссылки на литературу приводятся только на английском языке (даже если у источника нет англоязычной версии).


%\bibitem{1} D.G. Northcott,
%\href{https://doi.org/10.1017/CBO9780511565908}{\it Ideal Theory}, Cambridge University %Press, Cambridge, 1953.

%\bibitem{2} I.I. Ivanov,
%\href{https://doi.org/10.1002/cpa.20282}{\it Some properties of commutative algebras}, Siberian Electronic
%Mathematical Reports, {\bf 1} (2004), 142--143.

%\bibitem{3} A.V. Tetenov, M. Samuel, D.A. Vauilin,
%\href{https://doi.org/10.21538/0134-4889-2017-23-4-281-291}{\it On dendrites generated by polyhedral systems and their ramification points}, Trudy Inst. Mat. Mekh. UrO RAN, \textbf{23}:4 (2017), 281--291.

\bibitem{1} M. Holt, {\it The method of characteristics for steady supersonic rotation flow in three dimensional}, J. Fluid. Mech., {\bf 1} (1956), P.409--422. 

\bibitem{2} V. V. Rusanov, {\it General equations characteristics of gas dynamics}, Journal of Computational Mathematics and Mathematical Physics, {\bf 3} (1963), 508--527. 

\bibitem{3} K. M. Magomedov, {\it Method of characteristics for numerical modelling of gas flow}, Journal of Computational Mathematics and Mathematical Physics, {\bf 6}:2 (1966), 313--325.

\bibitem{4} K. M. Magomedov, A.S. Holodov, {\it About numerical schemes construction for hyperbolic equations based on characteristics relations}, Journal of Computational Mathematics and Mathematical Physics, {\bf 9}:2 (1969), 373--386. 

\bibitem{5} A. Wiin-Nielson, {\it On the application of trajectory methods in numerical forecasting}, Tellus, {\bf 11} (1959), 180--186.

\bibitem{6} O. Pironneau, {\it On the transport-diffusion algorithm and its applications to the Navier-Stokes equations}, Numerische Mathematik, {\bf 38}:3 (1982), 309--332.

\bibitem{7} E. Celledonia, B. K. Kometaa, O. Verdiera, {\it High order semi-Lagrangian methods for the incompressible Navier-Stokes equations}, Journal of Scientific Computing, {\bf 66} (2016), 91--115.

\bibitem{8} E. Dementyeva, E. Karepova, V. Shaidurov, {\it The semi-Lagrangian method for the Navier-Stokes problem for an incompressible fluid}, AIP Conference Proceedings, {\bf 1895} (2017), Article ID 110001.

\bibitem{9} A. Lozovskiy, M. Olshanskii, Y. Vassilevski, {\it A finite element scheme for the numerical solution of the Navier--Stokes/Biot coupled problem}, Russian Journal of Numerical Analysis and Mathematical Modelling, \textbf{37} (2022), 159--174. 

\bibitem{10} T. F. Russell, M. A. Celia, {\it An overview of research on Eulerian--Lagrangian localized adjoint methods (ELLAM)}, Advances in Water Resources, {\bf 25} (2002), 1215--1231

\bibitem{11} T. Arbogast, C. Huang, C. Hung, {\it A fully conservative Eulerian-Lagrangian stream-tube method for advection-diffusion problems}, SIAM Journal on Scientific Computing, {\bf 34}: 4 (2012), 447--478

\bibitem{12} M. Tolstykh, G. Goyman, R. Fadeev, V. Shashkin, {\it Implementation of SL-AV global atmosphere model with 10 km horizontal resolution}, Supercomputing. RuSCDays 2020. Communications in Computer and Information Science, Springer, Cham, {\bf 1331} (2020), 216--225. 

\bibitem{13} V. V. Shashkin, R. Y. Fadeev, M. A. Tolstykh, {\it Simulation of stratospheric processes with the SLAV072L96 atmospheric general circulation model}, Russ. Meteorol. Hydrol, {\bf 48} (2023), 487--499. 

\bibitem{14} N. Koh, M. Mohr, S. Eib, U. Rude, \textit{A massively parallel Eulerian-Lagrangian method for advection-dominated transport in viscous fluids}, SIAM Journal on Scientific Computing, \textbf{44}:3 (2022), 260--285.

\bibitem{15} M. Lentine, J. T. Gretarsson, R. Fedkiw, \textit{An unconditionally stable fully conservative semi-Lagrangian method}, Journal of Computational Physics, \textbf{230} (2011), 2857--2879.

\bibitem{16} V. V. Shaidurov, A. V. Vyatkin, E. V. Kuchunova, \textit{Semi-Lagrangian difference approximations with different stability requirements}, Russian Journal of Numerical Analysis and Mathematical Modelling, \textbf{33}:2 (2018), 123--135.

\bibitem{17} A. Iske, M. Kaser, \textit{Conservative semi-Lagrangian advection on adaptive unstructured meshes}, Numerical Methods for Partial Differential Equations, \textbf{20} (2004), 388--411. 

\bibitem{18} K. M. Terekhov, K. D. Nikitin, M. A. Olshanskii, and Y. V. Vassilevski, \textit{A semi-Lagrangian method on dynamically adapted octree meshes}, Russ. J. Numer. Anal. Math. Modelling, \textbf{30}:6 (2015), 363--380. 

\bibitem{19} M. Tavelli, W. Boscheri, \textit{A high-order parallel Eulerian-Lagrangian algorithm for advection-diffusion problems on unstructured meshes}, Internat. J. Numer. Methods Fluids, \textbf{91}:7 (2019), 332--347.

\bibitem{20} T. Arbogast, W. Wang, \textit{Convergence of a fully conservative volume corrected characteristic method for transport problems}, SIAM J. Numer. Anal, \textbf{48}:3 (2010), 797--823. 

\bibitem{21} C. Huang, T. Arbogast, J. Qiu, \textit{An Eulerian-Lagrangian WENO finite volume scheme for advection problems}, Journal of Computational Physics, \textbf{231}:11 (2012), 4028--4052. 

\bibitem{22} C. Huang, T. Arbogast, \textit{An implicit Eulerian--Lagrangian WENO3 scheme for nonlinear conservation laws}, J. Sci. Comput., \textbf{77 }(2018), 1084--1114.

\bibitem{23} B. K. Kometa, A. Tambue, N. Iqbal, \textit{Semi-Lagrangian discontinuous Galerkin methods for scalar hyperbolic conservation laws}, International Journal for Numerical Methods in Fluids, \textbf{94}:5 (2022), 482--503.

\bibitem{24} A. V. Vyatkin, V. V. Shaydurov, E. V. Kuchunova, \textit{Two Finite Volume Schemes for Advection Equation}, Lecture Notes in Computational Science and Engineering, \textbf{141} (2022), 542--561.

\bibitem{25} S. K. Godunov, V. S. Ryabenkii \textit{Difference scheme (introduction to theory)}, Nauka, Moscow, 1977.

\bibitem{26} A. V. Vyatkin, A. D. Malzev \textit{Conservative semi-Lagrangian method for continuity equation with various time steps in different parts of computational domain}, Siberian Aerospace Journal, \textbf{24}:2 (2022), 218-233.

\end{thebibliography}
\end{document}
