Compare commits
4 Commits
e885a649b2
...
4e0737ce02
| Author | SHA1 | Date | |
|---|---|---|---|
|
4e0737ce02
|
||
|
|
eff5d82eb6
|
||
|
|
55f7077896
|
||
|
|
fd0d028c97
|
489
main.tex
489
main.tex
@ -114,10 +114,10 @@
|
||||
противоположные знаки.
|
||||
\item \(a\) --- scalar
|
||||
|
||||
Первый конец интервала \([a,b]\).
|
||||
Первый конец интервала \([a;b]\).
|
||||
\item \(b\) --- scalar
|
||||
|
||||
Второй конец интервала \([a,b]\).
|
||||
Второй конец интервала \([a;b]\).
|
||||
\item \(xtol\) --- number, необязательный
|
||||
|
||||
Вычисленный корень \(x0\) будет удовлетворять
|
||||
@ -797,7 +797,7 @@ LU-разложение \cite[с. 259]{book:levitin}. Для получения
|
||||
последний --- с помощью сглаживания.
|
||||
\subsection{Интерполяционный полином Лагранжа}
|
||||
\subsubsection{Описание метода}
|
||||
Для известных пар значений \((x_i,y_i), x_i \in [a,b]; i = 1,2,3,\dots,n\)
|
||||
Для известных пар значений \((x_i,y_i), x_i \in [a;b]; i = 1,2,3,\dots,n\)
|
||||
строится полином \(P(x)\), удовлетворяющий условию \(P(x_i) = y_i\).
|
||||
|
||||
Полином \(P(x)\) определяется следующей формулой:
|
||||
@ -1164,6 +1164,11 @@ LU-разложение \cite[с. 259]{book:levitin}. Для получения
|
||||
минимизации выражения \(||ax-b||^2\), где \(a,b\) --- матрицы
|
||||
параметров; функция \textbf{nnls} отличается от нее
|
||||
ограничением на положительность искомых значений.
|
||||
|
||||
Подходят только для нахождения параметров многочленов,
|
||||
не буду рассматриваться в данной работе по причине
|
||||
сложной подготовки входных данных для решения задачи
|
||||
метода, а также их узкой специализации.
|
||||
\end{enumerate}
|
||||
|
||||
Функция \textbf{curve\_fit} имеет следующий набор параметров:
|
||||
@ -1268,7 +1273,7 @@ LU-разложение \cite[с. 259]{book:levitin}. Для получения
|
||||
\verb|False|.
|
||||
\item \(nan\_policy\) --- {\verb|"raise"|, \verb|"omit"|, \verb|None|}, необязательный
|
||||
|
||||
Определяет, как действовать, если входные данные содержат nan. Доступны следующие параметры (по умолчанию — \verb|None|):
|
||||
Определяет, как действовать, если входные данные содержат nan. Доступны следующие параметры (по умолчанию --- \verb|None|):
|
||||
\begin{itemize}
|
||||
\item \verb|"raise"|: выдает ошибку
|
||||
|
||||
@ -1338,7 +1343,7 @@ LU-разложение \cite[с. 259]{book:levitin}. Для получения
|
||||
будет рассматриваться как одномерный массив с одним элементом.
|
||||
\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|.
|
||||
Метод вычисления матрицы Якоби (матрица размером \(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|, необязательный
|
||||
|
||||
Есть два способа указать границы:
|
||||
@ -1408,41 +1413,495 @@ LU-разложение \cite[с. 259]{book:levitin}. Для получения
|
||||
\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"| учитывает вызовы функций в оценке Якобиана).
|
||||
\item Для \verb|"lm"|: \(100 \cdot n\), если \(jac\) является вызываемым, и \(100 \cdot n \cdot (n + 1)\) в противном случае (поскольку \verb|"lm"| учитывает вызовы функций в оценке якобиана).
|
||||
\end{itemize}
|
||||
\item \(diff\_step\) --- \verb|None| / array\_like, необязательный
|
||||
|
||||
|
||||
\item
|
||||
\end{enumerate}
|
||||
\section{Численное интегрирование}
|
||||
\subsection{Метод прямоугольников}
|
||||
\subsubsection{Описание метода}
|
||||
\subsubsection{Реализации метода в библиотеках numpy, scipy}
|
||||
|
||||
Определяет относительный размер шага для аппроксимации якобиана конечной разностью. Фактический шаг вычисляется как \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
|
||||
|
||||
\section{Численное решение обыкновенных дифференциальных уравнений}
|
||||
Массив, который необходимо интегрировать.
|
||||
\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 из callables, необязательный
|
||||
|
||||
События для отслеживания. Если \verb|None| (по умолчанию),
|
||||
события отслеживаться не будут. Событие происходит, когда
|
||||
какая-либо функция, переданная в этом параметре, равна 0.
|
||||
Каждая функция должна иметь сигнатуру \(event(t, y)\) и
|
||||
возвращать float. Решатель найдет точное значение \(t\), при
|
||||
котором \(event(t, y(t)) = 0\), используя алгоритм поиска
|
||||
корня. По умолчанию будут найдены все нули. Решатель ищет
|
||||
смену знака на каждом шаге, поэтому, если в течение одного
|
||||
шага происходит несколько пересечений нуля, события могут быть
|
||||
пропущены. Кроме того, каждая функция \(event\) может иметь
|
||||
следующие атрибуты:
|
||||
\begin{itemize}
|
||||
\item \(terminal\): bool, необязательный
|
||||
|
||||
Определяет, следует ли прекратить интегрирование, если
|
||||
произойдет это событие. По умолчанию --- \verb|False|.
|
||||
\item \(direction\): float, необязательный
|
||||
|
||||
Направление пересечения нуля. Если направление
|
||||
положительное, событие сработает только при переходе
|
||||
от отрицательного к положительному и наоборот, если
|
||||
направление отрицательное. Если 0, то любое
|
||||
направление вызовет появление событие. По умолчанию
|
||||
--- 0.
|
||||
\end{itemize}
|
||||
Вы можете назначить атрибуты, например
|
||||
\verb|event.terminal = True|, любой функции в Python.
|
||||
\item \(vectorized\) --- bool, необязательный
|
||||
|
||||
Определяет, можно ли вызвать \(fun\) векторизованным образом.
|
||||
По умолчанию --- \verb|False|.
|
||||
Если \(vectorized\) имеет значение \verb|False|, \(fun\) всегда
|
||||
будет вызываться с \(y\) формы \((n,)\), где \verb|n = len(y0)|.
|
||||
|
||||
Если векторизация имеет значение \verb|True|, \(fun\) можно
|
||||
вызвать с помощью \(y\) формы \((n, k)\), где \(k\) --- целое
|
||||
число. В этом случае \(fun\) должно вести себя так, чтобы
|
||||
\verb|fun(t, y)[:, i] == fun(t, y[:, i])| (то есть каждый
|
||||
столбец возвращаемого массива является производной \(dt/dy\)
|
||||
соответствующего столбца \(y\)).
|
||||
|
||||
\(vectorized\)=\verb|True| позволяет быстрее аппроксимировать
|
||||
якобиан конечной разностью методом \verb|"Radau"|, но в
|
||||
некоторых случаях приводит к более медленному выполнению
|
||||
других методов, в том числе и для \verb|"Radau"| при некоторых
|
||||
условиях (например, малое значение \verb|len(y0)|).
|
||||
\item \(args\) --- tuple, необязательный
|
||||
|
||||
Дополнительные аргументы для передачи пользовательским функциям.
|
||||
Если заданы, дополнительные аргументы передаются всем
|
||||
пользовательским функциям. Так, если, например, \(fun\) имеет
|
||||
сигнатуру \(fun(t, y, a, b, c)\), то \(jac\) (если задан) и
|
||||
любые функции обработки событий должны иметь одинаковую
|
||||
сигнатуру, а \(args\) должен быть tuple длины 3.
|
||||
По умолчанию --- \verb|None|.
|
||||
\item \(**options\)
|
||||
|
||||
Опции, которые передаются выбранному решателю. Все опции,
|
||||
доступные для уже реализованных решателей, перечислены ниже.
|
||||
\begin{itemize}
|
||||
\item \(first_step\) --- float / None, необязательный
|
||||
|
||||
Начальный размер шага. По умолчанию установлено
|
||||
значение \verb|None|, что означает, что его должен
|
||||
выбирать решатель.
|
||||
\item \(max\_step\) --- float, optional
|
||||
|
||||
Максимально допустимый размер шага. По умолчанию
|
||||
используется \verb|np.inf|, то есть размер шага не
|
||||
ограничен и определяется исключительно решателем.
|
||||
\item \(rtol\), \(atol\) --- float / array\_like, необязательный
|
||||
|
||||
Относительные и абсолютные погрешности при
|
||||
вычислениях. Решатель сохраняет локальные оценки
|
||||
погрешности меньше, чем \(atol + rtol * abs(y)\).
|
||||
Здесь \(rtol\) контролирует относительную точность
|
||||
(количество правильных цифр), а \(atol\) ---
|
||||
абсолютную точность (количество правильных десятичных
|
||||
знаков). Чтобы достичь желаемого значения \(rtol\),
|
||||
установите значение \(atol <\)
|
||||
\verb|min(rtol * abs(y))|, чтобы значение \(rtol\)
|
||||
доминировало над допустимой ошибкой. Если \(atol\)
|
||||
больше, чем \verb|rtol * abs(y)|, количество
|
||||
правильных цифр не гарантируется. И наоборот, чтобы
|
||||
получить желаемый \(atol\), установите \(rtol\) так,
|
||||
чтобы \(atol <\) \verb|rtol * abs(y)|. Если значения
|
||||
\(y\) имеют разные масштабы, возможно, было бы
|
||||
полезно установить разные значения \(atol\) для
|
||||
разных компонентов, передав array\_like с формой
|
||||
\((n,)\) для \(atol\). Значения по умолчанию:
|
||||
\(1e-3\) для \(rtol\) и \(1e-6\) для \(atol\).
|
||||
\item \(jac\) --- array\_like / sparse\_matrix / callable / \verb|None|, необязательный
|
||||
|
||||
Матрица Якоби правой части системы по y, необходимая
|
||||
для методов \verb|"Radau"|. Матрица Якоби имеет форму
|
||||
\((n, n)\) и ее элемент \((i, j)\) равен
|
||||
\verb|d f_i / d y_j.|. Есть три способа определения
|
||||
матрицы Якоби:
|
||||
\begin{itemize}
|
||||
\item Если array\_like или sparse\_matrix, матрица
|
||||
Якоби считается постоянной.
|
||||
\item Если callable, предполагается, что она
|
||||
зависит как от \(t\), так и от \(y\); при
|
||||
необходимости ее значение будет равно
|
||||
возвращаемому значению \(jac(t, y)\). Для
|
||||
метода \verb|"Radau"| возвращаемое значение
|
||||
может быть разреженной матрицей.
|
||||
\item Если \verb|None| (по умолчанию), матрица
|
||||
Якоби будет аппроксимироваться конечными
|
||||
разностями.
|
||||
\end{itemize}
|
||||
Обычно рекомендуется явно указывать матрицу Якоби,
|
||||
а не полагаться на конечно-разностное приближение.
|
||||
\item \(jac\_sparsity\) --- array\_like / sparse matrix / \verb|None|, необязательный
|
||||
|
||||
Определяет структуру разреженности матрицы Якоби для
|
||||
конечно-разностного приближения. Его форма должна быть
|
||||
\((n, n)\). Этот аргумент игнорируется, если
|
||||
\(jac \ne\) \verb|None|. Если матрица Якоби имеет
|
||||
лишь несколько ненулевых элементов в каждой строке,
|
||||
обеспечение разреженной структуры значительно ускорит
|
||||
вычисления. Нулевая запись означает, что
|
||||
соответствующий элемент матрицы Якоби всегда равен
|
||||
нулю. Если \verb|None| (по умолчанию), матрица Якоби
|
||||
считается плотной.
|
||||
\end{itemize}
|
||||
\end{enumerate}
|
||||
\chapter{Экспериментальное исследование возможностей библиотек}
|
||||
|
||||
\chapter*{Заключение}
|
||||
|
||||
@ -25,4 +25,7 @@
|
||||
\bibitem{links:scipy_doc} Scipy API Reference [Электронный ресурс] -- URL:~\url{https://docs.scipy.org/doc/scipy/reference/index.html} (\LiteratureAccessDate[3]).
|
||||
\bibitem{links:tiobe_index} TIOBE. Официальный сайт проекта
|
||||
[Электронный ресурс] -- URL:~\url{https://www.tiobe.com/tiobe-index/} (\LiteratureAccessDate).
|
||||
\bibitem{journal:cartwright} Simpson's Rule Cumulative Integration with MS Excel and Irregularly-spaced Data / Cartwright, Kenneth V. // Journal of Mathematical Sciences and Mathematics Education. -- 2017. -- Vol.~12~(2) -- P.~1-9.
|
||||
\bibitem{article:fehlberg} Classical Fifth-, Sixth-, Seventh-, and Eighth-Order Runge-Kutta Formulas with Stepsize Control / Fehlberg E. // NASA technical report 287 -- 1968.
|
||||
-- P.~82
|
||||
\end{thebibliography}
|
||||
|
||||
Loading…
Reference in New Issue
Block a user