[content] Add trapezoidal rule, simpson rule
This commit is contained in:
AVAtarMod
parent
fd0d028c97
commit
55f7077896
GPG Key ID:
43198AE4D0774328
195
main.tex
195
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, необязательный
|
||||
|
||||
Определяет относительный размер шага для аппроксимации якобиана конечной разностью. Фактический шаг вычисляется как \verb|x * diff_step|. Если \verb|None| (по умолчанию), то \(diff_step\) принимается за обычную "оптимальную" степень машинного эпсилона для используемой схемы конечных разностей.
|
||||
\item \(tr\_solver\) --- {\verb|None|, \verb|"exact"|, \verb|"lsmr"|}, необязательный
|
||||
|
||||
\item
|
||||
Метод решения подзадач доверительной области, применим только для методов \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}
|
||||
\section{Численное интегрирование}
|
||||
\subsection{Метод прямоугольников}
|
||||
\subsubsection{Описание метода}
|
||||
\subsubsection{Реализации метода в библиотеках numpy, scipy}
|
||||
Функция возвращает результат в виде экземпляра класса \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{Метод Эйлера}
|
||||
|
||||
Loading…
Reference in New Issue
Block a user