diff --git a/appendix.tex b/appendix.tex index 3e1ab1f..af0953b 100644 --- a/appendix.tex +++ b/appendix.tex @@ -14,6 +14,51 @@ \addcontentsline{toc}{chapter}{Приложения} \appendix -% \section{Скрипты установки БД для компонента <<Хранение~данных>>} -% \label{script:storage} -% \inputcode{../software/architecture/first-level/component-storage/CSt.db-script.sql}{frame=none,language=sql} + +\section{Результаты вывода программы} +\label{output_program} +\begin{figure}[H] + \centering + \includegraphics[width=0.6\textwidth]{assets/Thomas} + \caption{Вывод программы для решения СЛУ методом прогонки} +\end{figure} +\begin{figure}[H] + \centering + \includegraphics[width=0.6\textwidth]{assets/Inverted} + \caption{Вывод программы для решения СЛУ методом обратной матрицы} +\end{figure} +\begin{figure}[H] + \centering + \includegraphics[width=0.6\textwidth]{assets/lagrange} + \caption{График интерполяции функции методом Лагранжа} +\end{figure} +\begin{figure}[H] + \centering + \includegraphics[width=0.6\textwidth]{assets/newton} + \caption{График интерполяции функции методом Ньютона} +\end{figure} +\begin{figure}[H] + \centering + \includegraphics[width=0.6\textwidth]{assets/CubicSpline} + \caption{График сплайн-интерполяции функции с помощью \textbf{CubicSpline}} +\end{figure} +\begin{figure}[H] + \centering + \includegraphics[width=0.6\textwidth]{assets/PchipInterpolator} + \caption{График сплайн-интерполяции функции с помощью \textbf{PchipInterpolator}} +\end{figure} +\begin{figure}[H] + \centering + \includegraphics[width=0.6\textwidth]{assets/CubicHermiteSpline} + \caption{График сплайн-интерполяции функции с помощью \textbf{CubicHermiteSpline}} +\end{figure} +\begin{figure}[H] + \centering + \includegraphics[width=0.6\textwidth]{assets/Akima1DInterpolator} + \caption{График сплайн-интерполяции функции с помощью \textbf{Akima1DInterpolator}} +\end{figure} + + +\section{Листинг программы} +\label{program_code} +\inputcode{code/main.py}{frame=none,language=python} diff --git a/assets/Akima1DInterpolator.png b/assets/Akima1DInterpolator.png new file mode 100644 index 0000000..c53ca94 Binary files /dev/null and b/assets/Akima1DInterpolator.png differ diff --git a/assets/CubicHermiteSpline.png b/assets/CubicHermiteSpline.png new file mode 100644 index 0000000..90ea256 Binary files /dev/null and b/assets/CubicHermiteSpline.png differ diff --git a/assets/CubicSpline.png b/assets/CubicSpline.png new file mode 100644 index 0000000..8049602 Binary files /dev/null and b/assets/CubicSpline.png differ diff --git a/assets/Gauss.png b/assets/Gauss.png new file mode 100644 index 0000000..3cff322 Binary files /dev/null and b/assets/Gauss.png differ diff --git a/assets/Inverted.png b/assets/Inverted.png new file mode 100644 index 0000000..313e359 Binary files /dev/null and b/assets/Inverted.png differ diff --git a/assets/PchipInterpolator.png b/assets/PchipInterpolator.png new file mode 100644 index 0000000..915e26b Binary files /dev/null and b/assets/PchipInterpolator.png differ diff --git a/assets/Thomas.png b/assets/Thomas.png new file mode 100644 index 0000000..4df8cef Binary files /dev/null and b/assets/Thomas.png differ diff --git a/assets/bisect.png b/assets/bisect.png new file mode 100644 index 0000000..6b2f736 Binary files /dev/null and b/assets/bisect.png differ diff --git a/assets/curve_fit.png b/assets/curve_fit.png new file mode 100644 index 0000000..c356ba2 Binary files /dev/null and b/assets/curve_fit.png differ diff --git a/assets/lagrange.png b/assets/lagrange.png new file mode 100644 index 0000000..b1ac81e Binary files /dev/null and b/assets/lagrange.png differ diff --git a/assets/least_squares.png b/assets/least_squares.png new file mode 100644 index 0000000..c2fa77e Binary files /dev/null and b/assets/least_squares.png differ diff --git a/assets/newton.png b/assets/newton.png new file mode 100644 index 0000000..519df52 Binary files /dev/null and b/assets/newton.png differ diff --git a/assets/python-logo.png b/assets/python-logo.png new file mode 100644 index 0000000..346f493 Binary files /dev/null and b/assets/python-logo.png differ diff --git a/files/vyatsu_logo.png b/assets/vyatsu_logo.png similarity index 100% rename from files/vyatsu_logo.png rename to assets/vyatsu_logo.png diff --git a/code/main.py b/code/main.py new file mode 100644 index 0000000..9a2141f --- /dev/null +++ b/code/main.py @@ -0,0 +1,429 @@ +import scipy.integrate as sitg +import scipy.interpolate as sitp +import scipy.optimize as sopt +import scipy.linalg as salg +import math as m +import numpy as np + +import matplotlib.pyplot as plt + + +def create_subplot(): + return plt.subplots(layout='constrained')[1] + + +def plt_append(sp, x: list[float], y: list[float], label: str, format: str): + sp.plot(x, y, format, label=label) + + +def generate_array(min, max, density=10): + point_count = int(m.fabs(max-min)*density) + x = np.linspace(min, max, point_count) + return list(x.tolist()) + + +class NonLinear: + bisect_exp = "x**2 * np.sin(x)" + newton_exp = "np.sin(x) * np.sqrt(np.abs(x))" + + @staticmethod + def slice_array(range: list[float], val_min, val_max): + def index_search(range: list[float], val): + i = 0 + for v in range: + if v >= val: + return i + i += 1 + return -1 + + index_l = index_search( + range, val_min) if val_min is not None else range.index(min(range)) + index_r = index_search( + range, val_max) if val_max is not None else range.index(max(range)) + return range[index_l:index_r+1] + + @staticmethod + def bisect(x, x_min, x_max): + def f(x): return eval(NonLinear.bisect_exp) + y = f(np.array(x)) + root = sopt.bisect(f, x_min, x_max) + solution = root[0] if root is tuple else root + return list(y), (float(solution), float(f(solution))) + + @staticmethod + def plot_bisect(): + bounds = 0, 6 + split_val = 1 + x1 = generate_array(bounds[0], bounds[1]) + x2 = NonLinear.slice_array(x1, split_val, None) + + sp = create_subplot() + + sol1 = NonLinear.bisect(x1, bounds[0], bounds[1]) + sol2 = NonLinear.bisect(x2, split_val, bounds[1]) + + plt_append( + sp, x1, sol1[0], f"Исходные данные (y={NonLinear.bisect_exp})", "-b") + plt_append( + sp, *(sol1[1]), f"bisect на [{bounds[0]},{bounds[1]}]", "or") + plt_append( + sp, *(sol2[1]), f"bisect на [{split_val},{bounds[1]}]", "og") + + sp.set_title("scipy.optimize.bisect") + sp.legend(loc='lower left') + + @staticmethod + def newton(x, x0): + def f(x): return eval(NonLinear.bisect_exp) + y = f(np.array(x)) + root = sopt.newton(f, x0) + solution = root[0] if root is tuple else root + return list(y), (float(solution), float(f(solution))) + + @staticmethod + def plot_newton(): + bounds = -2, 7 + split_l, split_r = 2, 5 + x1 = generate_array(bounds[0], bounds[1]) + x2 = NonLinear.slice_array(x1, split_l, split_r) + x0_1, x0_2 = 1/100, 4 + sp = create_subplot() + + sol1 = NonLinear.newton(x1, x0_1) + sol2 = NonLinear.newton(x2, x0_2) + + plt_append( + sp, x1, sol1[0], f"Исходные данные (y={NonLinear.newton_exp})", "-b") + plt_append( + sp, *(sol1[1]), f"newton на отрезке [{bounds[0]},{bounds[1]}]", "or") + plt_append( + sp, *(sol2[1]), f"newton на отрезке [{split_l},{bounds[1]}]", "og") + + sp.set_title("scipy.optimize.newton") + sp.legend(loc='lower left') + + @staticmethod + def plot(method: str = "all"): + if method in ["bisect", "all"]: + NonLinear.plot_bisect() + if method in ["newton", "all"]: + NonLinear.plot_newton() + plt.ylabel("y") + plt.xlabel("x") + plt.show() + + +class SLE: + gauss_data = ([[13, 2], [3, 4]], [1, 2]) + invmatrix_data = ([[13, 2], [3, 4]], [1, 2]) + tridiagonal_data = ([[4, 5, 6, 7, 8, 9], + [2, 2, 2, 2, 2, 0]], + [1, 2, 2, 3, 3, 3]) + + @staticmethod + def var_str(index): + return f"x{index+1}" + + @staticmethod + def print_solution(data: list[float]): + print(" ", end='') + for i, val in enumerate(data[:-1]): + print(f"{SLE.var_str(i)} = {round(val,3)}, ", end='') + print(f"{SLE.var_str(len(data)-1)} = {round(data[-1],3)}") + + @staticmethod + def print_data(data: tuple[list[list[float]], list[float]], tridiagonal: bool = False): + if tridiagonal: + new_data = [] + new_len = len(data[0][0]) + zipped = list(zip(*tuple(data[0]))) + zipped[len(zipped)-1] = (zipped[len(zipped)-1] + [0], zipped[len(zipped)-2][1]) + complement_to = new_len - len(zipped[0]) + for i, val in enumerate(zipped): + zero_r = complement_to - i + if zero_r <= 0: + zero_r = 0 + mid_val = list(reversed(val[1:])) + list(val) + mid_end = len(mid_val) if zero_r > 0 else len( + mid_val) + (complement_to - i) + mid_beg = len(mid_val) - (new_len - zero_r) if zero_r > 0 else 0 + mid_beg = mid_beg if mid_beg >= 0 else 0 + zero_l = new_len - (zero_r + (mid_end - mid_beg)) + tmp = [0] * zero_l + \ + mid_val[mid_beg:mid_end] + [0] * zero_r + new_data.append(tmp) + data = (new_data, data[1]) + for i, val in enumerate(data[0]): + print(" ", end='') + for i_coef, coef in enumerate(val[:-1]): + if coef != 0: + print(f"({coef}{SLE.var_str(i_coef)}) + ", end='') + else: + print(f" {coef} + ", end='') + print(f"({val[-1]}{SLE.var_str(len(val)-1)})", end='') + print(f" = {data[1][i]}") + + @staticmethod + def gauss(system: list[list[float]], b: list[float]): + lup = salg.lu_factor(system) + solution = salg.lu_solve(lup, b) + return solution + + @staticmethod + def invmatrix(system: list[list[float]], b: list[float]): + m_inv = salg.inv(system) + solution = m_inv @ b + return solution + + @staticmethod + def tridiagonal(system: list[list[float]], b: list[float]): + solution = salg.solveh_banded(system, b, lower=True) + return solution + + @staticmethod + def print_gauss(): + print("Gauss method (LU decomposition)") + print(" Input system:") + SLE.print_data(SLE.gauss_data) + print(" Solution:") + SLE.print_solution(SLE.gauss(*SLE.gauss_data)) + + @staticmethod + def print_invmatrix(): + print("Inverted matrix method") + print(" Input system:") + SLE.print_data(SLE.invmatrix_data) + print(" Solution:") + SLE.print_solution(SLE.invmatrix(*SLE.invmatrix_data)) + + @staticmethod + def print_tridiagonal(): + print("Tridiagonal matrix method (Thomas algorithm)") + print(" Input system:") + SLE.print_data(SLE.tridiagonal_data, True) + print(" Solution:") + SLE.print_solution(SLE.tridiagonal(*SLE.tridiagonal_data)) + + @staticmethod + def print(method="all"): + if method in ["gauss", "all"]: + SLE.print_gauss() + if method in ["invmatrix", "all"]: + SLE.print_invmatrix() + if method in ["banded", "all"]: + SLE.print_tridiagonal() + + +class Approx: + function_exp = "np.sin(x) * np.sqrt(np.abs(x))" + least_sq_exp = "np.sin(x) * np.abs(x)" + + @staticmethod + def get_function_exp_der(*args): + function_der_exp = "(x * np.sin(x) + 2 * x**2 * np.cos(x)) / (2 * np.sqrt(np.abs(x)) ** 3)" + result = () + for i in args: + array = [] + for x in i: + array.append(eval(function_der_exp)) + + result = result + (array,) + return result + + @staticmethod + def generate_y(x_array, function): + result = [] + for x in x_array: + result.append(eval(function)) + return result + + @staticmethod + def lagrange(x, y): + return sitp.lagrange(x, y) + + @staticmethod + def get_approx_data(function=function_exp, bounds=[-6, 6]): + x1 = generate_array(bounds[0], bounds[1], 1/2) + x2 = generate_array(bounds[0], bounds[1], 1) + y1 = Approx.generate_y(x1, function) + y2 = Approx.generate_y(x2, function) + x_real = generate_array(bounds[0], bounds[1]) + y_real = Approx.generate_y(x_real, function) + return x1, x2, y1, y2, x_real, y_real + + @staticmethod + def plot_lagrange(): + x1, x2, y1, y2, x_real, y_real = Approx.get_approx_data() + + sp = create_subplot() + sol1 = np.polynomial.polynomial.Polynomial( + Approx.lagrange(x1, y1).coef[::-1]) + sol2 = np.polynomial.polynomial.Polynomial( + Approx.lagrange(x2, y2).coef[::-1]) + + plt_append( + sp, x_real, y_real, f"Исходные данные (y={Approx.function_exp})", "--b") + plt_append( + sp, x_real, sol1(np.array(x_real)), f"f1 = lagrange, кол-во точек = {len(x1)}", "-m") + plt_append( + sp, x_real, sol2(np.array(x_real)), f"f2 = lagrange, кол-во точек = {len(x2)}", "-r") + plt_append( + sp, x1, y1, f"Исходные точки для f1", ".m") + plt_append( + sp, x2, y2, f"Исходные точки для f2", ".r") + + sp.set_title("scipy.interpolate.lagrange") + sp.legend(loc='lower left') + + @staticmethod + def plot_spline(): + x1, x2, y1, y2, x_real, y_real = Approx.get_approx_data() + d1, d2 = Approx.get_function_exp_der(x1, x2) + + for interpolator in [sitp.CubicSpline, + sitp.PchipInterpolator, + sitp.CubicHermiteSpline, + sitp.Akima1DInterpolator]: + sp = create_subplot() + + if interpolator.__name__ != "CubicHermiteSpline": + args1 = x1, y1 + args2 = x2, y2 + else: + args1 = x1, y1, d1 + args2 = x2, y2, d2 + + sol1 = interpolator(*args1) + sol2 = interpolator(*args2) + plt_append( + sp, x_real, y_real, f"Исходные данные (y={Approx.function_exp})", "--b") + plt_append( + sp, x_real, sol1(np.array(x_real)), f"f1 = {interpolator.__name__}, кол-во точек = {len(x1)}", "-m") + plt_append( + sp, x_real, sol2(np.array(x_real)), f"f2 = {interpolator.__name__}, кол-во точек = {len(x2)}", "-r") + plt_append( + sp, x1, y1, f"Исходные точки для f1", ".m") + plt_append( + sp, x2, y2, f"Исходные точки для f2", ".r") + + sp.set_title(f"scipy.interpolate.{interpolator.__name__}") + sp.legend(loc='lower left') + + @staticmethod + def linear(x, a, b): + return a*x + b + + @staticmethod + def quadratic(x, a, b, c): + return a * (x**2) + (b*x) + c + + @staticmethod + def fract(x, a, b, c): + return x / (a * x + b) - c + + @staticmethod + def noise_y(y, rng): + diff = max(y) - min(y) + noise_coeff = diff*(10/100) + return y + (noise_coeff * rng.normal(size=len(y))) + + @staticmethod + def plot_least_squares_curvefit(): + rng = np.random.default_rng() + bounds = [3, 6] + x1, x2, y1, y2, x_real, y_real = Approx.get_approx_data( + Approx.least_sq_exp, bounds) + x_real = np.array(x_real) + + y_real = Approx.noise_y(y_real, rng) + base_functions = [Approx.linear, + Approx.quadratic, (Approx.fract, "x/(ax+b)")] + + sp = create_subplot() + plt_append( + sp, x_real, y_real, f"y={Approx.least_sq_exp} на [{bounds[0]};{bounds[1]}], с шумом", ".b") + for bf in base_functions: + if isinstance(bf, tuple): + bf, desc = bf[0], bf[1] + else: + bf, desc = bf, None + optimal_params, _ = sopt.curve_fit(bf, x_real, y_real) + desc_str = f" ({desc}) " if desc is not None else "" + plt_append( + sp, x_real, bf(np.array(x_real), *optimal_params), + f"МНК, вид функции - {bf.__name__}{desc_str}", "-") + sp.set_title(f"scipy.optimize.curve_fit") + sp.legend(loc='lower left') + + @staticmethod + def plot_least_squares(): + rng = np.random.default_rng() + + def exponential(x, a, b, c): + return np.sin(x) * np.sqrt(np.abs(x)) + exponential.str = "np.sin(x) * np.sqrt(np.abs(x))" + + def gen_y(x, a, b, c, noise=0., n_outliers=0): + y = exponential(x, a, b, c) + error = noise * rng.standard_normal(x.size) + outliers = rng.integers(0, x.size, n_outliers) + error[outliers] *= 10 + return y + error + + def loss(params, x, y): + return (exponential(x, params[0], params[1], params[2])) - y + params0 = np.array([0.1, 1, 0]) + bounds = [-5, 3] + params_real = (3, 1, 5) + x_approx = np.array(generate_array(bounds[0], bounds[1], 4)) + y_approx = np.array(gen_y(x_approx, *params_real, + noise=0.3, n_outliers=4)) + + params_lsq = sopt.least_squares( + loss, params0, loss='linear', args=(x_approx, y_approx)).x + params_soft_l1 = sopt.least_squares( + loss, params0, loss='soft_l1', args=(x_approx, y_approx),f_scale=0.1).x + params_cauchy = sopt.least_squares( + loss, params0, loss='cauchy', args=(x_approx, y_approx), f_scale=2).x + + x_real = np.array(generate_array(bounds[0], bounds[1])) + y_real = np.array(gen_y(x_real, *params_real, 0, 0)) + + sp = create_subplot() + sp.plot(x_real, y_real, "-b", + label=f"y={exponential.str} на [{bounds[0]};{bounds[1]}]") + sp.plot(x_approx, y_approx, ".r", label=f"Табличные значения с шумом") + sp.plot(x_real, gen_y(x_real, *params_lsq), color="green", + label=f"loss=\"linear\"", linestyle=(0, (5, 10))) + sp.plot(x_real, gen_y(x_real, *params_soft_l1), color="magenta", + label=f"loss=\"soft_l1\"", linestyle=(5, (5, 10))) + sp.plot(x_real, gen_y(x_real, *params_cauchy), color="black", + label=f"loss=\"cauchy\"", linestyle=(7, (5, 10))) + + sp.set_title(f"scipy.optimize.least_squares") + sp.legend(loc='lower left') + + @staticmethod + def plot(method: str = "all"): + if method in ["lagrange", "all"]: + Approx.plot_lagrange() + if method in ["spline", "all"]: + Approx.plot_spline() + if method in ["least_squares_curvefit", "all"]: + Approx.plot_least_squares_curvefit() + if method in ["least_squares", "all"]: + Approx.plot_least_squares() + plt.ylabel("y") + plt.xlabel("x") + plt.show() + + +def main(): + NonLinear.plot() + SLE.print() + Approx.plot() + + +if __name__ == "__main__": + main() diff --git a/conclusion.tex b/conclusion.tex new file mode 100644 index 0000000..df9600c --- /dev/null +++ b/conclusion.tex @@ -0,0 +1,15 @@ +\chapter*{Заключение} + +В ходе выполнения работы были выполнены поставленные задачи, цель +достигнута. В ходе исследования возможностей библиотек выяснилось, +что они не содержат реализаций численных методов для +решения систем нелинейных уравнений. + +В целом, в ходе работы с библиотеками стоит отметить положительные +стороны: гибкость интерфейсов, высокое качество решений, большое +количество доступной информации о результатах работы алгоритма. + +Из недостатков интерфейса библиотек можно отметить его неоднородность --- +некоторые численные методы реализованы с помощью классов, другие +с помощью функций, что может незначительно увеличить время на +освоение интерфейса библиотек и ее возможностей пользователями. diff --git a/experimental_research.tex b/experimental_research.tex new file mode 100644 index 0000000..adc8595 --- /dev/null +++ b/experimental_research.tex @@ -0,0 +1,64 @@ +\chapter{Экспериментальное исследование возможностей библиотек} +Для исследования возможностей библиотек была разработана программа +на языке Python. Она позволяет изучить работу наиболее популярных +численных методов: методов решения нелинейных уравнений и +СЛУ, а также аппроксимации функций. + +Ее структура состоит из классов \textbf{NonLinear}, \textbf{SLE}, +\textbf{Approx}. Они состоят только из статических методов, +пользовательских родительских классов не имеют. + +При запуске программа сначала выводит графики решений нелинейных +уравнений в отдельных окнах (рисунок \ref{bisect}, рисунок\ref{newton}), +входные и выходные данные для методов решения СЛУ в терминале, и +затем результаты аппроксимации, так же в виде графиков в отдельных +окнах. + +Результат работы метода Гаусса (вывод терминала) приведен на +рисунке \ref{gauss}. + +Вывод графиков осуществляется с помощью библиотеки \textbf{matplotlib}, +через функции \textbf{matplotlib.pyplot.plot} и \textbf{matplotlib.pyplot.subplots} \cite{links:matplotlib}. + +Все графики имеют заголовок, в котором написано название функции, +легенду в нижнем левом углу, в которой описаны данные графика. +\begin{figure}[ht] + \centering + \includegraphics[width=0.6\textwidth]{assets/bisect.png} + \caption{Результат исследования функции bisect} + \label{bisect} +\end{figure} + +\begin{figure}[ht] + \centering + \includegraphics[width=0.6\textwidth]{assets/newton.png} + \caption{Результат исследования функции newton} + \label{newton} +\end{figure} + + +\begin{figure}[ht] + \centering + \includegraphics[width=0.6\textwidth]{assets/Gauss.png} + \caption{Результат решения СЛУ методом Гаусса} + \label{gauss} +\end{figure} + +Для получения результатов исследования отдельных классов численных +методов есть следующие методы: +\begin{enumerate} + \item \textbf{NonLinear.plot} для вывода графиков результатов + решения нелинейных уравнений, \textbf{Approx.plot} --- для вывода + графиков решения задачи аппроксимации. + + \item \textbf{SLE.print} --- для вывода результатов решения + СЛУ в терминал +\end{enumerate} + +Данные методы самостоятельно вызываются при запуске программы +пользователем. + +Вывод терминала а также графики для остальных методов приведены +в приложении А. + +Код программы приведен в приложении Б. diff --git a/files/cpp-logo.png b/files/cpp-logo.png deleted file mode 100644 index 432d4eb..0000000 Binary files a/files/cpp-logo.png and /dev/null differ diff --git a/main.tex b/main.tex index ee5ad5e..08cf10b 100644 --- a/main.tex +++ b/main.tex @@ -1,5 +1,10 @@ +% TODO: Fix terminology usage of "якобиан" \input{vars} \input{config} +\sloppy + +\NewDocumentCommand{\MFArgs} + {}{x^{(p)}_1,x^{(p)}_2,x^{(p)}_3,\dots,x^{(p)}_n} \begin{document} \lstset{language=[11]C++} @@ -58,6 +63,7 @@ непрерывна, то обязательно найдется такое \(x_k \in (a,b)\), что \(F(x_k) = 0\) либо \(|F(x_k)| < \varepsilon\), где \(\varepsilon\) --- погрешность искомого решения. + \subsection{Метод деления отрезка пополам} Данный метод использует технику поиска решения, похожую на бинарный поиск. @@ -94,6 +100,7 @@ После применение метода на заданных входных данных получим множество решений уравнения \(Ans = \{x, x \in \textbf{R}\}, |Ans| \geq 1\). + \subsubsection{Реализации метода в библиотеках numpy, scipy} Библиотека scipy содержит функцию \textbf{bisect} из модуля \textbf{scipy.optimize} \cite{links:scipy_doc}, которая реализует @@ -108,10 +115,10 @@ противоположные знаки. \item \(a\) --- scalar - Первый конец интервала \([a,b]\). + Первый конец интервала \([a;b]\). \item \(b\) --- scalar - Второй конец интервала \([a,b]\). + Второй конец интервала \([a;b]\). \item \(xtol\) --- number, необязательный Вычисленный корень \(x0\) будет удовлетворять @@ -147,11 +154,1774 @@ \end{enumerate} \subsection{Метод Ньютона (метод касательных)} +Данный итерационный метод позволит найти единственное решение +уравнения (если оно существует) для каждого из выбранных начальных +приближений. +\subsubsection{Описание метода} +На итерации \(i, i = 1\dots k\), строится касательная к кривой +\(y = F(x)\) в точке \((x_i;F(x_i))\), затем находится \(x_{i+1}\) +--- пересечение данной касательной с осью \(Ox\). Процесс продолжается +пока не будет достигнута заданная точность (\(| F ( x_i ) |< \varepsilon\) +или \(| x_{i+1} - x_i | < \varepsilon\)). + +\subsubsection{Реализации метода в библиотеках numpy, scipy} +Библиотека scipy содержит функцию \textbf{newton} из модуля +\textbf{scipy.optimize} \cite{links:scipy_doc}, которая реализует +данный метод. + +Функция имеет следующие параметры (задаются в порядке перечисления): +\begin{enumerate} + \item \(f\) --- function + + Функция, корень которой требуется. Это должна быть функция + одной переменной вида \(f(x,a,b,c \dots)\), где \(a,b,c \dots\) + --- дополнительные аргументы, которые можно передать в параметре + \(args\). + \item \(x0\) --- float, sequence, или ndarray + + Начальная оценка корня, которая должна быть где-то рядом с + фактическим корнем. Если \(f\) не скалярная, то \(f\) должна + быть векторизована и возвращать последовательность или массив + той же формы, что и ее первый аргумент. + \item \(fprime\) --- callable, необязательный + + Производная функции, когда она доступна и удобна. Если это + \verb|None| (по умолчанию), то используется метод секущей. + \item \(args\) --- tuple, необязательный + + Дополнительные аргументы для использования при вызове функции. + \item \(tol\) --- float, необязательный + + Допустимая погрешность значения корня. Если + \(y=f(x),y \in \textbf{Z}\), рекомендуется большее значение + \(tol\), так как и действительная, и мнимая части \(x\) + вносят вклад в \(|x - x0|\). + \item \(maxiter\) --- int, необязательный + + Максимальное количество итераций. + \item \(fprime2\) --- callable, необязательный + + Производная функции второго порядка, если она доступна и удобна. + Если это \verb|None| (по умолчанию), то используется обычный + метод Ньютона или метод секущих. Если не \verb|None|, то + используется метод Галлея. + \item \(x1\) float, необязательный + + Еще одна оценка корня, которая должна быть где-то рядом с фактическим корнем. Используется, если \(fprime\) не указан. + \item \(rtol\) --- float, необязательный + + Допустимое отклонение (относительное) для прерывания работы. + \item \(full\_output\) --- bool, необязательный + + Если \(full\_output\) имеет значение \verb|False| (по умолчанию), + возвращается корень. Если \verb|True| и \(x0\) --- скаляр, + возвращаемое значение равно \verb|(x, r)|, где \(x\) --- это + корень, а \(r\) --- объект \verb|RootResults|. Если \verb|True| + и \(x0\) --- не скаляр, возвращаемое значение равно + \verb|(x, converged, zero_der)|, где: + \begin{itemize} + \item converged --- \verb|ndarray| из значений bool. Указывает, какие элементы сошлись успешно. + \item zero\_der --- \verb|ndarray| из значений bool. Указывает, какие элементы имеют нулевую производную. + \end{itemize} + \item \(disp\) --- bool, необязательный + + Если \verb|True| и алгоритм не сошелся, будет сгенерировано + исключение \verb|RuntimeError|, с сообщением, содержащим + количество итераций и текущее значение функции. В противном + случае статус сходимости записывается в возвращаемый объект + \verb|RootResults|. Игнорируется, если \verb|x0| не является + скалярным. Примечание: это не имеет ничего общего с + отображением, однако ключевое слово \verb|disp| нельзя + переименовать для сохранения обратной совместимости. +\end{enumerate} + \subsection{Метод простой итерации} -\chapter{Экспериментальное исследование возможностей библиотек} +Данный метод, как и предыдущий, не позволяет найти несколько решений +уравнения за исполнение алгоритма на единственном начальном приближении. +\subsubsection{Описание метода} +Уравнение \(F(x)=0\) приводим к виду \(x = \varphi(x)\), например +\(x = x-F(x)/M\), где \(M\) --- константа. + +Условие сходимости алгоритма: \(0<|\varphi'(x)|<1\). Исходя из него, +\(M\) определяется как +\begin{equation} + M=1.01 \cdot F'(x_0) + \label{formula:stab_simpleiter} +\end{equation} +где \(x_0\) --- +начальное приближение. + +Таким образом, для итерации \(i, i = 1\dots k,\) +\(x_i = \varphi(x_{i-1})\). + +Процесс продолжается, пока не будет достигнута заданная точность, в +данном случае достаточно будет выполнения условия: +\(|x_{i-1}-x_i|<\varepsilon\). +\subsubsection{Реализации метода в библиотеках numpy, scipy} +В библиотеках numpy, scipy не найдено реализаций данного метода. + +\section{Методы решения систем линейных алгебраических уравнений} +Система линейных алгебраических (далее, СЛУ) уравнений имеет вид +\begin{eqnarray} + \left\{ + \begin{aligned} + & a_{11}x_1 + a_{12}x_2 + a_{1n}x_n = b_1 \\ + & a_{21}x_1 + a_{22}x_2 + a_{2n}x_n = b_2 \\ + & \dots\dots\dots\dots\dots\dots\dots\dots \\ + & a_{n1}x_1 + a_{n2}x_2 + a_{nn}x_n = b_n \\ + \end{aligned} + \right. + \label{formula:SLE} +\end{eqnarray} +СЛУ также представима в матричной форме \(AX=B\), где \(A,X,B\) имеют +следующий вид: +\begin{eqnarray} + A = \left( + \begin{aligned} + & a_{11} \quad a_{12} \quad \dots \quad a_{1n} \\ + & a_{21} \quad a_{22} \quad \dots \quad a_{2n} \\ + & \ \vdots \ \qquad \vdots \quad \ \ddots \quad \ \vdots \\ + & a_{n1} \quad a_{n2} \quad \dots \quad a_{nn} \\ + \end{aligned} + \right), B = \left( + \begin{aligned} + & b_1 \\ + & b_2 \\ + & \ \vdots \\ + & b_n \\ + \end{aligned} + \right), X = \left( + \begin{aligned} + & x_1 \\ + & x_2 \\ + & \ \vdots \\ + & x_n \\ + \end{aligned} + \right) + \label{formula:eqn_matrix_system} +\end{eqnarray} + +Для решения таких систем существуют прямые и итерационные методы. +Прямые методы (к ним относятся "метод Гаусса", "метод обратной матрицы" +и "метод прогонки") позволяют получить решение за конечное количество +операций, точность которого ограничивается лишь погрешностью округления. +Итерационные методы (среди которых есть методы, такие как +"метод простой итерации" и "метод Зейделя") позволяют получить +приближенное решение с помощью последовательного приближения к точному. + +Для итерационных методов необходимо начальное приближение, которое +будет обозначено как \(x^{(0)}\). +\subsection{Метод Гаусса} +\subsubsection{Описание метода} +Для решения СЛУ система (\ref{formula:SLE}) приводится к +треугольному виду (\ref{formula:triag_eqn_system}) с помощью цепочки элементарных преобразований. +\begin{eqnarray} + \left\{ + \begin{aligned} + & a'_{11}x_1 + a'_{12}x_2 + a'_{1n}x_n = b'_1 \\ + & 0x_1 + a'_{22}x_2 + a'_{2n}x_n = b'_2 \\ + & \dots\dots\dots\dots\dots\dots\dots. \\ + & 0x_1 + 0x_2 + a'_{nn}x_n = b'_n \\ + \end{aligned} + \right. + \label{formula:triag_eqn_system} +\end{eqnarray} +Данный процесс называется прямым ходом, а нахождение неизвестных +\(x_n, x_{n-1}, \dots,x_1 \) --- обратным. +\subsubsection{Реализации метода в библиотеках numpy, scipy} +В библиотеке scipy реализован частный случай метода Гаусса --- +LU-разложение \cite[с. 259]{book:levitin}. Для получения решения +СЛУ необходимо задействовать две функции из модуля \textbf{scipy.linalg} \cite{links:scipy_doc}: +\begin{enumerate} + \item Для получения разложения используется функция + \textbf{lu\_factor}. + \item Для совершения обратного хода алгоритма используется + \textbf{lu\_solve}, которая принимает на вход разложение с + предыдущего этапа. +\end{enumerate} + +Функция \textbf{lu\_factor} имеет следующие параметры (задаются в порядке перечисления): +\begin{enumerate} + \item \(a\) --- (M, N) array\_like + + Матрица для разложения + \item \(overwrite\_a\) bool, необязательный + + Следует ли перезаписывать данные в A (может повысить + производительность). По умолчанию \verb|False|. + \item \(check\_finite\) bool, необязательный + + Проверять, содержит ли входная матрица только конечные числа. + Отключение может дать прирост производительности, но может + привести к проблемам (сбоям, незавершению), если входные данные + содержат бесконечности или NaN. По умолчанию \verb|True|. +\end{enumerate} + +Функция \textbf{lu\_solve} имеет следующие параметры (задаются в порядке перечисления): +\begin{enumerate} + \item \((lu, piv)\) --- tuple + + Факторизация матрицы коэффициентов a, полученная из + \textbf{lu\_factor}. + \item \(b\) --- array + + Правая сторона + \item \(trans\) --- {0, 1, 2}, необязательный + + Тип системы, которую необходимо решить: + + \begin{tabularx}{0.8\textwidth}{|X|X|} + \hline \(trans\) & вид системы \\ + \hline 0 & \(ax = b\) \\ + \hline 1 & \(a^T x = b\) \\ + \hline 2 & \(a^H x = b\) \\ + \hline + \end{tabularx}\\ + + По умолчанию \verb|0|. + \item \(overwrite\_b\) --- bool, необязательный + + Следует ли перезаписывать данные в \(b\) (может повысить производительность). По умолчанию \verb|False|. + \item \(check\_finite\) --- bool, необязательный + + Проверять, содержат ли входные матрицы только конечные числа. + Отключение может дать прирост производительности, но может + привести к проблемам (сбоям, незавершению), если входные данные + содержат бесконечности или NaN. По умолчанию \verb|True|. +\end{enumerate} + +\subsection{Метод обратной матрицы} +\subsubsection{Описание метода} +Исходная система (\ref{formula:SLE}) представляется в форме +\(AX=B\), тогда вектор неизвестных переменных \(X\) определяется по +формуле (\ref{formula:inv_m_method}). +\begin{equation} + X=A^{-1}B + \label{formula:inv_m_method} +\end{equation} +\subsubsection{Реализации метода в библиотеках numpy, scipy} +Отдельной функции для решения СЛУ не существует, вектор \(X\) можно +найти по формуле (\ref{formula:inv_m_method}). Для получения \(A^{-1}\) +существует функция \textbf{inv} в модуле \textbf{scipy.linalg} +\cite{links:scipy_doc} библиотеки scipy, и функция \textbf{inv} в модуле +\textbf{numpy.linalg} библиотеки numpy. Для перемножения \(A^{-1}\) +и \(B\) в языке Python есть оператор \verb|@|, начиная с версии \(3.5\) +\cite{links:numpy_doc}\cite{links:PEP465}. + +В результате для получения решения СЛУ необходимо выполнить выражение +\verb|inv(A) @ B|, используя одну из вышеописанных функций. + +Функция \textbf{inv} модуля \textbf{scipy.linalg} имеет следующие +параметры (задаются в порядке перечисления): +\begin{enumerate} + \item \(a\) --- array\_like + + Квадратная матрица, которую необходимо инвертировать. + \item \(overwrite\_a\) --- bool, необязательный + + Не запоминать состояние \(a\) (может улучшить + производительность). По умолчанию \verb|False|. + \item \(check\_finite\) --- bool, необязательный + + Проверять, содержат ли входные матрицы только конечные числа. + Отключение может дать прирост производительности, но может + привести к проблемам (сбоям, незавершению), если входные данные + содержат бесконечности или NaN. По умолчанию \verb|True|. +\end{enumerate} +\pagebreak + +Функция \textbf{inv} модуля \textbf{scipy.linalg} имеет следующие +параметры (задаются в порядке перечисления): +\begin{enumerate} + \item \(a\) --- аналогичен параметру \(a\) функции из модуля \textbf{scipy.linalg}. +\end{enumerate} + +\subsection{Метод прогонки} +Данный метод применяется для решения трехдиагональных матриц \hspace{1mm} вида +\begin{equation} + \left\{ + \begin{aligned} + & a_1 x_0 + b_1 x_1 + c_1 x_2 = d_1 \\ + & a_2 x_1 + b_2 x_2 + c_2 x_3 = d_2 \\ + & \dots\dots\dots\dots\dots\dots\dots \\ + & a_n x_{n-1} + b_n x_n + c_n x_{n+1} = d_n \\ + \end{aligned} + \right. + \label{formula:eqn_banded_system} +\end{equation} + +Является частным случаем метода Гаусса. +\subsubsection{Описание метода} +После исключения переменных ниже главной диагонали с помощью +элементарных преобразований, в каждом уравнении СЛУ остается +\(\leq 2\) неизвестных. В таком случае, формула обратного хода будет +следующей: \(x_i = U_i x_{i+1}+V_i, i = n,n-1,\dots,1\). После замены +\(i\) на \(i-1\) и подстановке выражения в общий вид уравнения из СЛУ +(\ref{formula:eqn_banded_system}) \(a_i x_{i-1} + b_i x_i + c_i x_{i+1}\) +получим следующее выражение: +\begin{equation} + x_i = - \frac{c_i}{a_i U_{i-1} + b_i} x_{i+1} + + \frac{d_i - a_i V_{i-1}}{a_i U_{i-1} + b_i} +\end{equation} +Из которого получим: +\begin{equation} + U_i = - \frac{c_i}{a_i U_{i-1} + b_i}, + V_i = \frac{d_i - a_i V_{i-1}}{a_i U_{i-1} + b_i} + \hspace{1cm} i = 1,2,3,\dots,n +\end{equation} +При этом \( c_n = 0; a_1 = 0\), в ходе выполнения алгоритма сначала +вычисляем \(U_i, V_i\), затем \(x_i, i =n,n-1,\dots,1\). +\pagebreak + +Данный метод в общем случае не устойчив, за исключением случаев, +когда матрица СЛУ обладает свойством диагонального преобладания +(условие \ref{formula:eqn_diag_dominant}) или она положительно +определенная \cite{links:bhatia}. +\begin{equation} + \sum_{i \ne j} |a_{ij}| < |a_{ii}|; \qquad i=1,2,3,\dots,n + \label{formula:eqn_diag_dominant} +\end{equation} + +\subsubsection{Реализации метода в библиотеках numpy, scipy} +В scipy для решения СЛУ вида (\ref{formula:eqn_banded_system}) +существует две функции в модуле \textbf{scipy.linalg}: +\textbf{solve\_banded} \cite{links:scipy_doc} и +\textbf{solveh\_banded} \cite{links:scipy_doc}, обе из которых могут +решать задачи для n-диагональных матриц. + +Различие между ними заключается в том, что \textbf{solve\_banded} +не использует метод прогонки, из-за низкой устойчивости метода в общем +случае, что позволяет найти решение даже если матрица не положительно +определенная или не обладает свойством диагонального преобладания. + +\textbf{solveh\_banded} реализует метод прогонки, но авторы библиотеки +указывают, что вводимая матрица должна быть положительно определенной +\cite{links:scipy_doc}. + +Функция \textbf{solve\_banded} имеет следующие параметры (задаются в +порядке перечисления): +\begin{enumerate} + \item \((l, u)\) --- (integer, integer) tuple + + Количество ненулевых нижних и верхних диагоналей. + \item \(ab\) --- (l + u + 1, M) array\_like + + Ленточная матрица. + \item \(b\) --- (M,) or (M, K) array\_like + + Правая сторона. + \item \(overwrite\_ab\) --- bool, необязательный + + Разрешить изменять данные в \(ab\) (может повысить + производительность). По умолчанию \verb|False|. + \item \(overwrite\_b\) --- bool, необязательный + + Разрешить изменять данные в \(b\) (может повысить + производительность). По умолчанию \verb|False|. + \item \(check\_finite\) --- bool, необязательный + + Проверять, содержат ли входные матрицы только конечные числа. + Отключение может дать прирост производительности, но может + привести к проблемам (сбоям, незавершению), если входные данные + содержат бесконечности или NaN. По умолчанию \verb|True|. +\end{enumerate} + +Функция \textbf{solveh\_banded} имеет несколько иной набор параметров +(задаются в порядке перечисления): +\begin{enumerate} + \item \(ab\) --- (u + 1, M) array\_like + + Ленточная матрица, \(u\) --- число верхних диагоналей. + \item \(b\) --- (M,) or (M, K) array\_like + + Правая сторона. + \item \(overwrite\_ab\) --- bool, необязательный + + Разрешить изменять данные в \(ab\) (может повысить + производительность). По умолчанию \verb|False|. + \item \(overwrite\_b\) --- bool, необязательный + + Разрешить изменять данные в \(b\) (может повысить + производительность). По умолчанию \verb|False|. + \item \(lower\) --- bool, необязательный + + Является ли матрица в нижней форме. (По умолчанию используется + верхняя форма), то есть \verb|False|. + \item \(check\_finite\) --- bool, необязательный + + Совпадает с параметром \(check\_finite\) функции + \textbf{solve\_banded}. +\end{enumerate} + +Обе функции принимают матрицу \(ab\) либо в верхней (\textbf{solve\_banded}), либо в нижней форме (\textbf{solveh\_banded} при +включенной опции \(lower\)). Например, для матрицы + +\begin{equation*} + \left( + \begin{aligned} + 5 && 2 && -1 && 0 && 0 \\ + 1 && 4 && 2 && -1 && 0 \\ + 0 && 1 && 3 && 2 && -1 \\ + 0 && 0 && 1 && 2 && 2 \\ + 0 && 0 && 0 && 1 && 1 + \end{aligned} + \right) +\end{equation*} +верхняя форма будет следующей: + +\begin{equation*} + \left( + \begin{aligned} + 0 && 0 && -1 && -1 && -1 \\ + 0 && 2 && 2 && 2 && 2 \\ + 5 && 4 && 3 && 2 && 1 \\ + 1 && 1 && 1 && 1 && 0 +\end{aligned} +\right) +\end{equation*} + +Так как данная матрица не эрмитова, и, следовательно, не положительно +определенна, нижняя форма для нее не существует. Если взять эрмитову +положительно определенную матрицу + +\begin{equation*} + \left( + \begin{aligned} + 4 && 2 && -1 && 0 && 0 && 0 \\ + 2 && 5 && 2 && -1 && 0 && 0 \\ + -1 && 2 && 6 && 2 && -1 && 0 \\ + 0 && -1 && 2 && 7 && 2 && -1 \\ + 0 && 0 && -1 && 2 && 8 && 2 \\ + 0 && 0 && 0 && -1 && 2 && 9 +\end{aligned} +\right) +\end{equation*} +то ее нижняя форма будет следующая: + +\begin{equation*} + \left( + \begin{aligned} + 4 && 5 && 6 && 7 && 8 && 9 \\ + 2 && 2 && 2 && 2 && 2 && 0 \\ + -1 && -1 && -1 && -1 && 0 && 0 +\end{aligned} +\right) +\end{equation*} + +\subsection{Метод простой итерации (метод Якоби)} +\subsubsection{Описание метода} +Для матрицы СЛУ (\ref{formula:eqn_matrix_system}) размеров +\(n \times n\), и начального приближения \(x^{(0)}\) приближенное +решение на итерации \(p, p = 1,2,3,\dots,k\) вычисляется по +следующей формуле: +\begin{equation*} + x^{(p+1)}_i = \frac{1}{a_{ii}} + (b_i - \sum_{j=1}^{i-1}a_{ij} x^{(p)}_j + - \sum_{j=i+1}^{n}a_{ij} x^{(p)}_j); + \quad i = 1,2,3,\dots,n +\end{equation*} + +Метод сходится при \(p \to \infty\), если матрица СЛУ обладает свойством +диагонального преобладания (выполняется условие +\ref{formula:eqn_diag_dominant}). Заданная точность достигается +при выполнении условия: +\begin{equation} + \max_{i \le i \le n} |x^{(p+1)}_i-x^{(p)}_i| < \varepsilon + \label{formula:precision_iter_sle} +\end{equation} +\subsubsection{Реализации метода в библиотеках numpy, scipy} +Реализаций данного метода в библиотеках numpy, scipy не найдено. + +\subsection{Метод Зейделя} +\subsubsection{Описание метода} +Для матрицы СЛУ (\ref{formula:eqn_matrix_system}) размеров +\(n \times n\), и начального приближения \(x^{(0)}\) приближенное +решение на итерации \(p, p = 1,2,3,\dots,k\) вычисляется по +следующей формуле: +\begin{equation*} + x^{(p+1)}_i = \frac{1}{a_{ii}} + (b_i - \sum_{j=1}^{i-1}a_{ij} x^{(p+1)}_j + - \sum_{j=i+1}^{n}a_{ij} x^{(p)}_j); + \quad i = 1,2,3,\dots,n +\end{equation*} + +Данный метод, в отличие от предыдущего, использует уже найденные +компоненты этой же итерации с м\'eньшим индексом. + +Сходимость и точность задаются условиями +(\ref{formula:eqn_diag_dominant}) и (\ref{formula:precision_iter_sle}). +\subsubsection{Реализации метода в библиотеках numpy, scipy} +Реализаций данного алгоритма в библиотеках numpy, scipy не найдено. + +\section{Численные методы решения систем нелинейных уравнений} +Система нелинейных уравнений имеет следующий вид: +\begin{equation} + \begin{aligned} + & F_1(x_1,x_2,\dots,x_n) = 0 \\ + & F_2(x_1,x_2,\dots,x_n) = 0 \\ + & \dots\dots\dots\dots\dots\dots \\ + & F_n(x_1,x_2,\dots,x_n) = 0 \\ + \end{aligned} + \label{formula:SnLE} +\end{equation} + +Для решения данной системы далее будет описано несколько итерационных +методов. Для каждого их них требуется начальное приближение +\(X^{(0)} = \{x^{(0)}_1,x^{(0)}_2,\dots x^{(0)}_n\}\). + +\subsection{Метод простой итерации (метод Якоби) для систем + нелинейных уравнений} +\subsubsection{Описание метода} +Для каждого из уравнений системы (\ref{formula:SnLE}) составляется +уравнение \(f_i: x_i = x_i - F_i(x)/M_i\), где \(M_i\) определяется из +условия сходимости уравнения (\ref{formula:stab_simpleiter}). +В результате получим систему +\begin{equation} + \begin{aligned} + & x_1 = f_1(x_1,x_2,\dots,x_n) \\ + & x_2 = f_2(x_1,x_2,\dots,x_n) \\ + & \dots\dots\dots\dots\dots\dots \\ + & x_n = f_n(x_1,x_2,\dots,x_n) \\ + \end{aligned} + \label{formula:SnLE-Jacobi-fsys} +\end{equation} + +Для получения приближенного решения уравнения применяется формула +\begin{equation} + x^{(p+1)}_{i} = f_i(x^{(p)}_1,x^{(p)}_2,\dots,x^{(p)}_n); + \quad i=1,2,3,\dots,n +\end{equation} +где \(p, p = 1,2,3,\dots,k\) --- номер итерации. + +Заданная точность \(\varepsilon\) достигается выполнением +следующего условия: +\begin{equation*} + \max_i |x^{(p+1)}_i - x^{(p)}_i | < \varepsilon; \quad + i = 1,2,3,\dots,n +\end{equation*} +\subsubsection{Реализации метода в библиотеках numpy, scipy} +Реализаций данного метода в библиотеках numpy, scipy не найдено. + +\subsection{Метод Зейделя для систем нелинейных уравнений} +Данный метод является модификацией предыдущего; отличие состоит в +условии сходимости и формуле получения приближенного решения. + +Ниже будут описаны только вышеперечисленные различия. +\subsubsection{Описание метода} +Формула вычисления приближенного решения данного алгоритма следующая: +\begin{equation*} + \begin{aligned} + & x^{(p+1)}_{1} = f_1(x^{(p)}_1,x^{(p)}_2,\dots,x^{(p)}_n) \\ + & x^{(p+1)}_{2} = f_2(x^{(p)}_1,x^{(p)}_2,\dots,x^{(p)}_n) \\ + & \dots\dots\dots\dots\dots\dots\dots\dots\dots \\ + & x^{(p+1)}_{n} = f_n(x^{(p)}_1,x^{(p)}_2,\dots,x^{(p)}_n) \\ + \end{aligned} +\end{equation*} +где \(p, p = 1,2,3,\dots,k\) --- номер итерации. + +Данный алгоритм более требователен к точности начального приближения. + +Сходимость метода зависит от характера функций исходной системы, +для определения которого необходимо вычислить значения матрицы +\begin{equation} + F' = \left( + \begin{aligned} + & f^{'}_{11} \ \ f^{'}_{12} \ \ f^{'}_{13} \ \ \dots \ \ f^{'}_{1n} \\ + & f^{'}_{21} \ \ f^{'}_{22} \ \ f^{'}_{23} \ \ \dots \ \ f^{'}_{2n} \\ + & \dots \ \ \dots \ \ \dots \ \ \dots \ \ \dots \\ + & f^{'}_{n1} \ \ f^{'}_{n2} \ \ f^{'}_{n3} \ \ \dots \ \ f^{'}_{nn} \\ + \end{aligned} + \right) +\end{equation} +где \(f^{'}_{ij} = \frac{\partial f_i}{\partial x_j}\). + +Сходимость метода обеспечивается выполнением следующего условия: +\vspace{-5mm} +\begin{equation*} + |f^{'}_{i1}| + |f^{'}_{i2}| + |f^{'}_{i3}| + \dots |f^{'}_{in}| < 1; + \quad i = 1,2,3,\dots,n +\end{equation*} + +\subsubsection{Реализации метода в библиотеках numpy, scipy} +Реализаций данного метода в библиотеках numpy, scipy не обнаружено. + +\subsection{Метод Ньютона решения систем нелинейных уравнений} +\subsubsection{Описание метода} +На каждой итерации \(p, p = 1,2,3,\dots,k\) вычисляются значения +\(\Delta x^{(p)}_1,\) \( \Delta x^{(p)}_2, \Delta x^{(p)}_3, \dots, +\Delta x^{(p)}_n \). + +Для этого исходная система уравнений раскладывается в ряд Тейлора по +\(\Delta x^{(p)}_1, \Delta x^{(p)}_2, \Delta x^{(p)}_3, \dots, +\Delta x^{(p)}_n \). Сохранив линейные по данным значениям части, получим +СЛУ: + +\begin{equation} + \begin{aligned} + & \splitdfrac{\frac{\partial F_1(\MFArgs)}{\partial x_1} \Delta x^{(p)}_1 + + \frac{\partial F_1(\MFArgs)}{\partial x_2} \Delta x^{(p)}_2 + \dots +} + {\frac{\partial F_1(\MFArgs)}{\partial x_n} \Delta x^{(p)}_n = -F_1(\MFArgs)} \\ + & \splitdfrac{\frac{\partial F_2(\MFArgs)}{\partial x_1} \Delta x^{(p)}_1 + + \frac{\partial F_2(\MFArgs)}{\partial x_2} \Delta x^{(p)}_2 + \dots +} + {\frac{\partial F_2(\MFArgs)}{\partial x_n} \Delta x^{(p)}_n = -F_2(\MFArgs)} \\ + & \dots\dots\dots\dots\dots\dots\dots\dots\dots\dots\dots\dots\dots\dots\dots \\ + \end{aligned} + \label{formula:SnLE-Newton-sys} +\end{equation} +\begin{equation*} + \begin{aligned} + & \splitdfrac{\frac{\partial F_n(\MFArgs)}{\partial x_1} \Delta x^{(p)}_1 + + \frac{\partial F_n(\MFArgs)}{\partial x_2} \Delta x^{(p)}_2 + \dots +} + {\frac{\partial F_n(\MFArgs)}{\partial x_n} \Delta x^{(p)}_n = -F_n(\MFArgs)} \\ + \end{aligned} +\end{equation*} +Данную систему можно решить любым наиболее подходящим с учетом доступных +ресурсов методом. + +Решением СЛУ (\ref{formula:SnLE-Newton-sys}) будет вектор +\(\Delta X^{(p)} = (\MFArgs)\). После этого, приближенное решение +исходной задачи находится по формуле +\(X^{(p+1)} = X^{(p)} + \Delta X^{(p)}\). + +Заданная точность достигается выполнением условия +(\ref{formula:precision_iter_sle}). Стоит учитывать, что данный метод +так же, как и предыдущий, требователен к точности начального приближения. +\subsubsection{Реализации метода в библиотеках numpy, scipy} +Реализаций данного алгоритма в библиотеках numpy, scipy не найдено, +кроме того, для его работы требуется нахождение частных производных. + +Для поиска частной производной в scipy есть функция \textbf{derivative} +в модуле \textbf{scipy.misc}, но она отмечена устаревшей и будет +удалена в версии 1.12.0, поэтому реализация данного метода с помощью +применения функций общей направленности библиотеки рассматриваться +не будет. + +\section{Аппроксимация функций} +Иногда значения некоторой функциональной зависимости +\(y= \widetilde{y}(x) \) известны для отдельных пар значений +\((x_i,y_i)\). + +Задача восстановления аналитической функции \(\widetilde{y}\) +по отдельным парам значений называется аппроксимацией. +Для получения ее однозначного решения необходимо задать общий вид +функции, включающей коэффициенты, и затем эти коэффициенты определить +с помощью подходящего метода. + +Для определения коэффициентов существует два основных подхода --- +интерполяция (когда \(\widetilde{y}(x_i) = y_i\)), и сглаживание, когда +требуется лишь минимизировать отклонение от известных значений. + +Первые три метода решают поставленную задачу с помощью интерполяции, +последний --- с помощью сглаживания. +\subsection{Интерполяционный полином Лагранжа} +\subsubsection{Описание метода} +Для известных пар значений \((x_i,y_i), x_i \in [a;b]; i = 1,2,3,\dots,n\) +строится полином \(P(x)\), удовлетворяющий условию \(P(x_i) = y_i\). + +Полином \(P(x)\) определяется следующей формулой: +\begin{equation*} + P(x) = L_{n-1}(x) = \sum_{i=1}^{n}y_i\cdot l_i(x) +\end{equation*} +\(l_i(x)\) --- фундаментальные полиномы Лагранжа, которые удовлетворяют равенству (\ref{formula:ip-lagr-eq1}), и имеют следующий вид: +\begin{equation*} + l_i(x) = \frac{(x-x_1)\dots (x-x_{i-1})(x-x_{i+1})\dots(x-x_n)}{(x_i-x_1)\dots(x_i-x_{i-1})(x_i-x_{i+1})\dots(x_i-x_n)} +\end{equation*} +\begin{equation} + l_k(x_i) = \left\{ + \begin{aligned} + & 1, i = k \\ + & 0, i \ne k. \\ + \end{aligned} + \right. + \label{formula:ip-lagr-eq1} +\end{equation} + +Из формулы полинома видно, что его степень не превышает \(n-1\). +Данное свойство сохранится и для следующего метода. + +\subsubsection{Реализации метода в библиотеках numpy, scipy} +Данный метод реализован в библиотеке scipy в виде функции +\textbf{lagrange} в модуле \textbf{scipy.interpolate}. +Одной из ее особенностей является низкая устойчивость --- авторы не +рекомендуют запускать алгоритм на более чем 20 парах \((x,y)\) входных +значений \cite{links:scipy_doc}. + +Данная функция имеет следующие параметры: +\begin{enumerate} + \item \(x\) --- array\_like + + \(x\) представляет координаты \(x\) набора точек данных. + \item \(w\) --- array\_like + + \(w\) представляет координаты \(y\) набора точек данных, + т. е. \(f(x)\). +\end{enumerate} + +Данная функция возвращает полином (тип --- \verb|poly1d| из +библиотеки \textbf{numpy} \cite{links:numpy_doc}). + +\subsection{Интерполяционный полином Ньютона} +\subsubsection{Описание метода} +Интерполяционный полином Ньютона имеет вид +\begin{align*} + N_{n-1}(x) = \Delta^{0}(x_{1}) + \Delta^{1}(x_{1},x_{2})(x-x_{1}) + \Delta^{2}(x_{1},x_{2},x_{3})(x-x_{1})(x-x_{2}) + \\ + + \dots + \Delta^{n-1}(x_{1},x_{2},\dots,x_{n-1})(x-x_{1})(x-x_{2})\dots(x-x_{n-1}) +\end{align*} +где \(\Delta^{0}(x_{1})\dots\Delta^{n-1}(x_{1},x_{2},\dots,x_{n-1})\) --- +конечные разницы порядка \(0,\dots, n-1\), которые вычисляются следующим +образом: +\begin{nospaceflalign*} + & \Delta^{0}(x_i) = y_i & \\ + & \Delta^{1}(x_i,x_k) = \frac{\Delta^{0}(x_i) - \Delta^{0}(x_k)}{x_i-x_k} & \\ + & \Delta^{2}(x_i,x_j,x_k) = \frac{\Delta^{1}(x_i,x_j) - \Delta^{1}(x_j,x_k)}{x_i-x_k} +\end{nospaceflalign*} +и т.д. +\subsubsection{Реализации метода в библиотеках numpy, scipy} +Реализаций данного метода в в библиотеках numpy, scipy не найдено. + +\subsection{Сплайн-интерполяция} +Перед рассмотрением метода необходимо ввести понятие интерполяционного +сплайна. Это кусочно-заданная функция, каждый фрагмент которой является +непрерывной функцией, и на концах его значение совпадает с известными +значениями функции. + +Иногда выбирают фрагменты сплайна такого вида, что непрерывными будут и их +производные, вплоть до нужного порядка (например, если в роли сплайнов +выбрать полиномы 3 степени, то они будут соответствовать условиями +непрерывности вместе с 1 и 2 производными). В зависимости от цели, +на сплайны могут быть наложены и иные ограничения. Также, для упрощения +вычислений, чаще всего все фрагменты сплайна --- функции из одного +класса, например, многочленов фиксированной степени. +\subsubsection{Описание метода} +Между парами значений \((x_{l(j-1)},y_{l(j-1)}),(x_{l(j)},y_{l(j)})\), +где \(l \in \textbf{N}, 0 \le l_j \le n\) строятся +фрагменты сплайна \(f_1(x),\dots,f_n(x)\), соответствущие условиям +непрерывности. Комбинация таких функций образует исходную: +\begin{equation*} + \widetilde{y}(x) = \left\{ + \begin{aligned} + & f_1(x),\; x_0 \le x < x_l \\ + & f_2(x),\; x_{l} \le x < x_{2l} \\ + & \ \dots \qquad \quad \dots \\ + & f_n(x),\; x_{jl} \le x < x_n \\ + \end{aligned} + \right. +\end{equation*} +При этом +\begin{equation*} + \begin{aligned} + & f_1(x_0) = y_0, f_1(x_l) = y_l \\ + & f_2(x_l) = y_l, f_2(x_{2l}) = y_{2l} \\ + & \ \dots \qquad \quad \dots \\ + & f_n(x_{jl}) = y_{jl},f_n(x_n) = y_n \\ + \end{aligned} +\end{equation*} +Для полиномов 3 степени чаще всего \(l=3\), и в таком случае по \(l\) +точкам на каждом фрагменте строиться полином (например, интерполяционный +полином Ньютона). В общем случае \(l\) определяется исходя из условий +и ограничений поставленной перед исследователем задачи. +\subsubsection{Реализации метода в библиотеках numpy, scipy} +В библиотеке scipy существует множество реализаций данного метода, +которые отличаются используемым видом сплайна а также поддерживаемой +размерностью функции. В модуле +\textbf{scipy.interpolate} описаны следующие объекты +\cite{links:scipy_doc}: +\begin{enumerate} + \item класс \textbf{CubicSpline}; + \item класс \textbf{PchipInterpolator}; + \item класс \textbf{CubicHermiteSpline}; + \item класс \textbf{Akima1DInterpolator}; + \item классы \textbf{RectBivariateSpline}, + \textbf{LSQBivariateSpline}, используются для многомерной + интерполяции; + \item функция \textbf{interp1d}, которая может интерполировать таблично + заданную функцию разными методами, в т.ч. сплайн-интерполяцией. + Признана устаревшей и поэтому не будет рассмотрена; +\end{enumerate} -\chapter*{Заключение} +Классы \textbf{CubicSpline}, \textbf{PchipInterpolator}, +\textbf{CubicHermiteSpline}, \textbf{Akima1DInterpolator} +используются для интерполяции функций одного аргумента, и их +применение одинаково: для создания сплайна необходимо создать +экземпляр класса, для получения значений сплайна необходимо +вызвать созданный экземпляр с необходимыми входными данными +(то есть вызвать метод \verb|__call__|). + +Каждый из классов имеет разное количество параметров, которые +можно задать при его создании его экземпляра; метод получения +данных сплайна, в свою очередь, имеет одинаковое количество +параметров для всех классов, поэтому он имеет смысл описать его сейчас. + +Метод \verb|__call__| вышеописанных классов имеет 3 параметра: +\begin{enumerate} + \item \(x\) --- array\_like + + Точки, для которых будет вычислено значения сплайна. + \item \(nu\) --- int, необязательный + + Порядок производной сплайна, который необходимо использовать при + вычислении значений. Должен быть неотрицательным. По умолчанию + --- 0 (то есть вычисляется исходный сплайн). + \item \(extrapolate\) --- {bool, \verb|"periodic"|, + \verb|None|}, необязательный + + Если bool, определяет, следует ли экстраполировать точки, + которые выходят за пределы границ, установленные первым и + последним интервалами или возвращать \verb|NaN|. Если + \verb|"periodic"|, используется периодическая экстраполяция. + Если \verb|None| (по умолчанию), используйте метод \textbf{extrapolate} + этого класса. +\end{enumerate} +Данный метод возвращает массив значений сплайна, соответствующих значений +\(x\). + +Конструктор экземпляра класса \textbf{CubicSpline} имеет следующие +параметры: +\begin{enumerate} + \item \(x\) --- array\_like, одномерный + + Одномерный массив, содержащий значения независимой переменной. + Значения должны быть действительными, конечными числами и + находиться в строго возрастающем порядке. + \item \(y\) --- array\_like + + Массив, содержащий значения зависимой переменной. Он может + иметь произвольное количество измерений, но длина должна + совпадать с длиной \(x\). Значения должны быть конечными. + \item \(axis\) int, необязательный + + Ось, вдоль которой предполагается изменение y. Это означает, + что для \(x[i]\) соответствующие значения равны + \verb|np.take(y, i, axis=axis)|. По умолчанию --- 0. + \item \(bc\_type\) --- string / tuple размера = 2, необязательный + + Тип граничного условия. Два дополнительных уравнения, заданные + граничными условиями, необходимы для определения всех + коэффициентов многочленов на каждом отрезке. + + Если \(bc\_type\) --- строка, то указанное в параметре условие + будет применено к обоим концам сплайна. Доступные условия: + \begin{itemize} + \item \verb|"not-a-knot"| (по умолчанию): первый и второй + сегменты на конце кривой представляют собой один и тот + же полином. Это хороший вариант по умолчанию, когда нет + информации о граничных условиях. + \item \verb|"periodic"|: Предполагается, что интерполируемые + функции являются периодическими с периодом + \(x[-1] \ \mbox {---}\ x[0]\). + Первое и последнее значение y должны быть идентичными: + \(y[0] == y[-1]\). Это граничное условие приведет к + тому, что \(y'[0] == y'[-1]\) и \(y''[0] == y''[-1]\). + \item \verb|"clamped"|: Первая производная на концах кривых + равна нулю. Предполагая, что y одномерный, + \(bc\_type=((1, 0.0), (1, 0.0))\) --- то же самое + условие. + \item \verb|"natural"|: Вторая производная на концах кривой + равна нулю. Предполагая, что y одномерный, + \(bc\_type=((2, 0.0), (2, 0.0))\) --- то же самое + условие. + \end{itemize} + + Если \(bc\_type\) представляет собой кортеж из двух элементов, + первое и второе значения будут применены в начале и конце + кривой соответственно. Значения кортежа могут быть одной из + ранее упомянутых строк (кроме \verb|"periodic"|) или кортежем + \((order, deriv\_values)\), позволяющим указывать произвольные + производные на концах кривой: + + \begin{itemize} + \item \(order\): порядок производной --- 1 или 2. + \item \(deriv\_value\): array\_like + + Содержит значения производной, форма должна быть такой + же, как \(y\), за исключением измерения \(axis\). + Например, если \(y\) --- одномерный, то + \(deriv\_values\) должен быть скаляром. Если \(y\) + является трехмерным, имеет форму \((n0, n1, n2)\) и + \(axis=1\), то \(deriv\_values\) должно быть + двухмерным и иметь форму \((n0, n2)\). + \end{itemize} + \item \(extrapolate\) --- {bool, \verb|"periodic"|, + \verb|None|}, необязательный + + Если bool, определяет, следует ли экстраполировать точки, + которые выходят за пределы границ, установленные первым и + последним интервалами или возвращать \verb|NaN|. Если + \verb|"periodic"|, используется периодическая экстраполяция. + Если \verb|None| (по умолчанию), параметр имеет значение + \verb|"periodic"| если \(bc\_type\)=\verb|"periodic"| и + значение \verb|True| в противном случае. +\end{enumerate} + +Класс \textbf{PchipInterpolator} представляет из себя кубический +эрмитов сплайн, поэтому является дочерним классом +\textbf{CubicHermiteSpline}. Конструктор экземпляра класса имеет +следующие параметры: +\begin{enumerate} + \item \(x\) --- ndarray, вида (xlen,1), то есть одномерный + + Одномерный массив монотонно возрастающих действительных + значений. Не может включать повторяющиеся значения (иначе + сплайн будет слишком конкретизированный) + \item \(y\) --- ndarray, вида (...,xlen,...) + + N-мерный массив действительных значений. Длина \(y\) по оси + интерполяции должна быть равна длине \(x\). Используйте + параметр \(axis\), чтобы выбрать ось интерполяции. + \item \(axis\) int, необязательный + + См. описание одноименного метода родительского класса. + \item \(extrapolate\) --- bool, необязательный + + Определяет, следует ли экстраполировать точки, + которые выходят за пределы границ, установленные первым и + последним интервалами или возвращать \verb|NaN|. + По умолчанию \verb|None|, значение см. в пункте + \ref{list:CubicHermiteSpline-extrapolate} описания + конструктора родительского класса. +\end{enumerate} + +Конструктор экземпляра класса \textbf{CubicHermiteSpline} имеет следующие +параметры: +\begin{enumerate} + \item \(x\) --- array\_like, вида (n,), то есть одномерный + + Одномерный массив, содержащий значения независимой переменной. Значения должны быть действительными, конечными и находиться в строго возрастающем порядке. + \item \(y\) --- array\_like + + Массив, содержащий значения зависимой переменной. + Он может иметь произвольное количество измерений, но длина по + оси \(axis\) должна совпадать с длиной \(x\). + Значения должны быть конечными. + \item \(dydx\) --- array\_like + + Массив, содержащий производные зависимой переменной. Он может + иметь произвольное количество измерений, но длина по оси + \(axis\) должна совпадать с длиной \(x\). Значения должны быть + конечными. + \item \(axis\) --- int, необязательный + + Ось, вдоль которой предполагается изменение y. Это означает, что для \(x[i]\) соответствующие значения равны + \verb|np.take(y, i, axis=axis)|. + По умолчанию --- 0. + \item \(extrapolate\) --- {bool, \verb|"periodic"|, + \verb|None|}, необязательный + + Если bool, определяет, следует ли экстраполировать точки, + которые выходят за пределы границ, установленные первым и + последним интервалами или возвращать \verb|NaN|. Если + \verb|"periodic"|, используется периодическая экстраполяция. + Если \verb|None| (по умолчанию), используйте метод \textbf{extrapolate} + этого класса. + \label{list:CubicHermiteSpline-extrapolate} +\end{enumerate} + +Конструктор экземпляра класса \textbf{Akima1DInterpolator} имеет +следующие параметры: +\begin{enumerate} + \item \(x\) --- ndarray, вида (n, 1) + + Одномерный массив монотонно возрастающих действительных + значений. + \item \(y\) --- ndarray, вида (..., n, ...) + + N-мерный массив реальных значений. Длина \(y\) по оси + интерполяции должна быть равна длине \(x\). Используйте + параметр \(axis\), чтобы выбрать ось интерполяции. + \item \(axis\) --- int, необязательный + + Ось в массиве \(y\), соответствующая значениям координаты + \(x\). По умолчанию --- 0. +\end{enumerate} + +\subsection{Сглаживание. Метод наименьших квадратов} +С помощью метода наименьших квадратов (далее, МНК) осуществляется +интерполяция функции с помощью сглаживания. +\subsubsection{Описание метода} +Задача метода --- минимизировать отклонение от исходных данных: +\begin{equation} + S = \sum_{i=1}^{n}(\widetilde{y}(x_i,a)-y_i)^2 \rightarrow \min + \label{formula:mnk} +\end{equation} +где \(a\) --- неизвестные параметры; \(n\) --- количество известных +пар значений (\(x,y\)); \(x_i,y_i\) --- значения \(i\) пары. + +В качестве \(\widetilde{y}(x)\) необходимо выбрать функцию вида, +который лучше всего подходит под исходные данные. Этот процесс +называют выбором зависимости. + +Примеры зависимостей (\(a,b,c,d\) --- неизвестные параметры): +\begin{alignat*}{2} + & ax + b & \ & \text{--- линейная}; \\ + & ax^2 + bx + c & \ & \text{--- полиномиальная}; \\ + & a \ln(x) + b & \ & \text{--- логарифмическая}; \\ + & be^{ax} & \ & \text{--- экспоненциальная}; +\end{alignat*} + +В результате подстановки выбранной зависимости и известных данных +в (\ref{formula:mnk}) задача минимизации сводится к нахождению +экстремумов функции \(S(a)\), и выбору подходящего из них. +\subsubsection{Реализации метода в библиотеках numpy, scipy} +В модуле \textbf{scipy.optimize} библиотеки scipy существует +несколько реализаций данного метода, каждая из которых имеет свою +специализацию. +\begin{enumerate} + \item \textbf{curve\_fit} --- функция, предназначенная для поиска + значения параметра \(a\) произвольной зависимости \(\widetilde + {y}(x_i,a)\), при котором выполняется условие + (\ref{formula:mnk}). + \item \textbf{least\_squares} --- функция, которая отличается от + \textbf{curve\_fit} возможностью задать собственный критерий + близости, область допустимых решений; позволяет задать больше + параметров точности. + \item \textbf{leastsq} --- функция, которая является устаревшей + версией \textbf{least\_squares}, поэтому ее подробного + рассмотрения не будет. + \item \textbf{lsq\_linear} --- функция, предназначенная для + минимизации выражения \(||ax-b||^2\), где \(a,b\) --- матрицы + параметров; функция \textbf{nnls} отличается от нее + ограничением на положительность искомых значений. + + Подходят только для нахождения параметров многочленов, + не буду рассматриваться в данной работе по причине + сложной подготовки входных данных для решения задачи + метода, а также их узкой специализации. +\end{enumerate} + +Функция \textbf{curve\_fit} имеет следующий набор параметров: +\begin{enumerate} + \item \(f\) --- callable + + Модельная функция \(f(x,\dots)\). Она должен принимать независимую переменную в качестве первого аргумента, а параметры --- в качестве оставшихся. + \item \(xdata\) --- array\_like + + Значения независимой переменной. Обычно это должна быть последовательность длины \(M\) или массив \((k,M)\)-образной формы для функций с \(k\) предикторами, и каждый элемент должен быть конвертируемым в \(float\), если это объект, подобный массиву. + \item \(ydata\) --- array\_like + + Значения зависимой переменной, массив длиной M. + При найденных параметрах \(f\) совпадает с значением + \(f(xdata, \dots)\). + \item \(p0\) --- array\_like, необязательный + + Начальное значение параметров \(f\) (длины N). Если + \verb|None| (значение по умолчанию), то все начальные + значения будут равны \verb|1| (если количество параметров + функции можно определить с помощью интроспекции, в противном + случае генерируется исключение \verb|ValueError|). + \item \(sigma\) --- \verb|None| / sequence длины M / массив MxM, необязательный + + Задает неопределенность в \(ydata\). Если мы определим + остатки как \(r = ydata - f(xdata, *popt)\), то интерпретация + значения \(sigma\) будет зависеть от ее количества измерений: + \begin{itemize} + \item одномерная \(sigma\) должна содержать значения + стандартных отклонений ошибок в \(ydata\). Используется + при вычислении параметра + \verb|chisq = sum((r/sigma)**2)|; + \item двухмерная \(sigma\) должна содержать ковариационную + матрицу ошибок в \(ydata\). Используется + при вычислении параметра + \verb|chisq = r.T @ inv(sigma) @ r|. + \end{itemize} + \item \(absolute\_sigma\) --- bool, необязательный + + Если значение \verb|True|, \(sigma\) используется в абсолютном смысле, и предполагаемая ковариация параметра \(pcov\) отражает эти абсолютные значения. + + Если значение \verb|False| (по умолчанию), имеют значение + только относительные величины значений \(sigma\). + Ковариационная матрица возвращаемого параметра \(pcov\) + основана на масштабировании \(sigma\) постоянным + коэффициентом. Эта константа устанавливается требованием, чтобы + приведенный выше \(chisq\) для параметров + \(popt\) при использовании масштабированной \(sigma\) был равен + единице. Другими словами, \(sigma\) масштабируется так, чтобы + соответствовать выборочной дисперсии остатков после подгонки. + Математически, + \verb|pcov(absolute_sigma=False) = pcov(absolute_sigma=True)| + \verb| * chisq(popt)/(M-N)|. + \item \(check\_finite\) --- bool, необязательный + + Если \verb|True|, то входные массивы проверяются на содержание + значений \verb|NaN| или \verb|Inf|, и если они содержат, то + генерируется исключение \verb|ValueError|. + + Значение \verb|False| может привести к выдаче бессмысленных + результатов, если входные массивы содержат \verb|NaN|. + + По умолчанию установлено значение \verb|True|, если + \(nan\_policy\) не указано явно, и значение \verb|False| в + противном случае. + \item \(bounds\) --- tuple из двух array\_like / \verb|Bounds|, + необязательный + + Нижняя и верхняя границы параметров \(f\). По умолчанию нет ограничений. Есть два способа указать границы: + \begin{itemize} + \item через экземпляр класса \verb|Bounds|; + + \item \verb|tuple| из двух array\_like, где каждый элемент + кортежа должен быть либо массивом с длиной, равной + количеству параметров, либо скаляром (в этом случае + граница считается одинаковой для всех параметров). + Используйте \verb|np.inf| с соответствующим знаком, + чтобы отключить ограничения для всех или некоторых + параметров. + \end{itemize} + \item \(method\) --- \{\verb|"lm"|, \verb|"trf"|, \verb|"dogbox"|\}, необязательный + + Метод, используемый для вычисления параметров. См. описание + аналогичного параметра функции \textbf{least\_squares} для + получения более подробной информации. По умолчанию используется + \verb|"lm"| для задач без ограничений и \verb|"trf"|, если + указано значение \(bounds\). + \item \(jac\) --- callable / string / \verb|None|, необязательный + + Функция с сигнатурой \(jac(x, ...)\), которая вычисляет матрицу + Якоби модельной функции относительно параметров как плотную + array\_like структуру. Он будет масштабироваться в + соответствии с предоставленным параметром \(sigma\). Если + \verb|False| (по умолчанию), якобиан будет оцениваться + численно. Строковые ключевые слова для методов \verb|"trf"| и + \verb|"dogbox"| можно использовать для выбора + конечно-разностной схемы, см. описание + \textbf{least\_squares}. + \item \(full\_output\) --- boolean, необязательный + + Если значение = \verb|True|, эта функция возвращает дополнительную информацию: \(infodict\), \(mesg\) и \(ier\). По умолчанию + \verb|False|. + \item \(nan\_policy\) --- {\verb|"raise"|, \verb|"omit"|, \verb|None|}, необязательный + + Определяет, как действовать, если входные данные содержат \verb|NaN|. Доступны следующие параметры (по умолчанию --- \verb|None|): + \begin{itemize} + \item \verb|"raise"|: выдает ошибку + + \item \verb|"omit"|: выполняет вычисления, игнорируя + значения \verb|NaN|. + + \item \verb|None|: никакая специальная обработка \verb|NaN| + не выполняется (кроме того, что делается с помощью + \(check\_finite\)); поведение при наличии + \verb|NaN| зависит от реализации и может измениться. + + \end{itemize} + Обратите внимание: если значение указано явно (не \verb|None|), + для \(check\_finite\) будет установлено значение \verb|False|. + \item \(**kwargs\) + + Именованные параметры, которые передаются в \textbf{leastsq} + для \(method\)='lm' или \textbf{least\_squares} в противном + случае. +\end{enumerate} +Данная функция возвращает результат в виде +\((popt,pcov,infodict,mesg,ier)\), где +\begin{enumerate} + \item \(popt\) --- вычисленные параметры + \(f\), при которых выполняется условие + \(f(xdata, *popt) - ydata \rightarrow \min\). + \item \(pcov\) --- 2-D array + + Вычисленная приблизительная ковариация \(popt\). + \item \(infodict\) --- dict + + По ключу \(nfev\) + хранится значение количества вызовов функции. Методы \("trf"\) + и \("dogbox"\) не учитывают вызовы функций для численной + аппроксимации якобиана, в отличие от метода \("lm"\). + По ключам \(fvec\), \(fjac\), \(ipvt\), \(qtf\) можно получить + дополнительную информацию о работе алгоритма. + \item \(mesg\) --- str + + Строковое сообщение с информацией о решении. + \item \(ier\) --- int + + Целочисленный флаг. Если он равен 1, 2, 3 или 4, решение + найдено. В противном случае решение не было найдено. \(mesg\) + предоставляет дополнительную информацию о решении. +\end{enumerate} +Разработчики отмечают, что входные данные \(xdata\), \(ydata\) и +выходные данные \(f\) должны иметь формат \(float64\), иначе данная +функция может вернуть неверные результаты \cite{links:scipy_doc}. + +Функция \textbf{least\_squares} предназначена для решения следующей +задачи минимизации \(F(x)\), где +\verb|F(x) = 0.5 * sum(rho(f_i(x)**2), i = 0, ..., m - 1)| +при ограничениях \(lb \leq x \leq ub\). То есть, \(f\) должна принимать +на вход значение параметров \(a\) и возвращать значение \(f(x,a) - y\). + +Данная функция имеет следующий набор параметров: +\begin{enumerate} + \item \(fun\) --- callable + Функция, которая вычисляет вектор остатков, имеет сигнатуру + \(fun(x, *args, **kwargs)\), т.е. минимизация продолжается + относительно ее первого аргумента. Аргумент \(x\), + передаваемый этой функции, представляет собой N-мерный массив размерности \((N,K)\) (никогда не скаляр, даже для \(N=1\)). \(fun\) должна вернуть одномерный массив типа array \((m,)\) или скаляр. Если \(x \in \textbf{C}\) или функция \(fun\) возвращает комплексные числа, ее необходимо обернуть действительнозначной функцией от действительнозначных аргументов. + \item \(x_0\) --- array\_like размерности (N,K) / float + + Начальное значение \(x\). Если параметр имеет тип float, он + будет рассматриваться как одномерный массив с одним элементом. + \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|. + \item \(bounds\) --- tuple из 2 array\_like / \verb|Bounds|, необязательный + + Есть два способа указать границы: + \begin{itemize} + \item через экземпляр класса \verb|Bounds|; + + \item \verb|tuple| из двух array\_like, где каждый элемент + кортежа должен быть либо массивом с размерностью + \(x0\), либо скаляром (в этом случае + граница считается одинаковой для всех параметров). + Используйте \verb|np.inf| с соответствующим знаком, + чтобы отключить ограничения для всех или некоторых + параметров. + \end{itemize} + \item \(method\) --- \{\verb|"lm"|, \verb|"trf"|, \verb|"dogbox"|\}, необязательный + + Допустимое отклонение на завершение по норме градиента. + По умолчанию --- 1e-8. Точное состояние зависит от используемого + значения \(method\): + \begin{itemize} + \item Для \verb|"trf"|: \(norm(g\_scaled, ord=np.inf) < gtol\), где \(g\_scaled\) --- значение градиента, масштабированное для учета наличия границ. + \item Для \verb|"dogbox"|: \(norm(g_free, ord=np.inf) < gtol\), где \(g_free\) --- градиент по отношению к переменным, которые не находятся в оптимальном состоянии на границе. + \item Для \verb|"lm"|: максимальное абсолютное значение косинуса углов между столбцами якобиана и вектором невязки меньше \(gtol\), или вектор невязки равен нулю. + \end{itemize} + + Если \verb|None| и \(method\) не \verb|"lm"|, завершение по этому условию отключено. Если \(method\) --- \verb|"lm"|, этот допуск должен быть выше, чем машинный эпсилон. + \item \(x\_scale\) --- array\_like / \verb|"jac"| + + Характеристический масштаб каждой переменной. + Установка \(x\_scale\) эквивалентна переформулировке проблемы + в масштабированных переменных \(xs = x/x\_scale\). + Альтернативная точка зрения состоит в том, что размер + доверительной области по j-му измерению пропорционален + \(x\_scale[j]\). Улучшения сходимости можно достичь, установив + \(x\_scale\) таким образом, чтобы шаг заданного размера по + любой из масштабируемых переменных имел аналогичный эффект на + функцию стоимости. Если установлено значение \verb|"jac"|, + масштаб итеративно обновляется с использованием обратных норм + столбцов матрицы Якоби. + \item \(loss\) --- str / callable, необязательный + + Определяет функцию потерь. Допускаются следующие значения ключевых слов: + \begin{itemize} + \item \verb|"linear"| (по умолчанию) : \verb|rho(z) = z|. Дает стандартную задачу наименьших квадратов. + \item \verb|"soft_l1"| : \verb|rho(z) = 2 * ((1 + z)**0.5 - 1)|. Гладкая аппроксимация потерь l1 (абсолютное значение). Обычно хороший выбор для надежного метода наименьших квадратов. + \item \verb|"huber"| : \verb|rho(z) = z if z <= 1 else 2*z**0.5 - 1|. Работает аналогично \verb|"soft_l1"|. + \item \verb|"cauchy"| : \verb|rho(z) = ln(1 + z)|. Сильно ослабляет влияние выбросов, но может вызвать трудности в процессе оптимизации. + \item \verb|"arctan"| : \verb|rho(z) = arctan(z)|. Ограничивает максимальный убыток по одному остатку, имеет свойства, аналогичные \verb|"cauchy"|. + \end{itemize} + + Если параметр имеет тип \verb|callable|, он должен принимать + одномерный \verb|ndarray| \verb|z=f**2| и возвращать + array\_like с формой \((3, m)\), где строка 0 содержит значения функции, строка 1 содержит первые производные, а строка 2 содержит вторые производные. \(method\) \verb|lm| поддерживает только \verb|"linear"| \(loss\). + + \item \(f\_scale\) --- float, необязательный + + Значение границы между ошибками данных и выбросами в данных, + по умолчанию - 1,0. Значение функции потерь оценивается + следующим образом: \verb|rho_(f**2) = C**2 * rho(f**2 / C**2)|, + где C --- это \(f\_scale\), а \(rho\) определяется параметром + \(loss\). Этот параметр не имеет никакого эффекта при + \(loss\)=\verb|"linear"|, но для других значений потерь он имеет решающее + значение. + \item \(max\_nfev\) --- \verb|None| / int, необязательный + + Максимальное количество выполнений функции перед завершением работы. Если \verb|None| (по умолчанию), значение выбирается автоматически: + \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"| учитывает вызовы функций в оценке якобиана). + \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"|}, необязательный + + Метод решения подзадач доверительной области, применим только для методов \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 + + Массив, который необходимо интегрировать. + \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 из callable, необязательный + + События для отслеживания. Если \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} +\input{experimental_research} +\input{conclusion} \addcontentsline{toc}{chapter}{Заключение} \input{sources} diff --git a/presentation.tex b/presentation.tex index 3eb76aa..9b772dc 100644 --- a/presentation.tex +++ b/presentation.tex @@ -1,6 +1,4 @@ \documentclass{beamer} - - \usepackage[russian]{babel} \usepackage[utf8]{inputenc} \usepackage[outputdir=cache]{minted} @@ -14,7 +12,7 @@ \title[]{\cwtitle} \institute[]{ФГБОУ ВО «Вятский государственный университет»} \date{\null} -\author[ ]{Студент ПМИб-3301-52-00 \cwauthor\\ \and к.п.н. А.Н.~Соколова} +\author[ ]{Студент ПМИб-3301-52-00 \cwauthor \newline \and к.п.н. А.Н.~Соколова} \newcommand\frametitleSpec[1]{% \frametitle{#1} @@ -27,43 +25,150 @@ \begin{document} \begin{frame} - \centering\includegraphics[width=0.4\textwidth]{files/vyatsu_logo.png}\\ + \centering\includegraphics[width=0.4\textwidth]{assets/vyatsu_logo.png}\\ \titlepage \end{frame} \begin{frame} \frametitle{План доклада} - \tableofcontents \end{frame} \begin{frame} \frametitleSpec{Введение} - Активное внедрение компьютеров во всевозможные отрасли жизни человека привело к тому, что при решении прикладных задач требование к скорости и дешевизне разработки стало выше требований производительности и ресурсоемкости разрабатываемой программы. + Для языка Python разработаны библиотеки numpy и scipy + математической направленности. + Они включают в себя множество алгоритмов решения разнообразных задач, в том числе для интегрирования функций, осуществления операций над массивами и + матрицами. \begin{enumerate} - \item Проблема состоит в том, что на данный момент ОС, которая бы могла более эффективно использовать текущее АО устройств, не разработана. - \item Целью данной работы является рассмотрение архитектуры и разработка части прототипа данной ОС. + \item Проблема состоит в том, что на данный момент особенности + библиотек, полнота их возможностей с точки зрения решения + задач численными методами недостаточно исследованы. + \item \textbf{Целью} данной работы является исследование + вышеприведенных характеристик данных библиотек. + \item Для достижения цели курсового проекта необходимо выполнить + следующие \textbf{задачи}: \begin{itemize} - \item Изучить архитектуру и требования прототипа ОС. - \item Определить перечень сервисов обеспечения целостности и оптимизации экосистемы устройств. - \item Разработать БД и сервисы обеспечения целостности и оптимизации прототипа экосистемы устройств. + \item Изучить распространенные численные методы решения + основных классов задач. + \item Изучить документацию библиотек numpy и scipy на + предмет реализации рассмотренных методов. + \item Экспериментально исследовать возможности + реализаций рассмотренных численных методов данных + библиотек. \end{itemize} \end{enumerate} \end{frame} \begin{frame} - \frametitleSpec{О программной реализации} - Реализация сервисов написана на языке C++ в виде отдельных компонентов. + \frametitleSpec{Рассмотренные численные методы} + + Всего было рассмотрено 20 численных методов, каждый из которых + решает свой класс задач: + \begin{enumerate} + \item Решение нелинейных уравнений + \begin{itemize} + \item Метод деления отрезка пополам + \item Метод касательных (Ньютона) + \item Метод простой итерации + \end{itemize} + \item Решение СЛУ: + \begin{itemize} + \item Метод Гаусса + \item Метод обратной матрицы + \item Метод прогонки + \item Метод простой итерации + \item Метод Зейделя + \end{itemize} + \item Решение систем нелинейных уравнений, с помощью метода Ньютона, и модифицированных версий методов Зейделя и простой итерации. + \end{enumerate} + +\end{frame} +\begin{frame} + \begin{enumerate} + \item Аппроксимация функций + \begin{itemize} + \item Интерполяционные полиномы Лагранжа, Ньютона + \item Сплайн-интерполяция + \item Сглаживание. Метод наименьших квадратов + \end{itemize} + \item Численное интегрирование + \begin{itemize} + \item Метод трапеций + \item Метод парабол + \end{itemize} + \item Решение задачи Коши ОДУ + \begin{itemize} + \item Метод Эйлера, и его модифицированная версия + \item Метод Рунге-Кутта + \end{itemize} + \end{enumerate} + +\end{frame} +\begin{frame} + \frametitleSpec{Экспериментальное исследование библиотек} + Была разработана программа на языке Python, использующая + возможности scipy и numpy для решения задач. \begin{figure}[h] \centering - \includegraphics[width=0.25\textwidth]{files/cpp-logo.png} + \includegraphics[width=0.25\textwidth]{assets/python-logo.png} \end{figure} + + В программе были приведены примеры решения СЛУ, нелинейных уравнений + и задач аппроксимации функций с помощью описанных в данной работе + методов. Программа не интерактивная. \end{frame} +\begin{frame} + \frametitleSpec{Результаты работы программы} + \begin{figure} + \includegraphics[width=0.5\textwidth]{assets/bisect} + \end{figure} + \begin{figure} + \includegraphics[width=0.65\textwidth]{assets/Gauss} + \end{figure} + +\end{frame} +\begin{frame} + \begin{figure} + \includegraphics[width=1\textwidth]{assets/Thomas.png} + \end{figure} +\end{frame} +\begin{frame} + \begin{figure} + \includegraphics[width=1\textwidth]{assets/lagrange.png} + \end{figure} +\end{frame} +\begin{frame} + \begin{figure} + \includegraphics[width=1\textwidth]{assets/CubicSpline} + \end{figure} +\end{frame} +\begin{frame} + \begin{figure} + \includegraphics[width=1\textwidth]{assets/Akima1DInterpolator.png} + \end{figure} +\end{frame} +\begin{frame} + \begin{figure} + \includegraphics[width=1\textwidth]{assets/curve_fit.png} + \end{figure} +\end{frame} \begin{frame} \frametitleSpec{Заключение} - При реализации прототипа экосистемы был разработан один из возможных вариантов деления на подсистемы, и рассмотрены возможные пути реализации. - Из недостатков текущей реализации можно отметить отсутствие компонента, отвечающего за безопасность и разделения прав доступа пользователей к вычислительным ресурсам сервера. - Таким образом, задачи проекта выполнены, цель проекта достигнута. + В ходе выполнения работы были выполнены поставленные задачи, цель + достигнута. В ходе исследования возможностей библиотек выяснилось, + что они не содержат реализаций численных методов для + решения систем нелинейных уравнений. + + В целом, в ходе работы с библиотеками стоит отметить положительные + стороны: гибкость интерфейсов, высокое качество решений, большое + количество доступной информации о результатах работы алгоритма. + + Из недостатков интерфейса библиотек можно отметить его + неоднородность --- некоторые численные методы реализованы с помощью + классов, другие с помощью функций, что может незначительно + увеличить время на освоение интерфейса библиотек и ее возможностей + пользователями. \end{frame} \begin{frame} \begin{center} diff --git a/sources.tex b/sources.tex index 066532e..a9aeb6f 100644 --- a/sources.tex +++ b/sources.tex @@ -1,23 +1,28 @@ \newcommand{\LiteratureAccessDate}[1][1]{% - дата~обращения: \ifcase#1\or 03.08.2023% - \or 07.08.2023% - \or 09.08.2023% + дата~обращения: \ifcase#1% + \or 03.08.2023% 1 (default) + \or 07.08.2023% 2 + \or 09.08.2023% 3 + \or 19.08.2023% 4 + \or 17.10.2023% 5 \else\@ctrerr\fi - } +} \renewcommand\bibname{Библиографический список} \begin{thebibliography}{00} \addcontentsline{toc}{chapter}{Библиографический список} - \bibitem{book:nm-examples} Ахмадиев Ф.Г., Габбасов Ф.Г., Ермолаева Л.Б., Маланичев И.В. Численные методы. Примеры и задачи. Учебно-методическое пособие по курсам «Информатика» и «Вычислительная математика». -- Казань: - КГАСУ, 2017. -- 107 с. - \bibitem{book:bahvalov} Бахвалов~Н.~С., Жидков~Н.~П., - Кобельков~Г.~М. Численные методы. -- 7-е изд. -- М.: БИНОМ. Лаборатория знаний, - 2011. -- 636 с., c илл. -- (Классический университетский учебник). + \bibitem{book:nm-examples} Ахмадиев Ф.Г., Габбасов Ф.Г., Ермолаева Л.Б., Маланичев И.В. Численные методы. Примеры и задачи. Учебно-методическое пособие по курсам «Информатика» и «Вычислительная математика». -- Казань: КГАСУ, 2017. -- 107 с. + \bibitem{book:bahvalov} Бахвалов~Н.~С., Жидков~Н.~П., Кобельков~Г.~М. Численные методы. -- 7-е изд. -- М.: БИНОМ. Лаборатория знаний, 2011. -- 636 с., c илл. -- (Классический университетский учебник). + \bibitem{links:matplotlib} Документация модуля pyplot библиотеки matplotlib [Электронный ресурс] -- URL:~\url{https://matplotlib.org/stable/api/pyplot_summary.html} (\LiteratureAccessDate[5]). + \bibitem{book:levitin} Левитин~A.~В. Алгоритмы: введение в разработку и анализ. -- Пер.~с~англ. -- М.:Издательский~дом~"Вильяме", 2006. -- 576 с., с ил. \bibitem{book:lectures} Письменный~Д.~Т. Конспект лекций по высшей математике. 2 часть. -- М.: Рольф, 2000. -- 256 с., с илл. + \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 \bibitem{links:numpy} Numpy. Официальный сайт проекта [Электронный ресурс] -- URL:~\url{https://numpy.org/} (\LiteratureAccessDate[2]). \bibitem{links:numpy_doc} Numpy API Reference [Электронный ресурс] -- URL:~\url{https://numpy.org/doc/stable/reference/index.html} (\LiteratureAccessDate[3]). + \bibitem{links:PEP465} PEP 465 -- A dedicated infix operator for matrix multiplication [Электронный ресурс] -- URL:~\url{A dedicated infix operator for matrix multiplication} (\LiteratureAccessDate[4]). + \bibitem{links:bhatia} Positive definite matrices / R. Bhatia // Princeton Series in Applied Mathematics -- 2007. \bibitem{links:python} Python. Официальный сайт проекта [Электронный ресурс] -- URL:~\url{https://www.python.org/} (\LiteratureAccessDate). \bibitem{links:scipy} Scipy. Официальный сайт проекта [Электронный ресурс] -- URL:~\url{https://scipy.org/} (\LiteratureAccessDate[2]). \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{links:tiobe_index} TIOBE. Официальный сайт проекта [Электронный ресурс] -- URL:~\url{https://www.tiobe.com/tiobe-index/} (\LiteratureAccessDate). \end{thebibliography}