From 4e0737ce02a596d55f53d40ac254d7aef9b7ce5f Mon Sep 17 00:00:00 2001 From: AVAtarMod Date: Sat, 7 Oct 2023 22:05:30 +0300 Subject: [PATCH] [content] Fix dashes, terminology. Add ODE methods and theory --- main.tex | 320 ++++++++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 303 insertions(+), 17 deletions(-) diff --git a/main.tex b/main.tex index afd88e5..b9e55ee 100644 --- a/main.tex +++ b/main.tex @@ -1273,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"|: выдает ошибку @@ -1422,7 +1422,7 @@ LU-разложение \cite[с. 259]{book:levitin}. Для получения Метод решения подзадач доверительной области, применим только для методов \verb|"trf"| и \verb|"dogbox"|. \begin{itemize} - \item \verb|"exact"| подходит для не очень больших задач с плотными матрицами Якоби. Вычислительная сложность на итерацию сравнима с сингулярным разложением матрицы Якобиана. + \item \verb|"exact"| подходит для не очень больших задач с плотными матрицами Якоби. Вычислительная сложность на итерацию сравнима с сингулярным разложением матрицы якобиана. \end{itemize} \item \(tr\_options\) --- dict, необязательный @@ -1435,7 +1435,7 @@ LU-разложение \cite[с. 259]{book:levitin}. Для получения \end{itemize} \item \(jac\_sparsity\) --- {\verb|None|, array\_like, разреженная матрица}, необязательный - Определяет структуру разреженности матрицы Якобиана для + Определяет структуру разреженности матрицы якобиана для конечно-разностной оценки, ее форма должна быть (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})) \end{equation*} \subsubsection{Реализации метода в библиотеках numpy, scipy} -В библиотеке scipy найдена функция \textbf{trapezoid} в модуле +В библиотеке scipy найдена функция \textbf{trapezoid} в модуле \textbf{scipy.integrate}. Она имеет следующие параметры: \begin{enumerate} \item \(y\) --- array\_like - - Массив входных данных по которым будет вычислятся интеграл + + Массив входных данных по которым будет вычислятся интеграл \item \(x\) --- array\_like, необязательный - - Массив значений \(x\), соответсвующих \(y\). Если \verb|None| - (по умолчанию), то в роли \(x\) будет создан массив равноотстоящих - значений (\(h = dx\)) + + Массив значений \(x\), соответсвующих \(y\). Если \verb|None| + (по умолчанию), то в роли \(x\) будет создан массив равноотстоящих + значений (\(h = dx\)) \item \(dx\) --- scalar, необязательный - - Расстояние между соседними значениями \(x\). Значение используется - когда \(x\)=\verb|None|. По умолчанию --- 1. + + Расстояние между соседними значениями \(x\). Значение используется + когда \(x\)=\verb|None|. По умолчанию --- 1. \item \(axis\) --- int, необязательный - Ось, относительно которой производить вычисление интеграла. По - умолчанию --- \(-1\). + Ось, относительно которой производить вычисление интеграла. По + умолчанию --- \(-1\). \end{enumerate} -Данная функция возвращает число, если \(y\) одномерный и массив +Данная функция возвращает число, если \(y\) одномерный и массив размерности \((n-1)\) для \(n\)-мерного \(y\). \subsection{Метод парабол (Симпсона)} \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{Метод Эйлера} +Данный метод, как и его модифицированная версия, применяется для +ОДУ 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*{Заключение}