1930 lines
128 KiB
TeX
1930 lines
128 KiB
TeX
% 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}, обе из которых могут
|
||
решать задачи для n-диагональных матриц.
|
||
|
||
Различие между ними заключается в том, что \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{equation*}
|
||
\left(
|
||
\begin{aligned}
|
||
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{aligned}
|
||
\right)
|
||
\end{equation*}
|
||
верхняя форма будет следующей:
|
||
|
||
\begin{equation*}
|
||
\left(
|
||
\begin{aligned}
|
||
0 && 0 && -1 && -1 && -1 \\
|
||
0 && 2 && 2 && 2 && 2 \\
|
||
5 && 4 && 3 && 2 && 1 \\
|
||
1 && 1 && 1 && 1 && 0
|
||
\end{aligned}
|
||
\right)
|
||
\end{equation*}
|
||
|
||
Так как данная матрица не эрмитова, и, следовательно, не положительно
|
||
определенна, нижняя форма для нее не существует. Если взять эрмитову
|
||
положительно определенную матрицу
|
||
|
||
\begin{equation*}
|
||
\left(
|
||
\begin{aligned}
|
||
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{aligned}
|
||
\right)
|
||
\end{equation*}
|
||
то ее нижняя форма будет следующая:
|
||
|
||
\begin{equation*}
|
||
\left(
|
||
\begin{aligned}
|
||
4 && 5 && 6 && 7 && 8 && 9 \\
|
||
2 && 2 && 2 && 2 && 2 && 0 \\
|
||
-1 && -1 && -1 && -1 && 0 && 0
|
||
\end{aligned}
|
||
\right)
|
||
\end{equation*}
|
||
|
||
\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}
|
||
\input{experimental_research}
|
||
\input{conclusion}
|
||
\addcontentsline{toc}{chapter}{Заключение}
|
||
|
||
\input{sources}
|
||
\include{appendix}
|
||
\end{document}
|