Minor changes

[pde-project.git] / Explicit.tex

\chapter{Явная разностная схема}

\section{Построение явной разностной схемы}

Введем регулярную равномерную сетку, на которой будем рассматривать дискретный аналог задачи \eqref{eq:problem}.
\begin{equation}
\label{eq:grid}
\begin{array}{lll}
y = i h_y, & i = \overline{0, I}, & h_y = l_y / I; \\
z = j h_z, & j = \overline{0, J}, & h_z = l_z / J; \\
t = k h_t, & k = \overline{0, K}, & h_t = T / K.
\end{array}
\end{equation}

Здесь $h_y$, $h_z$, $h_t$ --- шаги сетки по пространственным и временной координате соответственно.

Аппроксимируем модель \eqref{eq:problem} простейшей явной разностной схемой, используя метод замены производных разностными отношениями.

\begin{equation}
  \left\{
    \label{eq:scheme:explicit}
    \begin{array}{l}
      \begin{array}{lr}
        \frac{1}{c^2} \frac{u_{ij}^{k+1} - 2u_{ij}^{k} + u_{ij}^{k-1}}{h_t^2} =
        \frac{u_{i+1j}^k - 2u_{ij}^k + u_{i-1j}^k}{h_y^2} + \frac{u_{ij+1}^k - 2u_{ij}^k + u_{ij-1}^k}{h_z^2},&
        i = \overline{1, I-1},\\
        &j=\overline{1, J-1},\\
        &k=\overline{1,K-1};\\
      \end{array}\\
      \begin{array}{rclrclr}
        u_{0j}^k &=& 0, & u_{Ij}^k &=& 0, &j=\overline{1,J-1},k=\overline{0,K};\\
        %&&&&&&k=\overline{0,K};\\
        u_{i0}^k &=& \sin\frac{\pi i}{I}\sin\frac{2\pi c k h_t}{\lambda}, & u_{iJ}^k &=& 0, &i=\overline{0,I},k=\overline{0,K};\\
        %&&&&&&k=\overline{0,K};\\
        u_{ij}^0 &=& 0, & \frac{u_{ij}^1 - u_{ij}^0}{h_t} &=& 0,&i=\overline{1,I-1},j=\overline{1,J-1}.\\
        %&&&&&&j=\overline{1,J-1}.
      \end{array}
    \end{array}
  \right.
\end{equation}

Здесь $u_{ij}^k$ --- дискретная функция, заданная на сетке \eqref{eq:grid}.

\section{Аппроксимация явной разностной схемой}
\label{sec:approx:explicit}

Рассмотрим общий вид краевой задачи относительно функции $u$.
\begin{equation}
  \label{eq:problem:common}
  Lu = f.
\end{equation}

Здесь $L$~--- линейный оператор, а $f$~--- составной функциональный объект.

Для решения \eqref{eq:problem:common} мы используем сеточную задачу
\begin{equation}
  \label{eq:scheme:common}
  L_hu_h = f_h.
\end{equation}

Поскольку решение $u_h$ не равно дискретному аналогу точного решения $[u]_h$, при подстановке в $[u]_h$ в \eqref{eq:scheme:common} 
в правой части возникнет невязка.
\begin{equation}
  \label{eq:discrepancy:common}
  \delta f_h = L_h[u]_h - f_h.
\end{equation}

Говорят, что разностная схема \eqref{eq:scheme:common} аппроксимирует задачу \eqref{eq:problem:common} на решении $u$, если для любой 
последовательности сгущающихся сеток выполняются условия
\begin{eqnarray}
  \label{eq:approx:1}
  f_h &\xrightarrow{h\rightarrow 0}& f;\\
  \label{eq:approx:2}
  \delta f_h &\xrightarrow{h\rightarrow 0}& 0.
\end{eqnarray}


Покажем, что схема \eqref{eq:scheme:explicit} аппроксимирует задачу \eqref{eq:problem} и найдем порядки аппроксимации. Для этого проверим
выполнение условий аппроксимации.

Рассмотрим нашу конкретную схему \eqref{eq:scheme:explicit}. Если ее представить правую часть в виде составного функционального объекта, получим
\begin{equation}
  f_h =
  \left(
    \begin{array}{l}
      f^1_h\\
      f^2_h\\
      f^3_h\\
      f^4_h\\
      f^5_h\\
      f^6_h\\
      f^7_h\\
    \end{array}
  \right)
  =
  \left(
    \begin{array}{l}
      0\\
      0\\
      0\\
      \sin\frac{\pi i}{I}\sin\frac{2\pi c k h_t}{\lambda}\\
      0\\
      0\\
      0\\
    \end{array}
  \right).
\end{equation}

Очевидно, что \eqref{eq:approx:1} выполняется в силу того, что при построении разностной схемы в качестве $f_h$ мы брали 
дискретный аналог правой части $[f]_h$ исходной задачи. Поэтому необходимо лишь доказать выполнение условия \eqref{eq:approx:2}.

Рассмотрим по очереди каждую часть составной невязки $\delta f_h^i, \forall i=\overline{1,7}$. Затем выберем из полученных невязок
максимальную~--- она и будет характеризовать аппроксимацию схемой \eqref{eq:scheme:explicit} уравнения \eqref{eq:problem}.
\begin{multline}
  \label{eq:discrepancy:1}
  \delta f^1_h = \left\{
    \frac{1}{c^2}\frac{u(y_i, z_j, t_{k+1}) - 2u(y_i, z_j, t_k) + u(y_i, z_j, t_{k-1})}{h_t^2} - \right. \\
    - \frac{u(y_{i+1}, z_j, t_{k}) - 2u(y_i, z_j, t_k) + u(y_{i-1}, z_j, t_{k})}{h_y^2} - \\
    \left. - \frac{u(y_{i}, z_{j+1}, t_{k}) - 2u(y_i, z_j, t_k) + u(y_{i}, z_{j-1}, t_{k})}{h_y^2} \right\}.
\end{multline}

Разложим каждое слагаемое в выражении невязки в окрестностях некоторой точки $(y_i, z_j, t_k)$ для $i = \overline{1, I-1}$, $j = \overline{1, J-1}$, $k = \overline{1, K-1}$. Для простоты введем обозначение 
$u = u(y_i, z_j, t_k)$.
\begin{multline}
  \label{eq:item:1}
    \frac{1}{c^2}\frac{u(y_i, z_j, t_{k+1}) - 2u(y_i, z_j, t_k) + u(y_i, z_j, t_{k-1})}{h_t^2} =\\
    = \frac{1}{c^2h_t^2}
    \left(
      u + u'_th_t + u''_{tt}\frac{h_t^2}{2} + u^{'''}_{ttt}\frac{h_t^3}{6} + u^{IV}_{tttt}\frac{h_t^4}{24} + O(h_t^5)
      - 2u + \right.\\
      \left. + u - u'_{t}h_t + u''_{tt}\frac{h_t^2}{2} - u^{'''}_{ttt}\frac{h_t^3}{6} + u^{IV}_{tttt}\frac{h_t^4}{24} + O(h_t^5)
    \right) =\\
    = \frac{1}{c^2}\frac{u''_{tt}h_t^2 + u^{IV}\frac{h_t^4}{12} + O(h_t^5)}{h_t^2} = \frac{u''_{tt} + u^{IV}_{tttt}\frac{h_t^2}{12} + O(h_t^3)}{c^2}.
\end{multline}

Аналогично разложим и преобразуем два других слагаемых
\begin{equation}
  \label{eq:item:2}
  \frac{u(y_{i+1}, z_j, t_{k}) - 2u(y_i, z_j, t_k) + u(y_{i-1}, z_j, t_{k})}{h_y^2} =
  u''_{yy} + u^{IV}_{yyyy}\frac{h_y^2}{12} + O(h_y^3).
\end{equation}

\begin{equation}
  \label{eq:item:3}
  \frac{u(y_{i}, z_{j+1}, t_{k}) - 2u(y_i, z_j, t_k) + u(y_{i}, z_{j-1}, t_{k})}{h_z^2} =
  u''_{zz} + u^{IV}_{zzzz}\frac{h_z^2}{12} + O(h_z^3).
\end{equation}

Подставив \eqref{eq:item:1}, \eqref{eq:item:2} и \eqref{eq:item:3} в \eqref{eq:discrepancy:1} получим выражение для локальной невязки
\begin{multline}
  \label{eq:discrepancy:1:local}
  \delta f_h^1|_{y_i, z_j, t_k} = \left\{
    \frac{1}{c^2}\left( u''_{tt} + u^{IV}_{tttt}\frac{h_t^2}{12} +O(h_t^3) \right)
    - u''_{yy} - u^{IV}_{yyyy}\frac{h_y^2}{12} + O(h_y^3)-\right.\\\left.
    - u''_{zz} - u^{IV}_{zzzz}\frac{h_z^2}{12} + O(h_z^3)
  \right\}_{y_i, z_j, t_k}.
\end{multline}

Вспомним первое уравнение задачи \eqref{eq:problem} и учтем, что $u''_{tt} = c^2 (u''_{yy} + u''_{zz})$. Перепишем \eqref{eq:discrepancy:1:local}, приведя подобные слагаемые
\begin{multline}
  \label{eq:discrepancy:1:final}
  \delta f_h^1|_{y_i, z_j, t_k} =\\= \left\{
    \frac{1}{c^2}u^{IV}_{tttt}\frac{h_t^2}{12}
    - u^{IV}_{yyyy}\frac{h_y^2}{12}
    - u^{IV}_{zzzz}\frac{h_z^2}{12}
    + O(h_t^3) + O(h_y^3) + O(h_z^3)
  \right\}_{y_i, z_j, t_k}.
\end{multline}

Выделим в \eqref{eq:discrepancy:1:final}  ненулевые слагаемые наименьших порядков относительно шагов дискретизации $h_t$, $h_y$ и $h_z$. 
Порядки этих слагаемых и являются порядками локальной аппроксимации в первом соотношении разностной схемы \eqref{eq:scheme:explicit}. 
В нашем случае, это соотношение дает нам порядок $O(h_t^2, h_y^2, h_z^2)$.

Для всех граничных и начальных условий (кроме начального условия второго порядка) схемы \eqref{eq:scheme:explicit} характерно то, что
их левые и правые части ведут себя в узлах сетки так же, как соответствующие условия задачи \eqref{eq:problem}. Следовательно, погрешность
аппроксимации для этих соотношений нулевая.

Рассмотрим последнее соотношение схемы \eqref{eq:scheme:explicit}, соответствующее начальному условию второго рода задачи 
\eqref{eq:problem}. Видно, что левая часть этого условия получена путем замены производных разностными отношениями, то есть она лишь
 аппроксимирует оригинальное условие. 

Исследуем локальную невязку в узле $(y_i, z_j, t_0)$.
\begin{multline}
  \label{eq:discrepancy:7}
  \delta f^7_h|_{y_i, z_j, t_0} =
  \left\{
    \frac{u(y_i, z_j, t_1) - u(y_i, z_j, t_0)}{h_t}
  \right\}
  =\\=
  \left\{
    \frac{u|_{t_0} + \frac{\partial u}{\partial t}|_{t_0}h_t + \frac{\partial^2 u}{\partial t^2}|_{t_0}\frac{h_t^2}{2} + O(h_t^3)}{h_t}
  \right\}.
\end{multline}

Обозначим $u = u(y_i, z_j, t_0)$ и перепишем \eqref{eq:discrepancy:7} с учетом начальных условий задачи \eqref{eq:problem}
\begin{multline}
  \label{eq:discrepancy:7:final}
  \delta f^7_h|_{y_i, z_j, t_0} =
  \left\{
    \frac{u''_{tt}\frac{h_t^2}{2} + O(h_t^3)}{h_t}
  \right\}_{y_i, z_j, t_0}
  =
  \left\{
    u''_{tt}\frac{h_t}{2} + O(h_t^2)
  \right\}_{y_i, z_j, t_0}
  .
\end{multline}

Как видно из полученного выражения для невязки, начальное условие второго рода имеет порядок аппроксимации $O(h_t)$.

Используя определение нормы в пространстве $F_h$ (см. \cite{samarsky}), запишем выражение для погрешности аппроксимации схемы \eqref{eq:scheme:explicit}.
\begin{equation}
  \label{eq:approx:explicit}
  \Vert \delta f_h \Vert_{f_h} = O(h_t, h_y^2, h_z^2).
\end{equation}

\section{Устойчивость явной разностной схемы}
\label{sec:stable:explicit}

Исследуем \eqref{eq:scheme:explicit} на устойчивость. Считаем, что на правую часть $f_h$ схемы \eqref{eq:scheme:common} подается возмущение, причиной которого могут являться погрешности вычислений, некоторая неточность модели и другое. При этом $L_h u_h$ по сути дела представляет собой отклик нашей дискретной задачи на поданное возмущение.

Задача \eqref{eq:scheme:common} называется \textit{устойчивой}, если 
\begin{itemize}
\item она имеет решение,
\item это решение единственно,
\item порядок отклика не превышает порядка возмущения, т.~е.
\begin{equation}
\label{eq:stability_condition:common}
\left\| u_h \right\|_{U_h} = O \left( \left\| f_h \right\| \right).
\end{equation}
\end{itemize}

\eqref{eq:stability_condition:common} означает, что для некоторой нормы из функционального пространства $U_h$ должна подыскаться такая константа $C$, что
\begin{equation}
\label{eq:stability_condition}
\left\| u_h \right\|_{U_h} \le C \left\| f_h \right\|.
\end{equation}

Это определение устойчивости. Если показать первые два утверждения не составляет труда (видно, что \eqref{eq:scheme:explicit} действительно имеет единственное решение в силу наличия достаточного числа краевых условий), то доказать \eqref{eq:stability_condition} уже нелегко.

Вместо этого мы воспользуемся необходимым признаком Неймана, обоснование которого можно найти, например, в \cite{godunov}. Этот признак утверждает, что если схема устойчива, то найдется такая мелкость сетки по времени $h_1 > 0$ и такая константа $C_1 > 0$, что для всех более мелких сеток и любых $\alpha$
\begin{equation}
\label{eq:neumann}
\left| \lambda(\alpha) \right| \le 1 + C_1 h_t,
\end{equation}
где $\lambda(\alpha)$ --- семейство собственных чисел оператора перехода от одного временного слоя к другому.

Следуя \cite{godunov}, положим
\begin{equation}
\label{eq:neumann:u}
u_{ij}^k = \lambda^k(\alpha, \beta) e^{\hat{i} (\alpha i + \beta j)}, 
\end{equation}
где $\hat{i}^2 = -1$.

Подставим это выражение в рекуррентное соотношение из \eqref{eq:scheme:explicit}.

\begin{multline}
\label{eq:u_to_recurr}
\frac{1}{c^2 h_t^2} \left( \lambda^{k+1} e^{\hat{i} (\alpha i + \beta j)} - 2 \lambda^{k} e^{\hat{i} (\alpha i + \beta j)} + \lambda^{k-1} e^{\hat{i} (\alpha i + \beta j)} \right) = \\
\frac{1}{h_y^2} \left( \lambda^{k} e^{\hat{i} (\alpha (i+1) + \beta j)} - 2 \lambda^{k} e^{\hat{i} (\alpha i + \beta j)} + \lambda^{k} e^{\hat{i} (\alpha (i-1) + \beta j)} \right) + \\
\frac{1}{h_y^2} \left( \lambda^{k} e^{\hat{i} (\alpha i + \beta (j-1))} - 2 \lambda^{k} e^{\hat{i} (\alpha i + \beta j)} + \lambda^{k} e^{\hat{i} (\alpha i + \beta (j-1))} \right)
\end{multline}

Разделим \eqref{eq:u_to_recurr} на $\lambda^{k-1} e^{\hat{i} (\alpha i + \beta j)}$, обозначим
\begin{equation}
\label{eq:gamma}
\gamma_y = c^2 \frac{h_t^2}{h_y^2}, \gamma_z = c^2 \frac{h_t^2}{h_z^2}.
\end{equation}

Получаем квадратное уравнение относительно $\lambda$.

\begin{equation}
\label{eq:lambda_eq:temp}
\lambda^2 - 2 \lambda + 1 = \lambda \gamma_y \left( e^{\hat{i} \alpha i} - 2 + e^{- \hat{i} \alpha i} \right) + \lambda \gamma_z \left( e^{\hat{i} \beta j} - 2 + e^{- \hat{i} \beta j} \right).
\end{equation}

После преобразований \eqref{eq:lambda_eq:temp} перепишется в виде
\begin{equation}
\label{eq:lambda_eq}
\lambda^2 + 2 \lambda \left( 2 \gamma_y \sin^2 \frac{\alpha i}{2} +  2 \gamma_z \sin^2 \frac{\beta j}{2} - 1 \right) + 1 = 0.
\end{equation}

Проанализируем корни уравнения \eqref{eq:lambda_eq}. Для начала заметим, что по теореме Виета произведение этих корней всегда равно $1$. Это никак не зависит от $\alpha$, $\beta$, $\gamma_y$ или $\gamma_z$. Дискриминант такого уравнения
\begin{equation}
\label{eq:explicit:discriminant}
D = 4 \left( 2 \gamma_y \sin^2 \frac{\alpha i}{2} +  2 \gamma_z \sin^2 \frac{\beta j}{2} - 1 \right)^2 - 4.
\end{equation}

Если \eqref{eq:explicit:discriminant} отрицателен, то уравнение \eqref{eq:lambda_eq} имеет пару комплексно-сопряженных корней. Т.~к. их произведение равно $1$, то и по модулю каждое из них тоже равно $1$. В этом случае условие \eqref{eq:neumann} выполняется. Оно также выполняется в случае, если дискриминант \eqref{eq:explicit:discriminant} равен нулю. В этом случае едиснтвенный корень уравнения \eqref{eq:lambda_eq} $\lambda = 1$. Положительный же дискриминант нас не устраивает, потому что в этом случае один из корней получается по модулю больше $1$, а другой --- меньше. Получаем соотношение, выполнение которого гарантирует выполнение \eqref{eq:neumann}.
\begin{equation}
\label{eq:cond}
4 \left( 2 \gamma_y \sin^2 \frac{\alpha i}{2} +  2 \gamma_z \sin^2 \frac{\beta j}{2} - 1 \right)^2 - 4 \le 0.
\end{equation}

Разрешаем \eqref{eq:cond} относительно $\gamma_y$ и $\gamma_z$.
\[\begin{array}{l}
\left( 2 \gamma_y \sin^2 \frac{\alpha i}{2} + 2 \gamma_z \sin^2 \frac{\beta j}{2} - 1 \right)^2 \le 1; \\ \\
-1 \le 2 \gamma_y \sin^2 \frac{\alpha i}{2} + 2 \gamma_z \sin^2 \frac{\beta j}{2} - 1 \le 1; \\ \\
2 \gamma_y \sin^2 \frac{\alpha i}{2} + 2 \gamma_z \sin^2 \frac{\beta j}{2} \le 2; \\ \\
\gamma_y \sin^2 \frac{\alpha i}{2} + \gamma_z \sin^2 \frac{\beta j}{2} \le 1.
\end{array}\]

Учитывая, что квадраты синусов за счёт произвольных $\alpha$ и $\beta$ принимают различные значения на отрезке $[0; 1]$, получаем

\begin{equation}
  \label{eq:explicit:neg_discr}
  \gamma_y + \gamma_z \le 1.
\end{equation}

Из \cite{godunov} известно, что для трехслойных задач с постоянными коэффициентами в разностном уравнении признак Неймана имеет достаточную силу, если нет возмущения в начальных условиях. В этой связи мы можем с высокой долей уверенности утверждать об устойчивости задачи \eqref{eq:scheme:explicit} при условиях \eqref{eq:explicit:neg_discr}.

\section{Моделирование процесса на ЭВМ с помощью явной разностной схемы}
\label{sec:explicit:emulation}

Выразим из разностного уравнения задачи \eqref{eq:scheme:explicit} значение $u_{ij}^{k+1}$, благо это единственное значение функции на новом временном слое, входящее в уравнение.

\begin{equation}
\label{eq:explicit:resolve}
u_{ij}^{k+1} =  c^2 \frac{h_t^2}{h_y^2} \left( u_{i+1j}^k - 2u_{ij}^k + u_{i-1j}^k \right) + c^2 \frac{h_t^2}{h_z^2} \left( u_{ij+1}^k - 2u_{ij}^k + u_{ij-1}^k \right) + 2 u_{ij}^k - u_{ij}^{k-1}.
\end{equation}

Из начальных и краевых условий задачи \eqref{eq:scheme:explicit} нам известны значения функции $u_{ij}^k$ на первых двух временных слоях сетки $k = 0$ и $k = 1$, а также ее значения на граничных слоях $i = 0$, $i = I$, $j = 0$ и $j = J$. Этого достаточно, чтобы по двум известным предыдущим временным слоям полностью определить новый временной слой. Таким образом можно посчитать любое сколь угодно большое количество временных слоев.

Изобразим наглядно результаты моделирования. Будем строить зависимость $E_x$ от $z$, поскольку именно она наиболее интересна и наиболее полно отражает качество решения. Наряду с численным решением для сравнения приведём также график сеточного аналога более точного аналитического решения.

Если намеренно нарушить условие \eqref{eq:explicit:neg_discr}, то схема теряет устойчивость. Это можно наблюдать на рисунке \ref{pic:explicit:1}.

\begin{figure}[!hbtp]
  \centering
  \includegraphics[width=\linewidth]{graph/solution/1}
  \caption{График решения по явной схеме. $I = 50$, $J = 200$, $T = 1 \cdot 10^{-14}$, $h_t = 1 \cdot 10^{-15}$}
  \label{pic:explicit:1}
\end{figure}

Если задать малое число шагов по $z$, то получим устойчивую, но достаточно грубую схему. Изобразим её на рисунке \ref{pic:explicit:2}.

\begin{figure}[!hbtp]
  \centering
  \includegraphics[width=\linewidth]{graph/solution/2}
  \caption{График решения по явной схеме. $I = 50$, $J = 50$, $T = 1 \cdot 10^{-14}$, $h_t = 1 \cdot 10^{-16}$}
  \label{pic:explicit:2}
\end{figure}

При измельчении сетки наблюдается улучшение качества численного решения (см.~рис.~\ref{pic:explicit:3}).

\begin{figure}[!hbtp]
  \centering
  \includegraphics[width=\linewidth]{graph/solution/3}
  \caption{График решения по явной схеме. $I = 50$, $J = 200$, $T = 1 \cdot 10^{-14}$, $h_t = 1 \cdot 10^{-16}$}
  \label{pic:explicit:3}
\end{figure}

Изобразим также развитие волнового процесса во времени на рисунках~\ref{pic:explicit:4}~и~\ref{pic:explicit:5}. Число шагов по $y$ можно уменьшить для ускорения вычислений: это слабо сказывается на качестве решения. В последнем случае можно увидеть, как волна дошла до стенки волновода, отразилась и полученная встречная волна интерферирует с прямой волной.

\begin{figure}[!hbtp]
  \centering
  \includegraphics[width=\linewidth]{graph/solution/4}
  \caption{График решения по явной схеме. $I = 20$, $J = 200$, $T = 4 \cdot 10^{-14}$, $h_t = 1 \cdot 10^{-18}$}
  \label{pic:explicit:4}
\end{figure}

\begin{figure}[!hbtp]
  \centering
  \includegraphics[width=\linewidth]{graph/solution/5}
  \caption{График решения по явной схеме. $I = 10$, $J = 200$, $T = 1 \cdot 10^{-13}$, $h_t = 1 \cdot 10^{-16}$}
  \label{pic:explicit:5}
\end{figure}

\section{Вычислительный эксперимент}
\label{sec:explicit:research}

Проверим полученные результаты с помощью вычислительного эксперимента. Подвергнем исследованию порядки аппроксимации, полученные в разделе \ref{sec:approx:explicit}. Для этого поступим следующим образом. Зафиксируем шаги по двум осям, и будем вычислять решение для различных мелкостей сетки по оставшейся оси. Сравнивая это решение в каждой его точке с аналитическим решением в соответсвующих точках, найдём модуль разницы между ними. Среди таких модулей выберем максимальный. Таким образом, получим зависимость погрешности нашего решения от мелкости сетки по определённой оси.

Для начала проделаем это для временной оси. Фиксируем $I = 50$, $J = 50$. Для верности возьмём $600$ элементов ряда Фурье из аналитического решения. Получаем график, изображённый на рисунке \ref{pic:explicit:research:1}. Видно, что он действительно напоминает прямую, как и должно быть, исходя из \eqref{eq:approx:explicit}.

\begin{figure}[!hbtp]
  \centering
  \includegraphics[width=\linewidth]{graph/research/1}
  \caption{Зависимость погрешности численного решения по явной схеме от мелкости сетки по времени}
  \label{pic:explicit:research:1}
\end{figure}

Теперь фиксируем $h_t = 1 \cdot 10^{-16}$. Получаем график для зависимости погрешности решения от мелкости сетки по оси $y$ (рис.~\ref{pic:explicit:research:2}). Он тоже похож на ожидаемый.

Наконец при аналогичных условиях получаем график для оси $z$ (рис.~\ref{pic:explicit:research:3}). Он никак не похож на ожидаемую параболу. Это объясняется тем, что кроме погрешности численного решения, свой вклад вносит также погрешность аналитического решения, кроме того, при увеличении погрешности большой вклад вносят неквадратичные элементы (старших порядков) из разложения \eqref{eq:approx:explicit}.

\begin{figure}[!hbtp]
  \centering
  \includegraphics[width=\linewidth]{graph/research/2}
  \caption{Зависимость погрешности численного решения по явной схеме от мелкости сетки по оси $y$}
  \label{pic:explicit:research:2}
\end{figure}

\begin{figure}[!hbtp]
  \centering
  \includegraphics[width=\linewidth]{graph/research/4}
  \caption{Зависимость погрешности численного решения по явной схеме от мелкости сетки по оси $z$}
  \label{pic:explicit:research:3}
\end{figure}