nm-coursework/main.tex

\input{vars}
\input{config}

\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 --- ndarray из значений bool. Указывает, какие элементы сошлись успешно.
            \item zero\_der --- 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{Сплайн-интерполяция}
\subsubsection{Описание метода}
Между парами значений \((x_{i-1},y_{i-1}),(x_i,y_i)\) строятся
функции-сплайны \(f_1(x),\dots,f_n(x)\), соответствущие условиям
непрерывности. Комбинация таких функций образует исходную:
\begin{equation*}
   \widetilde{y}(x) = \left\{
   \begin{aligned}
       & f_1(x),\; x_0 \le x < x_1     \\
       & f_2(x),\; x_1 \le x < x_2     \\
       & \ \dots \qquad \quad \dots    \\
       & f_n(x),\; x_{n-1} \le x < x_n \\
   \end{aligned}
   \right.
\end{equation*}

Иногда выбирают сплайны такого вида, что непрерывными будут и их
производные, вплоть до нужного порядка (например, если в роли сплайнов
выбрать полиномы 3 степени, то они будут соответствовать условиями
непрерывности вместе с 1 и 2 производными). В зависимости от цели,
на сплайны могут быть наложены и иные ограничения.
\subsubsection{Реализации метода в библиотеках numpy, scipy}
В библиотеке scipy существует множество реализаций данного метода,
которые отличаются используемым видом сплайна.
\subsection{Сглаживание. Метод наименьших квадратов}
\subsubsection{Описание метода}
\subsubsection{Реализации метода в библиотеках numpy, scipy}

\section{Численное интегрирование}
\subsection{Метод прямоугольников}
\subsubsection{Описание метода}
\subsubsection{Реализации метода в библиотеках numpy, scipy}

\subsection{Метод трапеций}
\subsubsection{Описание метода}
\subsubsection{Реализации метода в библиотеках numpy, scipy}

\subsection{Метод парабол (Симпсона)}
\subsubsection{Описание метода}
\subsubsection{Реализации метода в библиотеках numpy, scipy}


\section{Численное решение обыкновенных дифференциальных уравнений}
\subsection{Метод Эйлера}
\subsubsection{Описание метода}
\subsubsection{Реализации метода в библиотеках numpy, scipy}

\subsection{Модифицированный метод Эйлера}
\subsubsection{Описание метода}
\subsubsection{Реализации метода в библиотеках numpy, scipy}

\subsection{Метод Рунге-Кутта}
\subsubsection{Описание метода}
\subsubsection{Реализации метода в библиотеках numpy, scipy}


\chapter{Экспериментальное исследование возможностей библиотек}

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

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