[content] Fix dashes, terminology. Add ODE methods and theory

This commit is contained in:
AVAtarMod 2023-10-07 22:05:30 +03:00
parent eff5d82eb6
commit 4e0737ce02
Signed by: stud128245
GPG Key ID: 43198AE4D0774328

320
main.tex
View File

@ -1273,7 +1273,7 @@ LU-разложение \cite[с. 259]{book:levitin}. Для получения
\verb|False|. \verb|False|.
\item \(nan\_policy\) --- {\verb|"raise"|, \verb|"omit"|, \verb|None|}, необязательный \item \(nan\_policy\) --- {\verb|"raise"|, \verb|"omit"|, \verb|None|}, необязательный
Определяет, как действовать, если входные данные содержат nan. Доступны следующие параметры (по умолчанию \verb|None|): Определяет, как действовать, если входные данные содержат nan. Доступны следующие параметры (по умолчанию --- \verb|None|):
\begin{itemize} \begin{itemize}
\item \verb|"raise"|: выдает ошибку \item \verb|"raise"|: выдает ошибку
@ -1422,7 +1422,7 @@ LU-разложение \cite[с. 259]{book:levitin}. Для получения
Метод решения подзадач доверительной области, применим только для методов \verb|"trf"| и \verb|"dogbox"|. Метод решения подзадач доверительной области, применим только для методов \verb|"trf"| и \verb|"dogbox"|.
\begin{itemize} \begin{itemize}
\item \verb|"exact"| подходит для не очень больших задач с плотными матрицами Якоби. Вычислительная сложность на итерацию сравнима с сингулярным разложением матрицы Якобиана. \item \verb|"exact"| подходит для не очень больших задач с плотными матрицами Якоби. Вычислительная сложность на итерацию сравнима с сингулярным разложением матрицы якобиана.
\end{itemize} \end{itemize}
\item \(tr\_options\) --- dict, необязательный \item \(tr\_options\) --- dict, необязательный
@ -1435,7 +1435,7 @@ LU-разложение \cite[с. 259]{book:levitin}. Для получения
\end{itemize} \end{itemize}
\item \(jac\_sparsity\) --- {\verb|None|, array\_like, разреженная матрица}, необязательный \item \(jac\_sparsity\) --- {\verb|None|, array\_like, разреженная матрица}, необязательный
Определяет структуру разреженности матрицы Якобиана для Определяет структуру разреженности матрицы якобиана для
конечно-разностной оценки, ее форма должна быть (m, n). конечно-разностной оценки, ее форма должна быть (m, n).
Если якобиан имеет лишь несколько ненулевых элементов в Если якобиан имеет лишь несколько ненулевых элементов в
каждой строке, обеспечение разреженной структуры значительно каждой строке, обеспечение разреженной структуры значительно
@ -1517,28 +1517,28 @@ f\(x\), если \(f^{(4)}(x) = 0, x \in [a;b] \), в соответствии
I \approx \frac{h}{2} \sum_{i=0}^{n-1} (f(x_i) + f(x_{i+1})) I \approx \frac{h}{2} \sum_{i=0}^{n-1} (f(x_i) + f(x_{i+1}))
\end{equation*} \end{equation*}
\subsubsection{Реализации метода в библиотеках numpy, scipy} \subsubsection{Реализации метода в библиотеках numpy, scipy}
В библиотеке scipy найдена функция \textbf{trapezoid} в модуле В библиотеке scipy найдена функция \textbf{trapezoid} в модуле
\textbf{scipy.integrate}. Она имеет следующие параметры: \textbf{scipy.integrate}. Она имеет следующие параметры:
\begin{enumerate} \begin{enumerate}
\item \(y\) --- array\_like \item \(y\) --- array\_like
Массив входных данных по которым будет вычислятся интеграл Массив входных данных по которым будет вычислятся интеграл
\item \(x\) --- array\_like, необязательный \item \(x\) --- array\_like, необязательный
Массив значений \(x\), соответсвующих \(y\). Если \verb|None| Массив значений \(x\), соответсвующих \(y\). Если \verb|None|
(по умолчанию), то в роли \(x\) будет создан массив равноотстоящих (по умолчанию), то в роли \(x\) будет создан массив равноотстоящих
значений (\(h = dx\)) значений (\(h = dx\))
\item \(dx\) --- scalar, необязательный \item \(dx\) --- scalar, необязательный
Расстояние между соседними значениями \(x\). Значение используется Расстояние между соседними значениями \(x\). Значение используется
когда \(x\)=\verb|None|. По умолчанию --- 1. когда \(x\)=\verb|None|. По умолчанию --- 1.
\item \(axis\) --- int, необязательный \item \(axis\) --- int, необязательный
Ось, относительно которой производить вычисление интеграла. По Ось, относительно которой производить вычисление интеграла. По
умолчанию --- \(-1\). умолчанию --- \(-1\).
\end{enumerate} \end{enumerate}
Данная функция возвращает число, если \(y\) одномерный и массив Данная функция возвращает число, если \(y\) одномерный и массив
размерности \((n-1)\) для \(n\)-мерного \(y\). размерности \((n-1)\) для \(n\)-мерного \(y\).
\subsection{Метод парабол (Симпсона)} \subsection{Метод парабол (Симпсона)}
\subsubsection{Описание метода} \subsubsection{Описание метода}
@ -1602,20 +1602,306 @@ f\(x\), если \(f^{(4)}(x) = 0, x \in [a;b] \), в соответствии
Стоит учесть, что если ось, по которой необходимо интегрировать, имеет Стоит учесть, что если ось, по которой необходимо интегрировать, имеет
только две точки, то интегрирование происходит с помощью метода трапеций. только две точки, то интегрирование происходит с помощью метода трапеций.
\section{Численное решение обыкновенных дифференциальных уравнений} \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{Метод Эйлера} \subsection{Метод Эйлера}
Данный метод, как и его модифицированная версия, применяется для
ОДУ 1 порядка. Для решения задач Коши ОДУ порядка \(n\) необходимо
исходное уравнение свести к системе ОДУ 1 порядка с помощью замены переменных.
\subsubsection{Описание метода} \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} \subsubsection{Реализации метода в библиотеках numpy, scipy}
Реализаций данного метода в библиотеках numpy, scipy не найдено.
\subsection{Модифицированный метод Эйлера} \subsection{Модифицированный метод Эйлера}
\subsubsection{Описание метода} \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} \subsubsection{Реализации метода в библиотеках numpy, scipy}
Реализаций данного метода в библиотеках numpy, scipy не найдено.
\subsection{Метод Рунге-Кутта} \subsection{Метод Рунге-Кутта}
\subsubsection{Описание метода} \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} \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{Экспериментальное исследование возможностей библиотек}
\chapter*{Заключение} \chapter*{Заключение}