From 55f70778965ec9a6c1b5ec1c8c146904287dabcc Mon Sep 17 00:00:00 2001 From: AVAtarMod Date: Fri, 6 Oct 2023 19:08:54 +0300 Subject: [PATCH] [content] Add trapezoidal rule, simpson rule --- main.tex | 201 +++++++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 187 insertions(+), 14 deletions(-) diff --git a/main.tex b/main.tex index 99acdbe..afd88e5 100644 --- a/main.tex +++ b/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} имеет следующий набор параметров: @@ -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,26 +1413,194 @@ 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{Описание метода} -\subsubsection{Реализации метода в библиотеках numpy, scipy} +На отрезке \([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{Численное решение обыкновенных дифференциальных уравнений} \subsection{Метод Эйлера}