1103 lines
63 KiB
TeX
1103 lines
63 KiB
TeX
\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 --- 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{Сплайн-интерполяция}
|
||
Перед рассмотрением метода необходимо ввести понятие интерполяционного
|
||
сплайна. Это кусочно-заданная функция, каждый фрагмент которой является
|
||
непрерывной функцией, и на концах его значение совпадает с известными
|
||
значениями функции.
|
||
|
||
Иногда выбирают фрагменты сплайна такого вида, что непрерывными будут и их
|
||
производные, вплоть до нужного порядка (например, если в роли сплайнов
|
||
выбрать полиномы 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}, используется для многомерной
|
||
\item интерполяции;
|
||
\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} имеет следующие
|
||
параметры:
|
||
\begin{enumerate}
|
||
\item
|
||
\end{enumerate}
|
||
|
||
Конструктор экземпляра класса \textbf{CubicHermiteSpline} имеет следующие
|
||
параметры:
|
||
\begin{enumerate}
|
||
\item
|
||
\end{enumerate}
|
||
|
||
Конструктор экземпляра класса \textbf{Akima1DInterpolator} имеет
|
||
следующие параметры:
|
||
\begin{enumerate}
|
||
\item
|
||
\end{enumerate}
|
||
|
||
|
||
\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}
|