% TODO: Fix terminology usage of "якобиан"
\input{vars}
\input{config}
\sloppy

\NewDocumentCommand{\MFArgs}
  {}{x^{(p)}_1,x^{(p)}_2,x^{(p)}_3,\dots,x^{(p)}_n}

\begin{document}
\lstset{language=[11]C++}

%title-page
\include{titlepage}
\thispagestyle{empty}
\clearpage

\tableofcontents
\thispagestyle{empty}
\clearpage
\input{intro}

\chapter{Описание численных методов и возможностей библиотек numpy и scipy}
В данном разделе будут описаны численные методы и их доступные
реализации в библиотеках в частях, посвященным классам задач, которые
они решают.

Стоит учитывать, что scipy основан на numpy, поэтому при рассмотрении
возможностей данных библиотек часто может возникнуть ситуация,
когда искомый функционал содержится только в scipy, либо в numpy
или scipy одновременно.

\section{Численное решение нелинейных уравнений}
При решении некоторых практических задач или проведении исследований
может быть получена математическая модель, которая включает
непрерывную функцию \(F(x), x \in \textbf{R}\), и необходимо определить корни уравнения
\(F(x) = 0\). Если данное уравнение не имеет вид \(ax + b = 0\),
где \(a,b\) -- константы, то оно будет нелинейным.

Для решения нелинейных уравнений существует несколько методов, в данной работе будут рассмотрены итерационные.

Каждый из итерационных методов, перечисленных ниже, соответствует
следующему алгоритму из двух этапов \cite[с. 15]{book:nm-examples}.
\begin{enumerate}
   \item Отыскание приближенного значения корня или содержащего
         его отрезка.
   \item Уточнения значения до некоторой степени точности.
\end{enumerate}
Начальное приближение определяется, исходя из физических соображений
решений похожих задач или графических методов. Если ни один из этих
способов не доступен или не позволяет получить начальное приближение,
удовлетворяющее требованиям, то применяют следующий алгоритм
отыскания начального приближения:
\begin{enumerate}
   \item Производится поиск двух близкорасположенных значений
         \(a\) и \(b\) таких, что \(F(a) \cdot F(b) < 0\), при этом
         \(F(x)\) должна быть всюду определена на отрезке \([a;b]\).
   \item В качестве начального приближения первой итерации принимается
         значение \(x_0 \in [a;b]\), обычно это середина данного
         отрезка.
\end{enumerate}

Так как выполняется условие \(F(a) \cdot F(b) < 0\) и \(F(x)\)
непрерывна, то обязательно найдется такое \(x_k \in (a,b)\), что
\(F(x_k) = 0\) либо \(|F(x_k)| < \varepsilon\), где \(\varepsilon\)
--- погрешность искомого решения.

\subsection{Метод деления отрезка пополам}
Данный метод использует технику поиска решения, похожую на бинарный
поиск.

\subsubsection{Описание метода}
Дано начальное приближение \(x_0 = (a+b)/2\) при
\(F(a) \cdot F(b) < 0\). Для поиска решения уравнения \(x_k\)
применяем следующий алгоритм:
\begin{enumerate}
   \item Рассмотрим отрезки \([a;x_i], [x_i;b]\), \(i = 0 \dots k\)
         --- номер итерации. На первой итерации \(i = 0\).
         \label{list:hls_begin}
   \item Из рассмотренных отрезков берем те, что удовлетворяют условию
         \(F(a) \cdot F(b) < 0\), где \(a,b\) --- границы отрезка.
         \label{list:hls_test}
   \item Для каждого из взятых в п.\ref{list:hls_test} отрезков
         вычисляем их длину \(l\). Если \(l < \varepsilon\),
         тогда дальнейшее выполнение данного алгоритма для данного
         отрезка прекращается. За решение уравнения принимается
         число \((a+b)/2\), округленное с учетом заданной погрешности.
         Если для решения задачи достаточно любого одного решения,
         то работа алгоритма прекращается.
   \item Для каждого из взятых в п.\ref{list:hls_test} отрезков
         устанавливаем значения
         \(a,b,x_{i+1}\). Для левого отрезка эти значения будут равны
         \(a = a,b = x_i,x_{i+1} = (a+x_i)/2\), для правого ---
         \(a = x_i,b = b,x_{i+1} = (b+x_i)/2\).
         \label{list:hls_prepare}
   \item Для каждого из взятых отрезков переходим к
         п.\ref{list:hls_begin}, с увеличением номера итерации
         на \(1\) и установленными относительно взятого
         отрезка значениями из п.\ref{list:hls_prepare}.
\end{enumerate}

После применение метода на заданных входных данных получим множество
решений уравнения \(Ans = \{x, x \in \textbf{R}\}, |Ans| \geq 1\).

\subsubsection{Реализации метода в библиотеках numpy, scipy}
Библиотека scipy содержит функцию \textbf{bisect} из модуля
\textbf{scipy.optimize} \cite{links:scipy_doc}, которая реализует
данный метод.

Функция имеет следующие параметры (задаются в порядке перечисления):
\begin{enumerate}
   \item \(f\) --- function

         Функция Python, возвращающая число.\(f\) должна быть
         непрерывной, а \(f(a)\) и \(f(b)\) должны иметь
         противоположные знаки.
   \item \(a\) --- scalar

         Первый конец интервала \([a;b]\).
   \item \(b\) --- scalar

         Второй конец интервала \([a;b]\).
   \item \(xtol\) --- number, необязательный

         Вычисленный корень \(x0\) будет удовлетворять
         \verb|np.allclose(x, x0,| \verb|atol=xtol,|
         \verb|rtol=rtol)|, где \(x\) --- точный корень.
         Параметр должен быть положительным.
   \item \(rtol\) ---  number, необязательный

         Вычисленный корень \(x0\) будет удовлетворять
         \verb|np.allclose(x, x0,| \verb|atol=xtol, rtol=rtol)|,
         где \(x\) --- точный корень. Параметр не может быть
         меньше значения по умолчанию \verb|4*np.finfo(float).eps|.
   \item \(maxiter\) --- int, необязательный

         Если сходимость не достигается в итерациях \(maxiter\),
         возникает ошибка. Должен быть \(\geq 0\).
   \item \(args\) --- tuple, необязательный

         Содержит дополнительные аргументы для функции \(f\).
         \(f\) вызывается с помощью \verb|apply(f, (x)+args)|.
   \item \(full\_output\) --- bool, необязательный

         Если \(full\_output\) имеет значение \verb|False|,
         возвращается корень. Если \(full\_output\) имеет значение
         \verb|True|, возвращаемое значение равно \verb|(x, r)|, где
         \(x\) --- это корень, а \(r\) --- объект \verb|RootResults|.
   \item \(disp\) --- bool, необязательный

         Если \verb|True|, будет сгенерировано исключение
         \verb|RuntimeError|, если алгоритм не сошелся. В противном
         случае статус сходимости записывается в возвращаемый объект
         \verb|RootResults|.
\end{enumerate}

\subsection{Метод Ньютона (метод касательных)}
Данный итерационный метод позволит найти единственное решение
уравнения (если оно существует) для каждого из выбранных начальных
приближений.
\subsubsection{Описание метода}
На итерации \(i, i = 1\dots k\), строится касательная к кривой
\(y = F(x)\) в точке \((x_i;F(x_i))\), затем находится \(x_{i+1}\)
--- пересечение данной касательной с осью \(Ox\). Процесс продолжается
пока не будет достигнута заданная точность (\(| F ( x_i ) |< \varepsilon\)
или \(| x_{i+1} - x_i | < \varepsilon\)).

\subsubsection{Реализации метода в библиотеках numpy, scipy}
Библиотека scipy содержит функцию \textbf{newton} из модуля
\textbf{scipy.optimize} \cite{links:scipy_doc}, которая реализует
данный метод.

Функция имеет следующие параметры (задаются в порядке перечисления):
\begin{enumerate}
   \item \(f\) --- function

         Функция, корень которой требуется. Это должна быть функция
         одной переменной вида \(f(x,a,b,c \dots)\), где \(a,b,c \dots\)
         --- дополнительные аргументы, которые можно передать в параметре
         \(args\).
   \item \(x0\) --- float, sequence, или ndarray

         Начальная оценка корня, которая должна быть где-то рядом с
         фактическим корнем. Если \(f\) не скалярная, то \(f\) должна
         быть векторизована и возвращать последовательность или массив
         той же формы, что и ее первый аргумент.
   \item \(fprime\) --- callable, необязательный

         Производная функции, когда она доступна и удобна. Если это
         \verb|None| (по умолчанию), то используется метод секущей.
   \item \(args\) --- tuple, необязательный

         Дополнительные аргументы для использования при вызове функции.
   \item \(tol\) --- float, необязательный

         Допустимая погрешность значения корня. Если
         \(y=f(x),y \in \textbf{Z}\), рекомендуется большее значение
         \(tol\), так как и действительная, и мнимая части \(x\)
         вносят вклад в \(|x - x0|\).
   \item \(maxiter\) --- int, необязательный

         Максимальное количество итераций.
   \item \(fprime2\) --- callable, необязательный

         Производная функции второго порядка, если она доступна и удобна.
         Если это \verb|None| (по умолчанию), то используется обычный
         метод Ньютона или метод секущих. Если не \verb|None|, то
         используется метод Галлея.
   \item \(x1\) float, необязательный

         Еще одна оценка корня, которая должна быть где-то рядом с фактическим корнем. Используется, если \(fprime\) не указан.
   \item \(rtol\) ---  float, необязательный

         Допустимое отклонение (относительное) для прерывания работы.
   \item \(full\_output\) --- bool, необязательный

         Если \(full\_output\) имеет значение \verb|False| (по умолчанию),
         возвращается корень. Если \verb|True| и \(x0\) --- скаляр,
         возвращаемое значение равно \verb|(x, r)|, где \(x\) --- это
         корень, а \(r\) --- объект \verb|RootResults|. Если \verb|True|
         и \(x0\) --- не скаляр, возвращаемое значение равно
         \verb|(x, converged, zero_der)|, где:
         \begin{itemize}
            \item converged --- \verb|ndarray| из значений bool. Указывает, какие элементы сошлись успешно.
            \item zero\_der --- \verb|ndarray| из значений bool. Указывает, какие элементы имеют нулевую производную.
         \end{itemize}
   \item \(disp\) --- bool, необязательный

         Если \verb|True| и алгоритм не сошелся, будет сгенерировано
         исключение \verb|RuntimeError|, с сообщением, содержащим
         количество итераций и текущее значение функции. В противном
         случае статус сходимости записывается в возвращаемый объект
         \verb|RootResults|. Игнорируется, если \verb|x0| не является
         скалярным. Примечание: это не имеет ничего общего с
         отображением, однако ключевое слово \verb|disp| нельзя
         переименовать для сохранения обратной совместимости.
\end{enumerate}

\subsection{Метод простой итерации}
Данный метод, как и предыдущий, не позволяет найти несколько решений
уравнения за исполнение алгоритма на единственном начальном приближении.
\subsubsection{Описание метода}
Уравнение \(F(x)=0\) приводим к виду \(x = \varphi(x)\), например
\(x = x-F(x)/M\), где \(M\) --- константа.

Условие сходимости алгоритма: \(0<|\varphi'(x)|<1\). Исходя из него,
\(M\) определяется как
\begin{equation}
   M=1.01 \cdot F'(x_0)
   \label{formula:stab_simpleiter}
\end{equation}
где \(x_0\) ---
начальное приближение.

Таким образом, для итерации \(i, i = 1\dots k,\)
\(x_i = \varphi(x_{i-1})\).

Процесс продолжается, пока не будет достигнута заданная точность, в
данном случае достаточно будет выполнения условия:
\(|x_{i-1}-x_i|<\varepsilon\).
\subsubsection{Реализации метода в библиотеках numpy, scipy}
В библиотеках numpy, scipy не найдено реализаций данного метода.

\section{Методы решения систем линейных алгебраических уравнений}
Система линейных алгебраических (далее, СЛУ) уравнений имеет вид
\begin{eqnarray}
   \left\{
   \begin{aligned}
       & a_{11}x_1 + a_{12}x_2 + a_{1n}x_n = b_1  \\
       & a_{21}x_1 + a_{22}x_2 + a_{2n}x_n = b_2  \\
       & \dots\dots\dots\dots\dots\dots\dots\dots \\
       & a_{n1}x_1 + a_{n2}x_2 + a_{nn}x_n = b_n  \\
   \end{aligned}
   \right.
   \label{formula:SLE}
\end{eqnarray}
СЛУ также представима в матричной форме \(AX=B\), где \(A,X,B\) имеют
следующий вид:
\begin{eqnarray}
   A = \left(
   \begin{aligned}
          & a_{11} \quad a_{12} \quad \dots \quad a_{1n}            \\
          & a_{21} \quad a_{22} \quad \dots \quad a_{2n}            \\
          & \ \vdots \ \qquad \vdots \quad  \ \ddots \quad \ \vdots \\
          & a_{n1} \quad a_{n2} \quad \dots \quad a_{nn}            \\
      \end{aligned}
   \right), B = \left(
   \begin{aligned}
          & b_1      \\
          & b_2      \\
          & \ \vdots \\
          & b_n      \\
      \end{aligned}
   \right), X = \left(
   \begin{aligned}
          & x_1      \\
          & x_2      \\
          & \ \vdots \\
          & x_n      \\
      \end{aligned}
   \right)
   \label{formula:eqn_matrix_system}
\end{eqnarray}

Для решения таких систем существуют прямые и итерационные методы.
Прямые методы (к ним относятся "метод Гаусса", "метод обратной матрицы"
и "метод прогонки") позволяют получить решение за конечное количество
операций, точность которого ограничивается лишь погрешностью округления.
Итерационные методы (среди которых есть методы, такие как
"метод простой итерации" и "метод Зейделя") позволяют получить
приближенное решение с помощью последовательного приближения к точному.

Для итерационных методов необходимо начальное приближение, которое
будет обозначено как \(x^{(0)}\).
\subsection{Метод Гаусса}
\subsubsection{Описание метода}
Для решения СЛУ система (\ref{formula:SLE}) приводится к
треугольному виду (\ref{formula:triag_eqn_system}) с помощью цепочки элементарных преобразований.
\begin{eqnarray}
   \left\{
   \begin{aligned}
       & a'_{11}x_1 + a'_{12}x_2 + a'_{1n}x_n = b'_1 \\
       & 0x_1 + a'_{22}x_2 + a'_{2n}x_n = b'_2       \\
       & \dots\dots\dots\dots\dots\dots\dots.        \\
       & 0x_1 + 0x_2 + a'_{nn}x_n = b'_n             \\
   \end{aligned}
   \right.
   \label{formula:triag_eqn_system}
\end{eqnarray}
Данный процесс называется прямым ходом, а нахождение неизвестных
\(x_n, x_{n-1}, \dots,x_1 \) --- обратным.
\subsubsection{Реализации метода в библиотеках numpy, scipy}
В библиотеке scipy реализован частный случай метода Гаусса ---
LU-разложение \cite[с. 259]{book:levitin}. Для получения решения
СЛУ необходимо задействовать две функции из модуля \textbf{scipy.linalg} \cite{links:scipy_doc}:
\begin{enumerate}
   \item Для получения разложения используется функция
         \textbf{lu\_factor}.
   \item Для совершения обратного хода алгоритма используется
         \textbf{lu\_solve}, которая принимает на вход разложение с
         предыдущего этапа.
\end{enumerate}

Функция \textbf{lu\_factor} имеет следующие параметры (задаются в порядке перечисления):
\begin{enumerate}
   \item \(a\) --- (M, N) array\_like

         Матрица для разложения
   \item \(overwrite\_a\) bool, необязательный

         Следует ли перезаписывать данные в A (может повысить
         производительность). По умолчанию \verb|False|.
   \item \(check\_finite\) bool, необязательный

         Проверять, содержит ли входная матрица только конечные числа.
         Отключение может дать прирост производительности, но может
         привести к проблемам (сбоям, незавершению), если входные данные
         содержат бесконечности или NaN. По умолчанию \verb|True|.
\end{enumerate}

Функция \textbf{lu\_solve} имеет следующие параметры (задаются в порядке перечисления):
\begin{enumerate}
   \item \((lu, piv)\) --- tuple

         Факторизация матрицы коэффициентов a, полученная из
         \textbf{lu\_factor}.
   \item \(b\) --- array

         Правая сторона
   \item \(trans\) --- {0, 1, 2}, необязательный

         Тип системы, которую необходимо решить:

         \begin{tabularx}{0.8\textwidth}{|X|X|}
            \hline \(trans\) & вид системы   \\
            \hline 0         & \(ax = b\)    \\
            \hline 1         & \(a^T x = b\) \\
            \hline 2         & \(a^H x = b\) \\
            \hline
         \end{tabularx}\\

         По умолчанию \verb|0|.
   \item \(overwrite\_b\) --- bool, необязательный

         Следует ли перезаписывать данные в \(b\) (может повысить производительность). По умолчанию \verb|False|.
   \item \(check\_finite\) --- bool, необязательный

         Проверять, содержат ли входные матрицы только конечные числа.
         Отключение может дать прирост производительности, но может
         привести к проблемам (сбоям, незавершению), если входные данные
         содержат бесконечности или NaN. По умолчанию \verb|True|.
\end{enumerate}

\subsection{Метод обратной матрицы}
\subsubsection{Описание метода}
Исходная система (\ref{formula:SLE}) представляется в форме
\(AX=B\), тогда вектор неизвестных переменных \(X\) определяется по
формуле (\ref{formula:inv_m_method}).
\begin{equation}
   X=A^{-1}B
   \label{formula:inv_m_method}
\end{equation}
\subsubsection{Реализации метода в библиотеках numpy, scipy}
Отдельной функции для решения СЛУ не существует, вектор \(X\) можно
найти по формуле (\ref{formula:inv_m_method}). Для получения \(A^{-1}\)
существует функция \textbf{inv} в модуле \textbf{scipy.linalg}
\cite{links:scipy_doc} библиотеки scipy, и функция \textbf{inv} в модуле
\textbf{numpy.linalg} библиотеки numpy. Для перемножения \(A^{-1}\)
и \(B\) в языке Python есть оператор \verb|@|, начиная с версии \(3.5\)
\cite{links:numpy_doc}\cite{links:PEP465}.

В результате для получения решения СЛУ необходимо выполнить выражение
\verb|inv(A) @ B|, используя одну из вышеописанных функций.

Функция \textbf{inv} модуля \textbf{scipy.linalg} имеет следующие
параметры (задаются в порядке перечисления):
\begin{enumerate}
   \item \(a\) --- array\_like

         Квадратная матрица, которую необходимо инвертировать.
   \item \(overwrite\_a\) --- bool, необязательный

         Не запоминать состояние \(a\) (может улучшить
         производительность). По умолчанию \verb|False|.
   \item \(check\_finite\) --- bool, необязательный

         Проверять, содержат ли входные матрицы только конечные числа.
         Отключение может дать прирост производительности, но может
         привести к проблемам (сбоям, незавершению), если входные данные
         содержат бесконечности или NaN. По умолчанию \verb|True|.
\end{enumerate}
\pagebreak

Функция \textbf{inv} модуля \textbf{scipy.linalg} имеет следующие
параметры (задаются в порядке перечисления):
\begin{enumerate}
   \item \(a\) --- аналогичен параметру \(a\) функции из модуля \textbf{scipy.linalg}.
\end{enumerate}

\subsection{Метод прогонки}
Данный метод применяется для решения трехдиагональных матриц \hspace{1mm} вида
\begin{equation}
   \left\{
   \begin{aligned}
       & a_1 x_0 + b_1 x_1 + c_1 x_2 = d_1         \\
       & a_2 x_1 + b_2 x_2 + c_2 x_3 = d_2         \\
       & \dots\dots\dots\dots\dots\dots\dots       \\
       & a_n x_{n-1} + b_n x_n + c_n x_{n+1} = d_n \\
   \end{aligned}
   \right.
   \label{formula:eqn_banded_system}
\end{equation}

Является частным случаем метода Гаусса.
\subsubsection{Описание метода}
После исключения переменных ниже главной диагонали с помощью
элементарных преобразований, в каждом уравнении СЛУ остается
\(\leq 2\) неизвестных. В таком случае, формула обратного хода будет
следующей: \(x_i = U_i x_{i+1}+V_i, i = n,n-1,\dots,1\). После замены
\(i\) на \(i-1\) и подстановке выражения в общий вид уравнения из СЛУ
(\ref{formula:eqn_banded_system}) \(a_i x_{i-1} + b_i x_i + c_i x_{i+1}\)
получим следующее выражение:
\begin{equation}
   x_i = - \frac{c_i}{a_i U_{i-1} + b_i} x_{i+1} +
   \frac{d_i - a_i V_{i-1}}{a_i U_{i-1} + b_i}
\end{equation}
Из которого получим:
\begin{equation}
   U_i = - \frac{c_i}{a_i U_{i-1} + b_i},
   V_i = \frac{d_i - a_i V_{i-1}}{a_i U_{i-1} + b_i}
   \hspace{1cm} i = 1,2,3,\dots,n
\end{equation}
При этом \( c_n = 0; a_1 = 0\), в ходе выполнения алгоритма сначала
вычисляем \(U_i, V_i\), затем \(x_i, i =n,n-1,\dots,1\).
\pagebreak

Данный метод в общем случае не устойчив, за исключением случаев,
когда матрица СЛУ обладает свойством диагонального преобладания
(условие \ref{formula:eqn_diag_dominant}) или она положительно
определенная \cite{links:bhatia}.
\begin{equation}
   \sum_{i \ne j} |a_{ij}| < |a_{ii}|; \qquad i=1,2,3,\dots,n
   \label{formula:eqn_diag_dominant}
\end{equation}

\subsubsection{Реализации метода в библиотеках numpy, scipy}
В scipy для решения СЛУ вида (\ref{formula:eqn_banded_system})
существует две функции в модуле \textbf{scipy.linalg}:
\textbf{solve\_banded} \cite{links:scipy_doc} и
\textbf{solveh\_banded} \cite{links:scipy_doc}.

Различие между ними заключается в том, что \textbf{solve\_banded}
не использует метод прогонки, из-за низкой устойчивости метода в общем
случае, что позволяет найти решение даже если матрица не положительно
определенная или не обладает свойством диагонального преобладания.

\textbf{solveh\_banded} реализует метод прогонки, но авторы библиотеки
указывают, что вводимая матрица должна быть положительно определенной
\cite{links:scipy_doc}.

Функция \textbf{solve\_banded} имеет следующие параметры (задаются в
порядке перечисления):
\begin{enumerate}
   \item \((l, u)\) --- (integer, integer) tuple

         Количество ненулевых нижних и верхних диагоналей.
   \item \(ab\) --- (l + u + 1, M) array\_like

         Ленточная матрица.
   \item \(b\) --- (M,) or (M, K) array\_like

         Правая сторона.
   \item \(overwrite\_ab\) --- bool, необязательный

         Разрешить изменять данные в \(ab\) (может повысить
         производительность). По умолчанию \verb|False|.
   \item \(overwrite\_b\) --- bool, необязательный

         Разрешить изменять данные в \(b\) (может повысить
         производительность). По умолчанию \verb|False|.
   \item \(check\_finite\) --- bool, необязательный

         Проверять, содержат ли входные матрицы только конечные числа.
         Отключение может дать прирост производительности, но может
         привести к проблемам (сбоям, незавершению), если входные данные
         содержат бесконечности или NaN. По умолчанию \verb|True|.
\end{enumerate}

Функция \textbf{solveh\_banded} имеет несколько иной набор параметров
(задаются в порядке перечисления):
\begin{enumerate}
   \item \(ab\) --- (u + 1, M) array\_like

         Ленточная матрица, \(u\) --- число верхних диагоналей.
   \item \(b\) --- (M,) or (M, K) array\_like

         Правая сторона.
   \item \(overwrite\_ab\) --- bool, необязательный

         Разрешить изменять данные в \(ab\) (может повысить
         производительность). По умолчанию \verb|False|.
   \item \(overwrite\_b\) --- bool, необязательный

         Разрешить изменять данные в \(b\) (может повысить
         производительность). По умолчанию \verb|False|.
   \item \(lower\) --- bool, необязательный

         Является ли матрица в нижней форме. (По умолчанию используется
         верхняя форма), то есть \verb|False|.
   \item \(check\_finite\) --- bool, необязательный

         Совпадает с параметром \(check\_finite\) функции
         \textbf{solve\_banded}.
\end{enumerate}

Обе функции принимают матрицу \(ab\) либо в верхней (\textbf{solve\_banded}), либо в нижней форме (\textbf{solveh\_banded} при
включенной опции \(lower\)). Например, для матрицы

\begin{tabular}[htpb]{ccccc}
   5 & 2 & -1 & 0  & 0  \\
   1 & 4 & 2  & -1 & 0  \\
   0 & 1 & 3  & 2  & -1 \\
   0 & 0 & 1  & 2  & 2  \\
   0 & 0 & 0  & 1  & 1  \\
\end{tabular}\\
верхняя форма будет следующей:

\begin{tabular}[htpb]{ccccc}
   0 & 0 & -1 & -1 & -1 \\
   0 & 2 & 2  & 2  & 2  \\
   5 & 4 & 3  & 2  & 1  \\
   1 & 1 & 1  & 1  & 0  \\
\end{tabular}

Так как данная матрица не эрмитова, и, следовательно, не положительно
определенна, описание нижней формы для нее неуместно. Если взять эрмитову
положительно определенную матрицу

\begin{tabular}[htpb]{cccccc}
   4  & 2  & -1 & 0  & 0  & 0  \\
   2  & 5  & 2  & -1 & 0  & 0  \\
   -1 & 2  & 6  & 2  & -1 & 0  \\
   0  & -1 & 2  & 7  & 2  & -1 \\
   0  & 0  & -1 & 2  & 8  & 2  \\
   0  & 0  & 0  & -1 & 2  & 9  \\
\end{tabular}\\
то ее нижняя форма будет следующая:

\begin{tabular}[htpb]{cccccc}
   4  & 5  & 6  & 7  & 8 & 9 \\
   2  & 2  & 2  & 2  & 2 & 0 \\
   -1 & -1 & -1 & -1 & 0 & 0 \\
\end{tabular}

\subsection{Метод простой итерации (метод Якоби)}
\subsubsection{Описание метода}
Для матрицы СЛУ (\ref{formula:eqn_matrix_system}) размеров
\(n \times n\), и начального приближения \(x^{(0)}\) приближенное
решение на итерации \(p, p = 1,2,3,\dots,k\) вычисляется по
следующей формуле:
\begin{equation*}
   x^{(p+1)}_i = \frac{1}{a_{ii}}
   (b_i - \sum_{j=1}^{i-1}a_{ij} x^{(p)}_j
   - \sum_{j=i+1}^{n}a_{ij} x^{(p)}_j);
   \quad i = 1,2,3,\dots,n
\end{equation*}

Метод сходится при \(p \to \infty\), если матрица СЛУ обладает свойством
диагонального преобладания (выполняется условие
\ref{formula:eqn_diag_dominant}). Заданная точность достигается
при выполнении условия:
\begin{equation}
   \max_{i \le i \le n} |x^{(p+1)}_i-x^{(p)}_i| < \varepsilon
   \label{formula:precision_iter_sle}
\end{equation}
\subsubsection{Реализации метода в библиотеках numpy, scipy}
Реализаций данного метода в библиотеках numpy, scipy не найдено.

\subsection{Метод Зейделя}
\subsubsection{Описание метода}
Для матрицы СЛУ (\ref{formula:eqn_matrix_system}) размеров
\(n \times n\), и начального приближения \(x^{(0)}\) приближенное
решение на итерации \(p, p = 1,2,3,\dots,k\) вычисляется по
следующей формуле:
\begin{equation*}
   x^{(p+1)}_i = \frac{1}{a_{ii}}
   (b_i - \sum_{j=1}^{i-1}a_{ij} x^{(p+1)}_j
   - \sum_{j=i+1}^{n}a_{ij} x^{(p)}_j);
   \quad i = 1,2,3,\dots,n
\end{equation*}

Данный метод, в отличие от предыдущего, использует уже найденные
компоненты этой же итерации с м\'eньшим индексом.

Сходимость и точность задаются условиями
(\ref{formula:eqn_diag_dominant}) и (\ref{formula:precision_iter_sle}).
\subsubsection{Реализации метода в библиотеках numpy, scipy}
Реализаций данного алгоритма в библиотеках numpy, scipy не найдено.

\section{Численные методы решения систем нелинейных уравнений}
Система нелинейных уравнений имеет следующий вид:
\begin{equation}
   \begin{aligned}
       & F_1(x_1,x_2,\dots,x_n) = 0     \\
       & F_2(x_1,x_2,\dots,x_n) = 0     \\
       & \dots\dots\dots\dots\dots\dots \\
       & F_n(x_1,x_2,\dots,x_n) = 0     \\
   \end{aligned}
   \label{formula:SnLE}
\end{equation}

Для решения данной системы далее будет описано несколько итерационных
методов. Для каждого их них требуется начальное приближение
\(X^{(0)} = \{x^{(0)}_1,x^{(0)}_2,\dots x^{(0)}_n\}\).

\subsection{Метод простой итерации (метод Якоби) для систем
   нелинейных уравнений}
\subsubsection{Описание метода}
Для каждого из уравнений системы (\ref{formula:SnLE}) составляется
уравнение \(f_i: x_i = x_i - F_i(x)/M_i\), где \(M_i\) определяется из
условия сходимости уравнения (\ref{formula:stab_simpleiter}).
В результате получим систему
\begin{equation}
   \begin{aligned}
       & x_1 = f_1(x_1,x_2,\dots,x_n)   \\
       & x_2 = f_2(x_1,x_2,\dots,x_n)   \\
       & \dots\dots\dots\dots\dots\dots \\
       & x_n = f_n(x_1,x_2,\dots,x_n)   \\
   \end{aligned}
   \label{formula:SnLE-Jacobi-fsys}
\end{equation}

Для получения приближенного решения уравнения применяется формула
\begin{equation}
   x^{(p+1)}_{i} = f_i(x^{(p)}_1,x^{(p)}_2,\dots,x^{(p)}_n);
   \quad i=1,2,3,\dots,n
\end{equation}
где \(p, p = 1,2,3,\dots,k\) --- номер итерации.

Заданная точность \(\varepsilon\) достигается выполнением
следующего условия:
\begin{equation*}
   \max_i |x^{(p+1)}_i - x^{(p)}_i | < \varepsilon; \quad
   i = 1,2,3,\dots,n
\end{equation*}
\subsubsection{Реализации метода в библиотеках numpy, scipy}
Реализаций данного метода в библиотеках numpy, scipy не найдено.

\subsection{Метод Зейделя для систем нелинейных уравнений}
Данный метод является модификацией предыдущего; отличие состоит в
условии сходимости и формуле получения приближенного решения.

Ниже будут описаны только вышеперечисленные различия.
\subsubsection{Описание метода}
Формула вычисления приближенного решения данного алгоритма следующая:
\begin{equation*}
   \begin{aligned}
       & x^{(p+1)}_{1} = f_1(x^{(p)}_1,x^{(p)}_2,\dots,x^{(p)}_n) \\
       & x^{(p+1)}_{2} = f_2(x^{(p)}_1,x^{(p)}_2,\dots,x^{(p)}_n) \\
       & \dots\dots\dots\dots\dots\dots\dots\dots\dots            \\
       & x^{(p+1)}_{n} = f_n(x^{(p)}_1,x^{(p)}_2,\dots,x^{(p)}_n) \\
   \end{aligned}
\end{equation*}
где \(p, p = 1,2,3,\dots,k\) --- номер итерации.

Данный алгоритм более требователен к точности начального приближения.

Сходимость метода зависит от характера функций исходной системы,
для определения которого необходимо вычислить значения матрицы
\begin{equation}
   F' = \left(
   \begin{aligned}
       & f^{'}_{11} \ \ f^{'}_{12} \ \ f^{'}_{13} \ \ \dots \ \  f^{'}_{1n} \\
       & f^{'}_{21} \ \ f^{'}_{22} \ \ f^{'}_{23} \ \ \dots \ \  f^{'}_{2n} \\
       & \dots      \ \      \dots \ \      \dots \ \ \dots \ \      \dots  \\
       & f^{'}_{n1} \ \ f^{'}_{n2} \ \ f^{'}_{n3} \ \ \dots \ \ f^{'}_{nn}  \\
   \end{aligned}
   \right)
\end{equation}
где \(f^{'}_{ij} = \frac{\partial f_i}{\partial x_j}\).

Сходимость метода обеспечивается выполнением следующего условия:
\vspace{-5mm}
\begin{equation*}
   |f^{'}_{i1}| + |f^{'}_{i2}| + |f^{'}_{i3}| + \dots |f^{'}_{in}| < 1;
   \quad i = 1,2,3,\dots,n
\end{equation*}

\subsubsection{Реализации метода в библиотеках numpy, scipy}
Реализаций данного метода в библиотеках numpy, scipy не обнаружено.

\subsection{Метод Ньютона решения систем нелинейных уравнений}
\subsubsection{Описание метода}
На каждой итерации \(p, p = 1,2,3,\dots,k\) вычисляются значения
\(\Delta x^{(p)}_1,\) \( \Delta x^{(p)}_2, \Delta x^{(p)}_3, \dots,
\Delta x^{(p)}_n \).

Для этого исходная система уравнений раскладывается в ряд Тейлора по
\(\Delta x^{(p)}_1, \Delta x^{(p)}_2, \Delta x^{(p)}_3, \dots,
\Delta x^{(p)}_n \). Сохранив линейные по данным значениям части, получим
СЛУ:

\begin{equation}
   \begin{aligned}
       & \splitdfrac{\frac{\partial F_1(\MFArgs)}{\partial x_1} \Delta x^{(p)}_1 +
         \frac{\partial F_1(\MFArgs)}{\partial x_2} \Delta x^{(p)}_2 + \dots +}
      {\frac{\partial F_1(\MFArgs)}{\partial x_n} \Delta x^{(p)}_n = -F_1(\MFArgs)}  \\
       & \splitdfrac{\frac{\partial F_2(\MFArgs)}{\partial x_1} \Delta x^{(p)}_1 +
         \frac{\partial F_2(\MFArgs)}{\partial x_2} \Delta x^{(p)}_2 + \dots +}
      {\frac{\partial F_2(\MFArgs)}{\partial x_n} \Delta x^{(p)}_n = -F_2(\MFArgs)}  \\
       & \dots\dots\dots\dots\dots\dots\dots\dots\dots\dots\dots\dots\dots\dots\dots \\
   \end{aligned}
   \label{formula:SnLE-Newton-sys}
\end{equation}
\begin{equation*}
   \begin{aligned}
       & \splitdfrac{\frac{\partial F_n(\MFArgs)}{\partial x_1} \Delta x^{(p)}_1 +
         \frac{\partial F_n(\MFArgs)}{\partial x_2} \Delta x^{(p)}_2 + \dots +}
      {\frac{\partial F_n(\MFArgs)}{\partial x_n} \Delta x^{(p)}_n = -F_n(\MFArgs)} \\
   \end{aligned}
\end{equation*}
Данную систему можно решить любым наиболее подходящим с учетом доступных
ресурсов методом.

Решением СЛУ (\ref{formula:SnLE-Newton-sys}) будет вектор
\(\Delta X^{(p)} = (\MFArgs)\). После этого, приближенное решение
исходной задачи находится по формуле
\(X^{(p+1)} = X^{(p)} + \Delta X^{(p)}\).

Заданная точность достигается выполнением условия
(\ref{formula:precision_iter_sle}). Стоит учитывать, что данный метод
так же, как и предыдущий, требователен к точности начального приближения.
\subsubsection{Реализации метода в библиотеках numpy, scipy}
Реализаций данного алгоритма в библиотеках numpy, scipy не найдено,
кроме того, для его работы требуется нахождение частных производных.

Для поиска частной производной в scipy есть функция \textbf{derivative}
в модуле \textbf{scipy.misc}, но она отмечена устаревшей и будет
удалена в версии 1.12.0, поэтому реализация данного метода с помощью
применения функций общей направленности библиотеки рассматриваться
не будет.

\section{Аппроксимация функций}
Иногда значения некоторой функциональной зависимости
\(y= \widetilde{y}(x) \) известны для отдельных пар значений
\((x_i,y_i)\).

Задача восстановления аналитической функции \(\widetilde{y}\)
по отдельным парам значений называется аппроксимацией.
Для получения ее однозначного решения необходимо задать общий вид
функции, включающей коэффициенты, и затем эти коэффициенты определить
с помощью подходящего метода.

Для определения коэффициентов существует два основных подхода ---
интерполяция (когда \(\widetilde{y}(x_i) = y_i\)), и сглаживание, когда
требуется лишь минимизировать отклонение от известных значений.

Первые три метода решают поставленную задачу с помощью интерполяции,
последний --- с помощью сглаживания.
\subsection{Интерполяционный полином Лагранжа}
\subsubsection{Описание метода}
Для известных пар значений \((x_i,y_i), x_i \in [a;b]; i = 1,2,3,\dots,n\)
строится полином \(P(x)\), удовлетворяющий условию \(P(x_i) = y_i\).

Полином \(P(x)\) определяется следующей формулой:
\begin{equation*}
   P(x) = L_{n-1}(x) = \sum_{i=1}^{n}y_i\cdot l_i(x)
\end{equation*}
\(l_i(x)\) --- фундаментальные полиномы Лагранжа, которые удовлетворяют равенству (\ref{formula:ip-lagr-eq1}), и имеют следующий вид:
\begin{equation*}
   l_i(x) = \frac{(x-x_1)\dots (x-x_{i-1})(x-x_{i+1})\dots(x-x_n)}{(x_i-x_1)\dots(x_i-x_{i-1})(x_i-x_{i+1})\dots(x_i-x_n)}
\end{equation*}
\begin{equation}
   l_k(x_i) =  \left\{
   \begin{aligned}
       & 1, i = k    \\
       & 0, i \ne k. \\
   \end{aligned}
   \right.
   \label{formula:ip-lagr-eq1}
\end{equation}

Из формулы полинома видно, что его степень не превышает \(n-1\).
Данное свойство сохранится и для следующего метода.

\subsubsection{Реализации метода в библиотеках numpy, scipy}
Данный метод реализован в библиотеке scipy в виде функции
\textbf{lagrange} в модуле \textbf{scipy.interpolate}.
Одной из ее особенностей является низкая устойчивость --- авторы не
рекомендуют запускать алгоритм на более чем 20 парах \((x,y)\) входных
значений \cite{links:scipy_doc}.

Данная функция имеет следующие параметры:
\begin{enumerate}
   \item \(x\) --- array\_like

         \(x\) представляет координаты \(x\) набора точек данных.
   \item \(w\) --- array\_like

         \(w\) представляет координаты \(y\) набора точек данных,
         т. е. \(f(x)\).
\end{enumerate}

Данная функция возвращает полином (тип --- \verb|poly1d| из
библиотеки \textbf{numpy} \cite{links:numpy_doc}).

\subsection{Интерполяционный полином Ньютона}
\subsubsection{Описание метода}
Интерполяционный полином Ньютона имеет вид
\begin{align*}
   N_{n-1}(x) = \Delta^{0}(x_{1}) + \Delta^{1}(x_{1},x_{2})(x-x_{1}) + \Delta^{2}(x_{1},x_{2},x_{3})(x-x_{1})(x-x_{2}) + \\
   + \dots + \Delta^{n-1}(x_{1},x_{2},\dots,x_{n-1})(x-x_{1})(x-x_{2})\dots(x-x_{n-1})
\end{align*}
где \(\Delta^{0}(x_{1})\dots\Delta^{n-1}(x_{1},x_{2},\dots,x_{n-1})\) ---
конечные разницы порядка \(0,\dots, n-1\), которые вычисляются следующим
образом:
\begin{nospaceflalign*}
   & \Delta^{0}(x_i) = y_i                                                               & \\
   & \Delta^{1}(x_i,x_k) = \frac{\Delta^{0}(x_i) - \Delta^{0}(x_k)}{x_i-x_k}             & \\
   & \Delta^{2}(x_i,x_j,x_k) = \frac{\Delta^{1}(x_i,x_j) - \Delta^{1}(x_j,x_k)}{x_i-x_k}
\end{nospaceflalign*}
и т.д.
\subsubsection{Реализации метода в библиотеках numpy, scipy}
Реализаций данного метода в в библиотеках numpy, scipy не найдено.

\subsection{Сплайн-интерполяция}
Перед рассмотрением метода необходимо ввести понятие интерполяционного
сплайна. Это кусочно-заданная функция, каждый фрагмент которой является
непрерывной функцией, и на концах его значение совпадает с известными
значениями функции.

Иногда выбирают фрагменты сплайна такого вида, что непрерывными будут и их
производные, вплоть до нужного порядка (например, если в роли сплайнов
выбрать полиномы 3 степени, то они будут соответствовать условиями
непрерывности вместе с 1 и 2 производными). В зависимости от цели,
на сплайны могут быть наложены и иные ограничения. Также, для упрощения
вычислений, чаще всего все фрагменты сплайна --- функции из одного
класса, например, многочленов фиксированной степени.
\subsubsection{Описание метода}
Между парами значений \((x_{l(j-1)},y_{l(j-1)}),(x_{l(j)},y_{l(j)})\),
где \(l \in \textbf{N}, 0 \le l_j \le n\) строятся
фрагменты сплайна \(f_1(x),\dots,f_n(x)\), соответствущие условиям
непрерывности. Комбинация таких функций образует исходную:
\begin{equation*}
   \widetilde{y}(x) = \left\{
   \begin{aligned}
       & f_1(x),\; x_0 \le x < x_l      \\
       & f_2(x),\; x_{l} \le x < x_{2l} \\
       & \ \dots \qquad \quad \dots     \\
       & f_n(x),\; x_{jl} \le x < x_n   \\
   \end{aligned}
   \right.
\end{equation*}
При этом
\begin{equation*}
   \begin{aligned}
       & f_1(x_0) = y_0, f_1(x_l) = y_l       \\
       & f_2(x_l) = y_l, f_2(x_{2l}) = y_{2l} \\
       & \ \dots \qquad \quad \dots           \\
       & f_n(x_{jl}) = y_{jl},f_n(x_n) = y_n  \\
   \end{aligned}
\end{equation*}
Для полиномов 3 степени чаще всего \(l=3\), и в таком случае по \(l\)
точкам на каждом фрагменте строиться полином (например, интерполяционный
полином Ньютона). В общем случае \(l\) определяется исходя из условий
и ограничений поставленной перед исследователем задачи.
\subsubsection{Реализации метода в библиотеках numpy, scipy}
В библиотеке scipy существует множество реализаций данного метода,
которые отличаются используемым видом сплайна а также поддерживаемой
размерностью функции. В модуле
\textbf{scipy.interpolate} описаны следующие объекты
\cite{links:scipy_doc}:
\begin{enumerate}
   \item класс \textbf{CubicSpline};
   \item класс \textbf{PchipInterpolator};
   \item класс \textbf{CubicHermiteSpline};
   \item класс \textbf{Akima1DInterpolator};
   \item классы \textbf{RectBivariateSpline},
         \textbf{LSQBivariateSpline}, используются для многомерной
         интерполяции;
   \item функция \textbf{interp1d}, которая может интерполировать таблично
         заданную функцию разными методами, в т.ч. сплайн-интерполяцией.
         Признана устаревшей и поэтому не будет рассмотрена;
\end{enumerate}


Классы \textbf{CubicSpline}, \textbf{PchipInterpolator},
\textbf{CubicHermiteSpline}, \textbf{Akima1DInterpolator}
используются для интерполяции функций одного аргумента, и их
применение одинаково: для создания сплайна необходимо создать
экземпляр класса, для получения значений сплайна необходимо
вызвать созданный экземпляр с необходимыми входными данными
(то есть вызвать метод \verb|__call__|).

Каждый из классов имеет разное количество параметров, которые
можно задать при его создании его экземпляра; метод получения
данных сплайна, в свою очередь, имеет одинаковое количество
параметров для всех классов, поэтому он имеет смысл описать его сейчас.

Метод \verb|__call__| вышеописанных классов имеет 3 параметра:
\begin{enumerate}
   \item \(x\) --- array\_like

         Точки, для которых будет вычислено значения сплайна.
   \item \(nu\) --- int, необязательный

         Порядок производной сплайна, который необходимо использовать при
         вычислении значений. Должен быть неотрицательным. По умолчанию
         --- 0 (то есть вычисляется исходный сплайн).
   \item \(extrapolate\) --- {bool, \verb|"periodic"|,
         \verb|None|}, необязательный

         Если bool, определяет, следует ли экстраполировать точки,
         которые выходят за пределы границ, установленные первым и
         последним интервалами или возвращать \verb|NaN|. Если
         \verb|"periodic"|, используется периодическая экстраполяция.
         Если \verb|None| (по умолчанию), используйте метод \textbf{extrapolate}
         этого класса.
\end{enumerate}
Данный метод возвращает массив значений сплайна, соответствующих значений
\(x\).

Конструктор экземпляра класса \textbf{CubicSpline} имеет следующие
параметры:
\begin{enumerate}
   \item \(x\) --- array\_like, одномерный

         Одномерный массив, содержащий значения независимой переменной.
         Значения должны быть действительными, конечными числами и
         находиться в строго возрастающем порядке.
   \item \(y\) --- array\_like

         Массив, содержащий значения зависимой переменной. Он может
         иметь произвольное количество измерений, но длина должна
         совпадать с длиной \(x\). Значения должны быть конечными.
   \item \(axis\) int, необязательный

         Ось, вдоль которой предполагается изменение y. Это означает,
         что для \(x[i]\) соответствующие значения равны
         \verb|np.take(y, i, axis=axis)|. По умолчанию --- 0.
   \item \(bc\_type\) --- string / tuple размера = 2, необязательный

         Тип граничного условия. Два дополнительных уравнения, заданные
         граничными условиями, необходимы для определения всех
         коэффициентов многочленов на каждом отрезке.

         Если \(bc\_type\) --- строка, то указанное в параметре условие
         будет применено к обоим концам сплайна. Доступные условия:
         \begin{itemize}
            \item \verb|"not-a-knot"| (по умолчанию): первый и второй
                  сегменты на конце кривой представляют собой один и тот
                  же полином. Это хороший вариант по умолчанию, когда нет
                  информации о граничных условиях.
            \item \verb|"periodic"|: Предполагается, что интерполируемые
                  функции являются периодическими с периодом
                  \(x[-1] \ \mbox {---}\ x[0]\).
                  Первое и последнее значение y должны быть идентичными:
                  \(y[0] == y[-1]\). Это граничное условие приведет к
                  тому, что \(y'[0] == y'[-1]\) и \(y''[0] == y''[-1]\).
            \item \verb|"clamped"|: Первая производная на концах кривых
                  равна нулю. Предполагая, что y одномерный,
                  \(bc\_type=((1, 0.0), (1, 0.0))\) --- то же самое
                  условие.
            \item \verb|"natural"|: Вторая производная на концах кривой
                  равна нулю. Предполагая, что y одномерный,
                  \(bc\_type=((2, 0.0), (2, 0.0))\) --- то же самое
                  условие.
         \end{itemize}

         Если \(bc\_type\) представляет собой кортеж из двух элементов,
         первое и второе значения будут применены в начале и конце
         кривой соответственно. Значения кортежа могут быть одной из
         ранее упомянутых строк (кроме \verb|"periodic"|) или кортежем
         \((order, deriv\_values)\), позволяющим указывать произвольные
         производные на концах кривой:

         \begin{itemize}
            \item \(order\): порядок производной --- 1 или 2.
            \item \(deriv\_value\): array\_like

                  Содержит значения производной, форма должна быть такой
                  же, как \(y\), за исключением измерения \(axis\).
                  Например, если \(y\) --- одномерный, то
                  \(deriv\_values\) должен быть скаляром. Если \(y\)
                  является трехмерным, имеет форму \((n0, n1, n2)\) и
                  \(axis=1\), то \(deriv\_values\) должно быть
                  двухмерным и иметь форму \((n0, n2)\).
         \end{itemize}
   \item \(extrapolate\) --- {bool, \verb|"periodic"|,
         \verb|None|}, необязательный

         Если bool, определяет, следует ли экстраполировать точки,
         которые выходят за пределы границ, установленные первым и
         последним интервалами или возвращать \verb|NaN|. Если
         \verb|"periodic"|, используется периодическая экстраполяция.
         Если \verb|None| (по умолчанию), параметр имеет значение
         \verb|"periodic"| если \(bc\_type\)=\verb|"periodic"| и
         значение \verb|True| в противном случае.
\end{enumerate}

Класс \textbf{PchipInterpolator} представляет из себя кубический
эрмитов сплайн, поэтому является дочерним классом
\textbf{CubicHermiteSpline}. Конструктор экземпляра класса имеет
следующие параметры:
\begin{enumerate}
   \item \(x\) --- ndarray, вида (xlen,1), то есть одномерный

         Одномерный массив монотонно возрастающих действительных
         значений. Не может включать повторяющиеся значения (иначе
         сплайн будет слишком конкретизированный)
   \item \(y\) --- ndarray, вида (...,xlen,...)

         N-мерный массив действительных значений. Длина \(y\) по оси
         интерполяции должна быть равна длине \(x\). Используйте
         параметр \(axis\), чтобы выбрать ось интерполяции.
   \item \(axis\) int, необязательный

         См. описание одноименного метода родительского класса.
   \item \(extrapolate\) --- bool, необязательный

         Определяет, следует ли экстраполировать точки,
         которые выходят за пределы границ, установленные первым и
         последним интервалами или возвращать \verb|NaN|.
         По умолчанию \verb|None|, значение см. в пункте
         \ref{list:CubicHermiteSpline-extrapolate} описания
         конструктора родительского класса.
\end{enumerate}

Конструктор экземпляра класса \textbf{CubicHermiteSpline} имеет следующие
параметры:
\begin{enumerate}
   \item \(x\) --- array\_like, вида (n,), то есть одномерный

         Одномерный массив, содержащий значения независимой переменной. Значения должны быть действительными, конечными и находиться в строго возрастающем порядке.
   \item \(y\) --- array\_like

         Массив, содержащий значения зависимой переменной.
         Он может иметь произвольное количество измерений, но длина по
         оси \(axis\) должна совпадать с длиной \(x\).
         Значения должны быть конечными.
   \item \(dydx\) --- array\_like

         Массив, содержащий производные зависимой переменной. Он может
         иметь произвольное количество измерений, но длина по оси
         \(axis\) должна совпадать с длиной \(x\). Значения должны быть
         конечными.
   \item \(axis\) --- int, необязательный

         Ось, вдоль которой предполагается изменение y. Это означает, что для \(x[i]\) соответствующие значения равны
         \verb|np.take(y, i, axis=axis)|.
         По умолчанию --- 0.
   \item \(extrapolate\) --- {bool, \verb|"periodic"|,
         \verb|None|}, необязательный

         Если bool, определяет, следует ли экстраполировать точки,
         которые выходят за пределы границ, установленные первым и
         последним интервалами или возвращать \verb|NaN|. Если
         \verb|"periodic"|, используется периодическая экстраполяция.
         Если \verb|None| (по умолчанию), используйте метод \textbf{extrapolate}
         этого класса.
         \label{list:CubicHermiteSpline-extrapolate}
\end{enumerate}

Конструктор экземпляра класса \textbf{Akima1DInterpolator} имеет
следующие параметры:
\begin{enumerate}
   \item \(x\) --- ndarray, вида (n, 1)

         Одномерный массив монотонно возрастающих действительных
         значений.
   \item \(y\) --- ndarray, вида (..., n, ...)

         N-мерный массив реальных значений. Длина \(y\) по оси
         интерполяции должна быть равна длине \(x\). Используйте
         параметр \(axis\), чтобы выбрать ось интерполяции.
   \item \(axis\) --- int, необязательный

         Ось в массиве \(y\), соответствующая значениям координаты
         \(x\). По умолчанию --- 0.
\end{enumerate}

\subsection{Сглаживание. Метод наименьших квадратов}
С помощью метода наименьших квадратов (далее, МНК) осуществляется
интерполяция функции с помощью сглаживания.
\subsubsection{Описание метода}
Задача метода --- минимизировать отклонение от исходных данных:
\begin{equation}
   S = \sum_{i=1}^{n}(\widetilde{y}(x_i,a)-y_i)^2 \rightarrow \min
   \label{formula:mnk}
\end{equation}
где \(a\) --- неизвестные параметры; \(n\) --- количество известных
пар значений (\(x,y\)); \(x_i,y_i\) --- значения \(i\) пары.

В качестве \(\widetilde{y}(x)\) необходимо выбрать функцию вида,
который лучше всего подходит под исходные данные. Этот процесс
называют выбором зависимости.

Примеры зависимостей (\(a,b,c,d\) --- неизвестные параметры):
\begin{alignat*}{2}
    & ax + b        & \  & \text{--- линейная};         \\
    & ax^2 + bx + c & \  & \text{--- полиномиальная};   \\
    & a \ln(x) + b  & \  & \text{--- логарифмическая};  \\
    & be^{ax}       & \  & \text{--- экспоненциальная};
\end{alignat*}

В результате подстановки выбранной зависимости и известных данных
в (\ref{formula:mnk}) задача минимизации сводится к нахождению
экстремумов функции \(S(a)\), и выбору подходящего из них.
\subsubsection{Реализации метода в библиотеках numpy, scipy}
В модуле \textbf{scipy.optimize} библиотеки scipy существует
несколько реализаций данного метода, каждая из которых имеет свою
специализацию.
\begin{enumerate}
   \item \textbf{curve\_fit} --- функция, предназначенная для поиска
         значения параметра \(a\) произвольной зависимости \(\widetilde
         {y}(x_i,a)\), при котором выполняется условие
         (\ref{formula:mnk}).
   \item \textbf{least\_squares} --- функция, которая отличается от
         \textbf{curve\_fit} возможностью задать собственный критерий
         близости, область допустимых решений; позволяет задать больше
         параметров точности.
   \item \textbf{leastsq} --- функция, которая является устаревшей
         версией \textbf{least\_squares}, поэтому ее подробного
         рассмотрения не будет.
   \item \textbf{lsq\_linear} --- функция, предназначенная для
         минимизации выражения \(||ax-b||^2\), где \(a,b\) --- матрицы
         параметров; функция \textbf{nnls} отличается от нее
         ограничением на положительность искомых значений.

         Подходят только для нахождения параметров многочленов,
         не буду рассматриваться в данной работе по причине
         сложной подготовки входных данных для решения задачи
         метода, а также их узкой специализации.
\end{enumerate}

Функция \textbf{curve\_fit} имеет следующий набор параметров:
\begin{enumerate}
   \item \(f\) --- callable

         Модельная функция \(f(x,\dots)\). Она должен принимать независимую переменную в качестве первого аргумента, а параметры --- в качестве оставшихся.
   \item \(xdata\) --- array\_like

         Значения независимой переменной. Обычно это должна быть последовательность длины \(M\) или массив \((k,M)\)-образной формы для функций с \(k\) предикторами, и каждый элемент должен быть конвертируемым в \(float\), если это объект, подобный массиву.
   \item \(ydata\) --- array\_like

         Значения зависимой переменной, массив длиной M.
         При найденных параметрах \(f\) совпадает с значением
         \(f(xdata, \dots)\).
   \item \(p0\) --- array\_like, необязательный

         Начальное значение параметров \(f\) (длины N). Если
         \verb|None| (значение по умолчанию), то все начальные
         значения будут равны \verb|1| (если количество параметров
         функции можно определить с помощью интроспекции, в противном
         случае генерируется исключение \verb|ValueError|).
   \item \(sigma\) --- \verb|None| /  sequence длины M / массив MxM, необязательный

         Задает неопределенность в \(ydata\). Если мы определим
         остатки как \(r = ydata - f(xdata, *popt)\), то интерпретация
         значения \(sigma\) будет зависеть от ее количества измерений:
         \begin{itemize}
            \item одномерная \(sigma\) должна содержать значения
                  стандартных отклонений ошибок в \(ydata\). Используется
                  при вычислении параметра
                  \verb|chisq = sum((r/sigma)**2)|;
            \item двухмерная \(sigma\) должна содержать ковариационную
                  матрицу ошибок в \(ydata\). Используется
                  при вычислении параметра
                  \verb|chisq = r.T @ inv(sigma) @ r|.
         \end{itemize}
   \item \(absolute\_sigma\) --- bool, необязательный

         Если значение \verb|True|, \(sigma\) используется в абсолютном смысле, и предполагаемая ковариация параметра \(pcov\) отражает эти абсолютные значения.

         Если значение \verb|False| (по умолчанию), имеют значение
         только относительные величины значений \(sigma\).
         Ковариационная матрица возвращаемого параметра \(pcov\)
         основана на масштабировании \(sigma\) постоянным
         коэффициентом. Эта константа устанавливается требованием, чтобы
         приведенный выше \(chisq\) для параметров
         \(popt\) при использовании масштабированной \(sigma\) был равен
         единице. Другими словами, \(sigma\) масштабируется так, чтобы
         соответствовать выборочной дисперсии остатков после подгонки.
         Математически,
         \verb|pcov(absolute_sigma=False) = pcov(absolute_sigma=True)|
         \verb| * chisq(popt)/(M-N)|.
   \item \(check\_finite\) --- bool, необязательный

         Если \verb|True|, то входные массивы проверяются на содержание
         значений \verb|NaN| или \verb|Inf|, и если они содержат, то
         генерируется исключение \verb|ValueError|.

         Значение \verb|False| может привести к выдаче бессмысленных
         результатов, если входные массивы содержат \verb|NaN|.

         По умолчанию установлено значение \verb|True|, если
         \(nan\_policy\) не указано явно, и значение \verb|False| в
         противном случае.
   \item \(bounds\) --- tuple из двух array\_like / \verb|Bounds|,
         необязательный

         Нижняя и верхняя границы параметров \(f\). По умолчанию нет ограничений. Есть два способа указать границы:
         \begin{itemize}
            \item через экземпляр класса \verb|Bounds|;

            \item \verb|tuple| из двух array\_like, где каждый элемент
                  кортежа должен быть либо массивом с длиной, равной
                  количеству параметров, либо скаляром (в этом случае
                  граница считается одинаковой для всех параметров).
                  Используйте \verb|np.inf| с соответствующим знаком,
                  чтобы отключить ограничения для всех или некоторых
                  параметров.
         \end{itemize}
   \item \(method\) --- \{\verb|"lm"|, \verb|"trf"|, \verb|"dogbox"|\}, необязательный

         Метод, используемый для вычисления параметров. См. описание
         аналогичного параметра функции \textbf{least\_squares} для
         получения более подробной информации. По умолчанию используется
         \verb|"lm"| для задач без ограничений и \verb|"trf"|, если
         указано значение \(bounds\).
   \item \(jac\) --- callable / string / \verb|None|, необязательный

         Функция с сигнатурой \(jac(x, ...)\), которая вычисляет матрицу
         Якоби модельной функции относительно параметров как плотную
         array\_like структуру. Он будет масштабироваться в
         соответствии с предоставленным параметром \(sigma\). Если
         \verb|False| (по умолчанию), якобиан будет оцениваться
         численно. Строковые ключевые слова для методов \verb|"trf"| и
         \verb|"dogbox"| можно использовать для выбора
         конечно-разностной схемы, см. описание
         \textbf{least\_squares}.
   \item \(full\_output\) --- boolean, необязательный

         Если значение = \verb|True|, эта функция возвращает дополнительную информацию: \(infodict\), \(mesg\) и \(ier\). По умолчанию
         \verb|False|.
   \item \(nan\_policy\) --- {\verb|"raise"|, \verb|"omit"|, \verb|None|}, необязательный

         Определяет, как действовать, если входные данные содержат \verb|NaN|. Доступны следующие параметры (по умолчанию --- \verb|None|):
         \begin{itemize}
            \item \verb|"raise"|: выдает ошибку

            \item  \verb|"omit"|: выполняет вычисления, игнорируя
                  значения \verb|NaN|.

            \item  \verb|None|: никакая специальная обработка \verb|NaN|
                  не выполняется (кроме того, что делается с помощью
                  \(check\_finite\)); поведение при наличии
                  \verb|NaN| зависит от реализации и может измениться.

         \end{itemize}
         Обратите внимание: если значение указано явно (не \verb|None|),
         для \(check\_finite\) будет установлено значение \verb|False|.
   \item \(**kwargs\)

         Именованные параметры, которые передаются в \textbf{leastsq}
         для \(method\)='lm' или \textbf{least\_squares} в противном
         случае.
\end{enumerate}
Данная функция возвращает результат в виде
\((popt,pcov,infodict,mesg,ier)\), где
\begin{enumerate}
   \item \(popt\) --- вычисленные параметры
         \(f\), при которых выполняется условие
         \(f(xdata, *popt) - ydata \rightarrow \min\).
   \item \(pcov\) --- 2-D array

         Вычисленная приблизительная ковариация \(popt\).
   \item  \(infodict\) --- dict

         По ключу \(nfev\)
         хранится значение количества вызовов функции. Методы \("trf"\)
         и \("dogbox"\) не учитывают вызовы функций для численной
         аппроксимации якобиана, в отличие от метода \("lm"\).
         По ключам \(fvec\), \(fjac\), \(ipvt\), \(qtf\) можно получить
         дополнительную информацию о работе алгоритма.
   \item \(mesg\) --- str

         Строковое сообщение с информацией о решении.
   \item \(ier\) --- int

         Целочисленный флаг. Если он равен 1, 2, 3 или 4, решение
         найдено. В противном случае решение не было найдено. \(mesg\)
         предоставляет дополнительную информацию о решении.
\end{enumerate}
Разработчики отмечают, что входные данные \(xdata\), \(ydata\) и
выходные данные \(f\) должны иметь формат \(float64\), иначе данная
функция может вернуть неверные результаты \cite{links:scipy_doc}.

Функция \textbf{least\_squares} предназначена для решения следующей
задачи минимизации \(F(x)\), где
\verb|F(x) = 0.5 * sum(rho(f_i(x)**2), i = 0, ..., m - 1)|
при ограничениях \(lb \leq x \leq ub\). То есть, \(f\) должна принимать
на вход значение параметров \(a\) и возвращать значение \(f(x,a) - y\).

Данная функция имеет следующий набор параметров:
\begin{enumerate}
   \item \(fun\) --- callable
         Функция, которая вычисляет вектор остатков, имеет сигнатуру
         \(fun(x, *args, **kwargs)\), т.е. минимизация продолжается
         относительно ее первого аргумента. Аргумент \(x\),
         передаваемый этой функции, представляет собой N-мерный массив размерности \((N,K)\) (никогда не скаляр, даже для \(N=1\)). \(fun\) должна вернуть одномерный массив типа array \((m,)\) или скаляр. Если \(x \in \textbf{C}\) или функция \(fun\) возвращает комплексные числа, ее необходимо обернуть действительнозначной функцией от действительнозначных аргументов.
   \item \(x_0\) --- array\_like размерности (N,K) / float

         Начальное значение \(x\). Если параметр имеет тип float, он
         будет рассматриваться как одномерный массив с одним элементом.
   \item \(jac\) --- \{\verb|"2-point"|, \verb|"3-point"|, \verb|"cs"|, callable\}, необязательный

         Метод вычисления матрицы Якоби (матрица размером \(m \times  n\), где элемент \((i, j)\) является частной производной \(f[i]\) по \(x[j]\)). Ключевые слова выбирают схему конечных разностей для числовой оценки. Схема \verb|"3-point"| более точная, но требует в два раза больше операций, чем \verb|"2-point"| (по умолчанию). Схема \verb|"cs"| использует комплексные шаги и, хотя потенциально является наиболее точной, она применима только тогда, когда \(fun\) правильно обрабатывает комплекснозначные входные данные и может быть аналитически продолжена на комплексную плоскость. \(method\) \verb|"lm"| всегда использует схему \verb|"2-point"|. Если параметр имеет тип \verb|callable|, он используется как \verb|jac(x, *args, **kwargs)| и должен возвращать хорошее приближение (или точное значение) для якобиана в виде array\_like (к выводу \(jac\) будет применен \verb|np.atleast_2d|), разреженной матрицы (предпочтительна \verb|csr_matrix| для производительности) или \verb|scipy.sparse.linalg.LinearOperator|.
   \item \(bounds\) --- tuple из 2 array\_like / \verb|Bounds|, необязательный

         Есть два способа указать границы:
         \begin{itemize}
            \item через экземпляр класса \verb|Bounds|;

            \item \verb|tuple| из двух array\_like, где каждый элемент
                  кортежа должен быть либо массивом с размерностью
                  \(x0\), либо скаляром (в этом случае
                  граница считается одинаковой для всех параметров).
                  Используйте \verb|np.inf| с соответствующим знаком,
                  чтобы отключить ограничения для всех или некоторых
                  параметров.
         \end{itemize}
   \item \(method\) --- \{\verb|"lm"|, \verb|"trf"|, \verb|"dogbox"|\}, необязательный

         Допустимое отклонение на завершение по норме градиента.
         По умолчанию --- 1e-8. Точное состояние зависит от используемого
         значения \(method\):
         \begin{itemize}
            \item Для \verb|"trf"|: \(norm(g\_scaled, ord=np.inf) < gtol\), где \(g\_scaled\) --- значение градиента, масштабированное для учета наличия границ.
            \item Для \verb|"dogbox"|: \(norm(g_free, ord=np.inf) < gtol\), где \(g_free\) --- градиент по отношению к переменным, которые не находятся в оптимальном состоянии на границе.
            \item Для \verb|"lm"|: максимальное абсолютное значение косинуса углов между столбцами якобиана и вектором невязки меньше \(gtol\), или вектор невязки равен нулю.
         \end{itemize}

         Если \verb|None| и \(method\) не \verb|"lm"|, завершение по этому условию отключено. Если \(method\) --- \verb|"lm"|, этот допуск должен быть выше, чем машинный эпсилон.
   \item \(x\_scale\) --- array\_like / \verb|"jac"|

         Характеристический масштаб каждой переменной.
         Установка \(x\_scale\) эквивалентна переформулировке проблемы
         в масштабированных переменных \(xs = x/x\_scale\).
         Альтернативная точка зрения состоит в том, что размер
         доверительной области по j-му измерению пропорционален
         \(x\_scale[j]\). Улучшения сходимости можно достичь, установив
         \(x\_scale\) таким образом, чтобы шаг заданного размера по
         любой из масштабируемых переменных имел аналогичный эффект на
         функцию стоимости. Если установлено значение \verb|"jac"|,
         масштаб итеративно обновляется с использованием обратных норм
         столбцов матрицы Якоби.
   \item \(loss\) --- str / callable, необязательный

         Определяет функцию потерь. Допускаются следующие значения ключевых слов:
         \begin{itemize}
            \item \verb|"linear"| (по умолчанию) : \verb|rho(z) = z|. Дает стандартную задачу наименьших квадратов.
            \item \verb|"soft_l1"| : \verb|rho(z) = 2 * ((1 + z)**0.5 - 1)|. Гладкая аппроксимация потерь l1 (абсолютное значение). Обычно хороший выбор для надежного метода наименьших квадратов.
            \item \verb|"huber"| : \verb|rho(z) = z if z <= 1 else 2*z**0.5 - 1|. Работает аналогично \verb|"soft_l1"|.
            \item \verb|"cauchy"| : \verb|rho(z) = ln(1 + z)|. Сильно ослабляет влияние выбросов, но может вызвать трудности в процессе оптимизации.
            \item \verb|"arctan"| : \verb|rho(z) = arctan(z)|. Ограничивает максимальный убыток по одному остатку, имеет свойства, аналогичные \verb|"cauchy"|.
         \end{itemize}

         Если параметр имеет тип \verb|callable|, он должен принимать
         одномерный \verb|ndarray| \verb|z=f**2| и возвращать
         array\_like с формой \((3, m)\), где строка 0 содержит значения функции, строка 1 содержит первые производные, а строка 2 содержит вторые производные. \(method\) \verb|lm| поддерживает только \verb|"linear"| \(loss\).

   \item \(f\_scale\) --- float, необязательный

         Значение границы между ошибками данных и выбросами в данных,
         по умолчанию - 1,0. Значение функции потерь оценивается
         следующим образом: \verb|rho_(f**2) = C**2 * rho(f**2 / C**2)|,
         где C --- это \(f\_scale\), а \(rho\) определяется параметром
         \(loss\). Этот параметр не имеет никакого эффекта при
         \(loss\)=\verb|"linear"|, но для других значений потерь он имеет решающее
         значение.
   \item \(max\_nfev\) --- \verb|None| / int, необязательный

         Максимальное количество выполнений функции перед завершением работы. Если \verb|None| (по умолчанию), значение выбирается автоматически:
         \begin{itemize}
            \item Для  \verb|"trf"| и  \verb|"dogbox"|: \(100 \cdot n\).

            \item Для \verb|"lm"|: \(100 \cdot n\), если \(jac\) является вызываемым, и \(100 \cdot n \cdot (n + 1)\) в противном случае (поскольку \verb|"lm"| учитывает вызовы функций в оценке якобиана).
         \end{itemize}
   \item \(diff\_step\) --- \verb|None| / array\_like, необязательный

         Определяет относительный размер шага для аппроксимации якобиана конечной разностью. Фактический шаг вычисляется как \verb|x * diff_step|. Если \verb|None| (по умолчанию), то \(diff_step\) принимается за обычную "оптимальную" степень машинного эпсилона для используемой схемы конечных разностей.
   \item \(tr\_solver\) --- {\verb|None|, \verb|"exact"|, \verb|"lsmr"|}, необязательный

         Метод решения подзадач доверительной области, применим только для методов \verb|"trf"| и \verb|"dogbox"|.
         \begin{itemize}
            \item \verb|"exact"| подходит для не очень больших задач с плотными матрицами Якоби. Вычислительная сложность на итерацию сравнима с сингулярным разложением матрицы якобиана.
         \end{itemize}
   \item \(tr\_options\) --- dict, необязательный

         Параметры ключевых слов передаются в решатель доверительной
         области.
         \begin{itemize}
            \item \verb|tr_solver="exact"|: \(tr\_options\) игнорируются.
            \item \verb|tr_solver="lsmr"|: опции для \textbf{scipy.sparse.linalg.lsmr}.
                  Кроме того, \(method\)=\verb|"trf"| поддерживает опцию \verb|'regularize'| (bool, по умолчанию --- \verb|True|), которая добавляет член регуляризации к нормальному уравнению, что улучшает сходимость, если якобиан имеет недостаточный ранг.
         \end{itemize}
   \item \(jac\_sparsity\) --- {\verb|None|, array\_like, разреженная матрица}, необязательный

         Определяет структуру разреженности матрицы якобиана для
         конечно-разностной оценки, ее форма должна быть (m, n).
         Если якобиан имеет лишь несколько ненулевых элементов в
         каждой строке, обеспечение разреженной структуры значительно
         ускорит вычисления. Нулевая запись означает, что
         соответствующий элемент якобиана тождественно равен нулю.
         Если предусмотрено, принудительно используется решатель
         доверительной области \verb|"lsmr"|. Если \verb|None|
         (по умолчанию), будет использоваться плотная разность.
         Не имеет эффекта при \(method\)=\verb|"lm"|.
   \item \(verbose\) --- {0, 1, 2}, необязательный

         Уровень детализации алгоритма:
         \begin{itemize}
            \item 0 (по умолчанию): бесшумная работа функции.
            \item 1: отображение отчета о завершении.
            \item 2: отображение хода выполнения во время итераций
                  (не поддерживается методом \verb|"lm"|).
         \end{itemize}
   \item \(args, kwargs\) --- tuple / dict, необязательный

         Дополнительные аргументы передаваемые в \(fun\) и \(jac\).
         Оба пусты по умолчанию. Сигнатура вызова для \(fun\) ---
         \(fun(x, *args, **kwargs)\), та же самая для \(jac\).
\end{enumerate}
Функция возвращает результат в виде экземпляра класса \verb|OptimizeResult|, у которого определены поля \(x, cost, fun, jac, grad, optimality, active\_mask, nfev, njev,status, \) \(message, success\).


\section{Численное интегрирование}
Задача численного интегрирования состоит в вычислении определенного
интеграла:
\begin{equation}
   I = \int_{a}^{b} f(x) \, dx
   \label{formula:nm-integral1}
\end{equation}

Для решения данной задачи выберем на отрезке \([a;b]\) \(n\)
различных узлов:
\(a = x_0 < x_1 < x_2 < ... < x_{n-1} < x_n = b\), затем по выбранным
узлам интерполируем \(f(x)\). Если в роли интерполяционной функции
выбрать полином \(P_m(x)\), то интеграл (\ref{formula:nm-integral1})
можно вычислить по формуле
\begin{equation*}
   I \approx \int_{a}^{b} P_m(x) \, dx \quad.
\end{equation*}

Данную формулу также называют квадратурной формулой интерполяционного
типа.

К наиболее распространенным методам относят метод прямоугольников,
метод трапеций и метод парабол (Симпсона). Первый метод использует
использует полиномы 0 степени, второй --- 1 степени, третий ---
2 степени.

В данной работе будет рассмотрены методы парабол и трапеций.

Метод парабол будет точнее других в большинстве случаев
(он позволяет находить точное решение для любых
f\(x\), если \(f^{(4)}(x) = 0, x \in [a;b] \), в соответствии с
формулой Джузеппе Пеано). Из его недостатков можно отметить низкую
точность на пилообразных функциях (т.е. значение которой резко
возрастают на отрезках малой длины).

При этом, если количество точек, по которым строится \(P_m(x)\), четно, то метод трапеций может оказаться удобнее, тем самым прекрасно дополняя
метод парабол.
\subsection{Метод трапеций}
\subsubsection{Описание метода}
На отрезке \([a;b]\) выбирается \(n\) узлов, на отрезках \([x_i;x_{i+1}]\) строятся интерполяционные многочлены 1 степени \(P_1(x)\).
Если используется метод Лагранжа, то \(P_1(x)\) определяется так:
\begin{equation*}
   P_1(x) = f(x_i) \frac{x-x_{i+1}}{x_i-x_{i+1}} + f(x_{i+1}) \frac{x-x_i}{x_{i+1}-x_i}
\end{equation*}
При этом
\begin{equation*}
   \int_{x_i}^{x_{i+1}} P_1(x) \, dx = \frac{1}{2} \left(f(x_i)+f(x_{i+1})\right)(x_{i+1} - x_i)
\end{equation*}
В случае равноотстоящих узлов
\(x_0, x_1 = x_0 + h, ..., x_n = x_0 + nh\) значение интеграла будет таким:
\begin{equation*}
   I \approx \frac{h}{2} \sum_{i=0}^{n-1} (f(x_i) + f(x_{i+1}))
\end{equation*}
\subsubsection{Реализации метода в библиотеках numpy, scipy}
В библиотеке scipy найдена функция \textbf{trapezoid} в модуле
\textbf{scipy.integrate}. Она имеет следующие параметры:
\begin{enumerate}
   \item \(y\) --- array\_like

         Массив входных данных по которым будет вычисляться интеграл
   \item \(x\) --- array\_like, необязательный

         Массив значений \(x\), соответствующих \(y\). Если \verb|None|
         (по умолчанию), то в роли \(x\) будет создан массив равноотстоящих
         значений (\(h = dx\))
   \item \(dx\) --- scalar, необязательный

         Расстояние между соседними значениями \(x\). Значение используется
         когда \(x\)=\verb|None|. По умолчанию --- 1.
   \item \(axis\) --- int, необязательный

         Ось, относительно которой производить вычисление интеграла. По
         умолчанию --- \(-1\).
\end{enumerate}

Данная функция возвращает число, если \(y\) одномерный и массив
размерности \((n-1)\) для \(n\)-мерного \(y\).
\subsection{Метод парабол (Симпсона)}
\subsubsection{Описание метода}
На отрезке \([a;b]\) выбирается \(2n+1\) узлов, \(n\) троек точек
\(x_{i-1},x_i,x_{i+1} (i = 1,3,5,...\, ,2n-1)\) с соответствующими
отрезками \([x_{i-1},x_{i+1}]\).

На каждом из отрезков интерполируем \(f(x)\) полиномом
2 степени \(P_2(x)\), например, через формулу Лагранжа:
\begin{equation*}
   \begin{split}
      P_2(x) = \ & f(x_{i-1})
      \frac{(x-x_i)(x-x_{i+1})}{(x_{i-1}-x_i)(x_{i-1}-x_{i+1})} + f(x_i)
      \frac{(x-x_{i-1})(x-x_{i+1})}{(x_i-x_{i-1})(x_i-x_{i+1})} + \\
      & + f(x_{i+1}) \frac{(x-x_{i-1})(x-x_i)}{(x_{i+1})(x_i+1-x_i)}
   \end{split}
\end{equation*}
При этом
\begin{equation*}
   \int_{x_{i-1}}^{x_{i+1}} P_2(x)\, dx = \frac{h}{3}
   \left( f(x_{i-1})+4f(x_i)+f(x_{i+1}) \right), \quad h = \frac{b-a}{N}
\end{equation*}
В результате, значение интеграла \(I\) будет вычисляться по формуле
\begin{equation*}
   I \approx  \frac{h}{3} \sum_{i=0}^{n-1}(f(x_{2i})+4f(x_{2i+1})+f(x_{2i+2}))
\end{equation*}
\subsubsection{Реализации метода в библиотеках numpy, scipy}
В библиотеке scipy есть реализация данного метода в виде
функции \textbf{simpson} модуля \textbf{scipy.integrate}.

Данная функция имеет следующие параметры:
\begin{enumerate}
   \item \(y\) --- array\_like

         Массив, который необходимо интегрировать.
   \item \(x\) --- array\_like, необязательный

         Если задано, точки, в которых производится выборка \(y\).
         \verb|None| по умолчанию.
   \item \(dx\) --- float, необязательный

         Расстояние между точками интегрирования по оси \(Ox\).
         Используется только тогда, когда \(x\)=\verb|None|.
         По умолчанию --- 1.
   \item \(axis\) --- int, необязательный

         Ось, по которой следует интегрировать. По умолчанию
         --- последняя ось (-1).
   \item \(even\) --- необязательный

         Устаревший параметр, будет удален в scipy \(1.13.0\).
         Отвечает за вычисление значения если количество точек четно.
         По умолчанию используется метод Симпсона для первых \(n-2\)
         отрезков с добавлением трехточечного параболического сегмента
         для последнего интервала с использованием уравнений, изложенных
         Картрайтом \cite{journal:cartwright}.
\end{enumerate}
Данная функция возвращает значение типа \(flot\) --- приблизительное
значение интеграла.

Стоит учесть, что если ось, по которой необходимо интегрировать, имеет
только две точки, то интегрирование происходит с помощью метода трапеций.

\section{Численное решение задачи Коши обыкновенных дифференциальных уравнений}
Обыкновенное дифференциальное уравнение (далее, ОДУ) --- это уравнение
вида \(F(x,y,y',y'',...,y^{(n)})\), где \(n\) --- порядок уравнения.
Решение ОДУ --- функция \(y=y(x)\), при подстановке которой в исходное
ОДУ получается верное тождество.

Так как для нахождения решения необходимо провести \(n\) интегрирований,
то общее решение ОДУ --- \(y=\varphi(x,C_1,C_2,...,C_n)\).

Частное решение, при котором будут найдены конкретные значения
\(C_1,...,C_n\), называется задачей Коши и для получения ее решения
необходимо, чтобы были заданы начальные условия:
\begin{equation*}
   y(x_0) = y_0, \quad y'(x_0) = y'_0, \ ...\ , y^{(n-1)}(x_0) = y^{(n-1)}_0
\end{equation*}

В данном разделе будут описаны разностные методы решения задачи Коши.
Процесс решения задачи с помощью таких методов состоит из следующих
этапов:
\begin{enumerate}
   \item Выборка узлов сетки --- дискретного множества точек из исходной
         области изменения аргумента
   \item Аппроксимация производных в узлах сетки конечно-разностными
         аналогами.
   \item Аппроксимация ОДУ системой разностных уравнений.
   \item Решение системы разностных уравнений.
\end{enumerate}
\subsection{Метод Эйлера}
Данный метод, как и его модифицированная версия, применяется для
ОДУ 1 порядка. Для решения задач Коши ОДУ порядка \(n\) необходимо
исходное уравнение свести к системе ОДУ 1 порядка с помощью замены переменных.
\subsubsection{Описание метода}
Дано ОДУ 1 порядка с начальными условиями:
\begin{equation*}
   y' = f(x,y), \qquad y(x_0) = y_0, x \in [a;b]
\end{equation*}

На отрезке \([a;b]\) выберем \(n\) точек:
\begin{equation}
   a = x_0 < x_1 < ... < x_n = b, \qquad x_{i+1} - x_i = h, i = 0,1,2,...,n-1
   \label{formula:euler}
\end{equation}

Затем получаем значение производной:
\begin{equation*}
   y'(x_i) \approx \frac{\Delta y}{\Delta x} = \frac{y_{i+1} - y_i}{h}
\end{equation*}
Исходя из (\ref{formula:euler}) получаем формулу Эйлера:
\begin{equation*}
   y_{i+1} = y_i + hf(x_i, y_i), \quad i = 0,1, ...,n-1
\end{equation*}
С помощью нее последовательно получаем значения \(y_1, y_2, ..., y_n\),
по которым можно провести интерполяцию, чтобы получить
\(y = \widetilde{y}(x)\).

Погрешность метода на каждой итерации --- \(O(h^2)\).
\subsubsection{Реализации метода в библиотеках numpy, scipy}
Реализаций данного метода в библиотеках numpy, scipy не найдено.

\subsection{Модифицированный метод Эйлера}
\subsubsection{Описание метода}
Итерационную формулу предыдущего метода разложим в ряд Тейлора:
\begin{equation}
   y_{i+1} = y_i + y'_i h + \frac{1}{2}y''h^2 + O(h^3)
   \label{formula:mod_euler}
\end{equation}
Приближенное значение \(y''\) вычисляется аналогично \(y'\):
\begin{equation*}
   y''_i = \frac{y'_{i+1} - y'_i}{h}
\end{equation*}
Подставляем данное выражение в (\ref{formula:mod_euler}), пренебрегая \(O(h^3)\):
\begin{equation*}
   y_{i+1} = y_i + \frac{h}{2} \left(f(x_i,y_i) + f(x_{i+1},y_{i+1})\right)
\end{equation*}
Данное формула получения \(y_{i+1}\) является неявной, поэтому находим
его приближенное значение в два шага --- сначала вычисляем значение
\(\widetilde{y}_{i+1}\), затем уточненное \(y_{i+1}\):
\begin{equation}
   \begin{aligned}
       & \widetilde{y}_{i+1} = y_i + hf(x_i,y_i)                                               \\
       & y_{i+1} = y_i + \frac{h}{2} \left(f(x_i,y_i) + f(x_{i+1}, \widetilde{y}_{i+1})\right)
   \end{aligned}
   \label{formula:mod_euler_iterations}
\end{equation}

Данный метод имеет меньшую погрешность, чем предыдущий --- \(O(h^3)\)
на каждой итерации.
\subsubsection{Реализации метода в библиотеках numpy, scipy}
Реализаций данного метода в библиотеках numpy, scipy не найдено.

\subsection{Метод Рунге-Кутта}
\subsubsection{Описание метода}
Формулы \ref{formula:mod_euler_iterations} приводим к следующему виду:
\begin{equation*}
   \begin{aligned}
       & y_{i+1} = y_i + (k_0 + k_1) / 2, \\
       & k_0 = hf(x_i,y_i),               \\
       & k_1 = hf(x_i + h,y_i + h)
   \end{aligned}
\end{equation*}
Данная формула представляет из себя метод Рунге-Кутта второго порядка.
От величины порядка зависит точность полученного решения. Наиболее часто
встречающийся метод четвертого порядка имеет следующий вид:
\begin{equation*}
   \begin{aligned}
       & y_{i+1} = y_i + (k_0 + 2k_1 + 2k_2 + k_3)/6 \\
       & k_0 = hf(x_i,y_i)                           \\
       & k_1 = hf(x_i + h/2,y_i + k_0/2)             \\
       & k_2 = hf(x_i + h/2,y_i + k_1/2)             \\
       & k_3 = hf(x_i + h,y_i + k_2)                 \\
   \end{aligned}
\end{equation*}
\subsubsection{Реализации метода в библиотеках numpy, scipy}
В библиотеке scipy реализована модификации явного метода Рунге-Кутта ---
явные методы Рунге-Кутта-Фельберга \cite{article:fehlberg}, в функции
\textbf{solve\_ivp} модуля \textbf{scipy.integrate}.
Метод Рунге-Кутта-Фельберга порядка \(n(m)\) нужно понимать как метод
Рунге-Кутта порядка \(n\) с погрешностью \(O(h^m)\).

Дополнительно, существуют классы для низкоуровневого управления
вычислениями:
\begin{enumerate}
   \item \textbf{RK23} --- метод Рунге-Кутта-Фельберга 2 порядка с погрешностью \(O(h^3)\).
   \item \textbf{RK45} --- метод Рунге-Кутта-Фельберга 4 порядка с погрешностью \(O(h^5)\).
   \item \textbf{DOP853} --- метод Рунге-Кутта 8 порядка.
\end{enumerate}

Функция \textbf{solve\_ivp} может найти решение задачи Коши для системы
ОДУ, и включает в себя реализации нескольких методов. Далее будут
описаны только те значения параметров, которые необходимы для решения
задачи исследуемым методом (Рунге-Кутта). Также стоит учесть, что при
описании параметров, вместо \(x\) будет использоваться \(t\), как и
принято в зарубежных источниках.
Данная функция имеет следующие параметры:
\begin{enumerate}
   \item \(fun\) --- callable

         Правая часть системы: производная \(\frac{dy(t)}{dt}\).
         Сигнатурой вызова является \(fun(t, y)\), где \(t\) --- скаляр,
         а \(y\) --- массив ndarray с \verb|len(y) = len(y0)|. \(fun\)
         должен  возвращать массив той же размерности, что и \(y\).
         См. \(vectorized\) для получения более детальной информации.
   \item \(t\_span\) --- пара значений float

         Интервал интегрирования \((t0, tf)\). Решатель (\(method\)) начинает выполнение с \(t=t0\) и осуществляет интегрирование, пока не выполнится условие \(t=tf\). И \(t0\), и \(tf\) должны быть числами с плавающей запятой или значениями, интерпретируемыми функцией преобразования чисел с плавающей запятой.
   \item \(y0\) --- array\_like формы (n,)

         Начальное состояние. Для задач на комплексной плоскости необходимо передавать комплексные \(y0\) (даже если начальное значение чисто вещественное).
   \item \(method\) --- string / \verb|OdeSolver|, необязательный

         Используемый метод интеграции:
         \begin{itemize}
            \item \verb|"RK45"| (по умолчанию): Явный метод Рунге-Кутта порядка 5(4). Погрешность контролируется в предположении точности метода четвертого порядка, но шаги выполняются с использованием формулы точности пятого порядка (проводится локальная экстраполяция). При включенном \(dense\_output\) используется интерполяционный полином четвертой степени. Может применяться на комплексной плоскости.

            \item \verb|"RK23"|: Явный метод Рунге-Кутта порядка 3(2). Погрешность контролируется в предположении точности метода второго порядка, но шаги выполняются с использованием формулы точности третьего порядка (проводится локальная экстраполяция). Для плотного вывода используется кубический полином Эрмита. Может применяться на комплексной плоскости.

            \item \verb|"DOP853"|: Явный метод Рунге-Кутта восьмого порядка. Является Python-реализацией алгоритма "DOP853", первоначально написанного на FORTRAN. При включенном \(dense\_output\) используется интерполяционный полином 7-го порядка с точностью до 7-го порядка. Может применяться на комплексной плоскости.

            \item \verb|"Radau"|: Неявный метод Рунге-Кутта семейства Radau IIA порядка 5. Погрешность контролируется с помощью встроенной формулы третьего порядка точности. Кубический полином, который удовлетворяет условиям коллокация, используется при включенном \(dense\_output\).
         \end{itemize}
         Явные методы Рунге-Кутта (\verb|"RK23"|, \verb|"RK45"|, \verb|"DOP853"|) следует использовать для нежестких уравнений, неявные методы (\verb|"Radau"|) --- для жестких. Среди методов Рунге-Кутта для решения с высокой точностью (низкие значения \(rtol\) и \(atol\)) рекомендуется \verb|"DOP853"|.

         Если не уверены, сначала попробуйте запустить \verb|"RK45"|. Если он делает необычно много итераций, расходится или терпит неудачу, ваша проблема, вероятно, будет сложной, и вам следует использовать \verb|"Radau"|.

         Вы также можете передать произвольный класс, производный от \(OdeSolver\), который реализует решатель.
   \item \(t\_eval\)--- array\_like / \verb|None|, необязательный

         Значения \(t\), для которых нужно сохранить вычисленные значения решения, должны быть отсортированы и находиться в пределах \(t\_span\). Если \verb|None| (по умолчанию), используются точки, выбранные решателем.
   \item \(dense\_output\) --- bool, необязательный

         Определяет, следует ли вычислять непрерывное решение. По умолчанию --- \verb|False|.
   \item \(events\) --- callable / list из callable, необязательный

         События для отслеживания. Если \verb|None| (по умолчанию),
         события отслеживаться не будут. Событие происходит, когда
         какая-либо функция, переданная в этом параметре, равна 0.
         Каждая функция должна иметь сигнатуру \(event(t, y)\) и
         возвращать float. Решатель найдет точное значение \(t\), при
         котором \(event(t, y(t)) = 0\), используя алгоритм поиска
         корня. По умолчанию будут найдены все нули. Решатель ищет
         смену знака на каждом шаге, поэтому, если в течение одного
         шага происходит несколько пересечений нуля, события могут быть
         пропущены. Кроме того, каждая функция \(event\) может иметь
         следующие атрибуты:
         \begin{itemize}
            \item \(terminal\): bool, необязательный

                  Определяет, следует ли прекратить интегрирование, если
                  произойдет это событие. По умолчанию --- \verb|False|.
            \item \(direction\): float, необязательный

                  Направление пересечения нуля. Если направление
                  положительное, событие сработает только при переходе
                  от отрицательного к положительному и наоборот, если
                  направление отрицательное. Если 0, то любое
                  направление вызовет появление событие. По умолчанию
                  --- 0.
         \end{itemize}
         Вы можете назначить атрибуты, например
         \verb|event.terminal = True|, любой функции в Python.
   \item \(vectorized\) --- bool, необязательный

         Определяет, можно ли вызвать \(fun\) векторизованным образом.
         По умолчанию --- \verb|False|.
         Если \(vectorized\) имеет значение \verb|False|, \(fun\) всегда
         будет вызываться с \(y\) формы \((n,)\), где \verb|n = len(y0)|.

         Если векторизация имеет значение \verb|True|, \(fun\) можно
         вызвать с помощью \(y\) формы \((n, k)\), где \(k\) --- целое
         число. В этом случае \(fun\) должно вести себя так, чтобы
         \verb|fun(t, y)[:, i] == fun(t, y[:, i])| (то есть каждый
         столбец возвращаемого массива является производной \(dt/dy\)
         соответствующего столбца \(y\)).

         \(vectorized\)=\verb|True| позволяет быстрее аппроксимировать
         якобиан конечной разностью методом \verb|"Radau"|, но в
         некоторых случаях приводит к более медленному выполнению
         других методов, в том числе и для \verb|"Radau"| при некоторых
         условиях (например, малое значение \verb|len(y0)|).
   \item \(args\) --- tuple, необязательный

         Дополнительные аргументы для передачи пользовательским функциям.
         Если заданы, дополнительные аргументы передаются всем
         пользовательским функциям. Так, если, например, \(fun\) имеет
         сигнатуру \(fun(t, y, a, b, c)\), то \(jac\) (если задан) и
         любые функции обработки событий должны иметь одинаковую
         сигнатуру, а \(args\) должен быть tuple длины 3.
         По умолчанию --- \verb|None|.
   \item \(**options\)

         Опции, которые передаются выбранному решателю. Все опции,
         доступные для уже реализованных решателей, перечислены ниже.
         \begin{itemize}
            \item \(first_step\) --- float / None, необязательный

                  Начальный размер шага. По умолчанию установлено
                  значение \verb|None|, что означает, что его должен
                  выбирать решатель.
            \item \(max\_step\) --- float, optional

                  Максимально допустимый размер шага. По умолчанию
                  используется \verb|np.inf|, то есть размер шага не
                  ограничен и определяется исключительно решателем.
            \item \(rtol\), \(atol\) --- float / array\_like, необязательный

                  Относительные и абсолютные погрешности при
                  вычислениях. Решатель сохраняет локальные оценки
                  погрешности меньше, чем \(atol + rtol * abs(y)\).
                  Здесь \(rtol\) контролирует относительную точность
                  (количество правильных цифр), а \(atol\) ---
                  абсолютную точность (количество правильных десятичных
                  знаков). Чтобы достичь желаемого значения \(rtol\),
                  установите значение \(atol <\)
                  \verb|min(rtol * abs(y))|, чтобы значение \(rtol\)
                  доминировало над допустимой ошибкой. Если \(atol\)
                  больше, чем \verb|rtol * abs(y)|, количество
                  правильных цифр не гарантируется. И наоборот, чтобы
                  получить желаемый \(atol\), установите \(rtol\) так,
                  чтобы \(atol <\) \verb|rtol * abs(y)|. Если значения
                  \(y\) имеют разные масштабы, возможно, было бы
                  полезно установить разные значения \(atol\) для
                  разных компонентов, передав array\_like с формой
                  \((n,)\) для \(atol\). Значения по умолчанию:
                  \(1e-3\) для \(rtol\) и \(1e-6\) для \(atol\).
            \item \(jac\) --- array\_like / sparse\_matrix / callable / \verb|None|, необязательный

                  Матрица Якоби правой части системы по y, необходимая
                  для методов \verb|"Radau"|. Матрица Якоби имеет форму
                  \((n, n)\) и ее элемент \((i, j)\) равен
                  \verb|d f_i / d y_j.|. Есть три способа определения
                  матрицы Якоби:
                  \begin{itemize}
                     \item Если array\_like или sparse\_matrix, матрица
                           Якоби считается постоянной.
                     \item Если callable, предполагается, что она
                           зависит как от \(t\), так и от \(y\); при
                           необходимости ее значение будет равно
                           возвращаемому значению \(jac(t, y)\). Для
                           метода \verb|"Radau"| возвращаемое значение
                           может быть разреженной матрицей.
                     \item Если \verb|None| (по умолчанию), матрица
                           Якоби будет аппроксимироваться конечными
                           разностями.
                  \end{itemize}
                  Обычно рекомендуется явно указывать матрицу Якоби,
                  а не полагаться на конечно-разностное приближение.
            \item \(jac\_sparsity\) --- array\_like / sparse matrix / \verb|None|, необязательный

                  Определяет структуру разреженности матрицы Якоби для
                  конечно-разностного приближения. Его форма должна быть
                  \((n, n)\). Этот аргумент игнорируется, если
                  \(jac \ne\) \verb|None|. Если матрица Якоби имеет
                  лишь несколько ненулевых элементов в каждой строке,
                  обеспечение разреженной структуры значительно ускорит
                  вычисления. Нулевая запись означает, что
                  соответствующий элемент матрицы Якоби всегда равен
                  нулю. Если \verb|None| (по умолчанию), матрица Якоби
                  считается плотной.
         \end{itemize}
\end{enumerate}
\chapter{Экспериментальное исследование возможностей библиотек}

\chapter*{Заключение}
\addcontentsline{toc}{chapter}{Заключение}

\input{sources}
\include{appendix}
\end{document}
