Merge pull request 'Добавление информации о методах (№1)' (#2) from dev into master

Reviewed-on: #2

6SST_NM.code-workspace

@@ -13,6 +13,12 @@
                "pdflatex"
             ]
          },
+         {
+            "name": "latexmk (xelatex)",
+            "tools": [
+               "xelatexmk"
+            ]
+         },
          {
             "name": "latexmk",
             "tools": [
@@ -34,12 +40,6 @@
                "pdflatex"
             ]
          },
-         {
-            "name": "tectonic",
-            "tools": [
-               "tectonic"
-            ]
-         }
       ],
       "latex-workshop.latex.tools": [
          {
@@ -50,6 +50,7 @@
                "-interaction=nonstopmode",
                "-file-line-error",
                "-pdf",
+               "-shell-escape",
                "-outdir=%OUTDIR%",
                "%DOC%"
             ],
@@ -81,6 +82,7 @@
             "command": "pdflatex",
             "args": [
                "-synctex=1",
+               // "-recorder",
                "-shell-escape",
                "-interaction=nonstopmode",
                "-file-line-error",
@@ -109,5 +111,6 @@
          }
       ],
       "latex-workshop.latex.rootFile.indicator": "\\begin{document}",
+      "latex-workshop.latex.recipe.default": "lastUsed",
    }
 }

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}

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()

conclusion.tex

@@ -0,0 +1,15 @@
+\chapter*{Заключение}
+
+В ходе выполнения работы были выполнены поставленные задачи, цель
+достигнута. В ходе исследования возможностей библиотек выяснилось,
+что они не содержат реализаций численных методов для
+решения систем нелинейных уравнений.
+
+В целом, в ходе работы с библиотеками стоит отметить положительные
+стороны: гибкость интерфейсов, высокое качество решений, большое
+количество доступной информации о результатах работы алгоритма.
+
+Из недостатков интерфейса библиотек можно отметить его неоднородность ---
+некоторые численные методы реализованы с помощью классов, другие
+с помощью функций, что может незначительно увеличить время на
+освоение интерфейса библиотек и ее возможностей пользователями.

config.tex

@@ -5,6 +5,7 @@
 \usepackage{threeparttable}
 \usepackage[labelsep=endash,tableposition=top,labelfont=md,textfont=md]{caption}
 \usepackage{indentfirst}
+\usepackage{float}
 \usepackage{rotating}
 \usepackage{setspace}
 \usepackage{array}
@@ -111,6 +112,73 @@
    tabsize=3,
    escapechar={|},
    emptylines=0,
+   extendedchars=true,
+   literate={а}{{\cyra}}1
+   {б}{{\cyrb}}1
+   {в}{{\cyrv}}1
+   {г}{{\cyrg}}1
+   {д}{{\cyrd}}1
+   {е}{{\cyre}}1
+   {ё}{\"{\cyre}}1
+   {ж}{{\cyrzh}}1
+   {з}{{\cyrz}}1
+   {и}{{\cyri}}1
+   {й}{{\cyrishrt}}1
+   {к}{{\cyrk}}1
+   {л}{{\cyrl}}1
+   {м}{{\cyrm}}1
+   {н}{{\cyrn}}1
+   {о}{{\cyro}}1
+   {п}{{\cyrp}}1
+   {р}{{\cyrr}}1
+   {с}{{\cyrs}}1
+   {т}{{\cyrt}}1
+   {у}{{\cyru}}1
+   {ф}{{\cyrf}}1
+   {х}{{\cyrh}}1
+   {ц}{{\cyrc}}1
+   {ч}{{\cyrch}}1
+   {ш}{{\cyrsh}}1
+   {щ}{{\cyrshch}}1
+   {ъ}{{\cyrhrdsn}}1
+   {ы}{{\cyrery}}1
+   {ь}{{\cyrsftsn}}1
+   {э}{{\cyrerev}}1
+   {ю}{{\cyryu}}1
+   {я}{{\cyrya}}1
+   {А}{{\CYRA}}1
+   {Б}{{\CYRB}}1
+   {В}{{\CYRV}}1
+   {Г}{{\CYRG}}1
+   {Д}{{\CYR96}}1
+   {Е}{{\CYRE}}1
+   {Ё}{{\"{\CYRE}}}1
+   {Ж}{{\CYRZH}}1
+   {З}{{\CYRZ}}1
+   {И}{{\CYRI}}1
+   {Й}{{\CYRISHRT}}1
+   {К}{{\CYRK}}1
+   {Л}{{\CYRL}}1
+   {М}{{\CYRM}}1
+   {Н}{{\CYRN}}1
+   {О}{{\CYRO}}1
+   {П}{{\CYRP}}1
+   {Р}{{\CYRR}}1
+   {С}{{\CYRS}}1
+   {Т}{{\CYRT}}1
+   {У}{{\CYRU}}1
+   {Ф}{{\CYRF}}1
+   {Х}{{\CYRH}}1
+   {Ц}{{\CYRC}}1
+   {Ч}{{\CYRCH}}1
+   {Ш}{{\CYRSH}}1
+   {Щ}{{\CYRSHCH}}1
+   {Ъ}{{\CYRHRDSN}}1
+   {Ы}{{\CYRERY}}1
+   {Ь}{{\CYRSFTSN}}1
+   {Э}{{\CYREREV}}1
+   {Ю}{{\CYRYU}}1
+   {Я}{{\CYRYA}}1
 }
 \newcommand\inputcode[3][]{
    {\bfseries #1\filename{#2}}:

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}
+
+Данные методы самостоятельно вызываются при запуске программы
+пользователем.
+
+Вывод терминала а также графики для остальных методов приведены 
+в приложении А.
+
+Код программы приведен в приложении Б.

main.tex

@@ -1,3 +1,4 @@
+% TODO: Fix terminology usage of "якобиан"
 \input{vars}
 \input{config}
 \sloppy
@@ -114,10 +115,10 @@
          противоположные знаки.
    \item \(a\) --- scalar
 
-         Первый конец интервала \([a,b]\).
+         Первый конец интервала \([a;b]\).
    \item \(b\) --- scalar
 
-         Второй конец интервала \([a,b]\).
+         Второй конец интервала \([a;b]\).
    \item \(xtol\) --- number, необязательный
 
          Вычисленный корень \(x0\) будет удовлетворять
@@ -485,7 +486,8 @@ LU-разложение \cite[с. 259]{book:levitin}. Для получения
 В scipy для решения СЛУ вида (\ref{formula:eqn_banded_system})
 существует две функции в модуле \textbf{scipy.linalg}:
 \textbf{solve\_banded} \cite{links:scipy_doc} и
-\textbf{solveh\_banded} \cite{links:scipy_doc}.
+\textbf{solveh\_banded} \cite{links:scipy_doc}, обе из которых могут 
+решать задачи для n-диагональных матриц.
 
 Различие между ними заключается в том, что \textbf{solve\_banded}
 не использует метод прогонки, из-за низкой устойчивости метода в общем
@@ -554,41 +556,57 @@ LU-разложение \cite[с. 259]{book:levitin}. Для получения
 Обе функции принимают матрицу \(ab\) либо в верхней (\textbf{solve\_banded}), либо в нижней форме (\textbf{solveh\_banded} при
 включенной опции \(lower\)). Например, для матрицы
 
-\begin{tabular}[htpb]{ccccc}
-   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{tabular}\\
+\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{tabular}[htpb]{ccccc}
-   0 & 0 & -1 & -1 & -1 \\
-   0 & 2 & 2  & 2  & 2  \\
-   5 & 4 & 3  & 2  & 1  \\
-   1 & 1 & 1  & 1  & 0  \\
-\end{tabular}
+\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{tabular}[htpb]{cccccc}
-   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{tabular}\\
+\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{tabular}[htpb]{cccccc}
-   4  & 5  & 6  & 7  & 8 & 9 \\
-   2  & 2  & 2  & 2  & 2 & 0 \\
-   -1 & -1 & -1 & -1 & 0 & 0 \\
-\end{tabular}
+\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{Описание метода}
@@ -797,7 +815,7 @@ LU-разложение \cite[с. 259]{book:levitin}. Для получения
 последний --- с помощью сглаживания.
 \subsection{Интерполяционный полином Лагранжа}
 \subsubsection{Описание метода}
-Для известных пар значений \((x_i,y_i), x_i \in [a,b]; i = 1,2,3,\dots,n\)
+Для известных пар значений \((x_i,y_i), x_i \in [a;b]; i = 1,2,3,\dots,n\)
 строится полином \(P(x)\), удовлетворяющий условию \(P(x_i) = y_i\).
 
 Полином \(P(x)\) определяется следующей формулой:
@@ -1164,6 +1182,11 @@ LU-разложение \cite[с. 259]{book:levitin}. Для получения
          минимизации выражения \(||ax-b||^2\), где \(a,b\) --- матрицы
          параметров; функция \textbf{nnls} отличается от нее
          ограничением на положительность искомых значений.
+
+         Подходят только для нахождения параметров многочленов,
+         не буду рассматриваться в данной работе по причине
+         сложной подготовки входных данных для решения задачи
+         метода, а также их узкой специализации.
 \end{enumerate}
 
 Функция \textbf{curve\_fit} имеет следующий набор параметров:
@@ -1268,7 +1291,7 @@ LU-разложение \cite[с. 259]{book:levitin}. Для получения
          \verb|False|.
    \item \(nan\_policy\) --- {\verb|"raise"|, \verb|"omit"|, \verb|None|}, необязательный
 
-         Определяет, как действовать, если входные данные содержат nan. Доступны следующие параметры (по умолчанию — \verb|None|):
+         Определяет, как действовать, если входные данные содержат \verb|NaN|. Доступны следующие параметры (по умолчанию --- \verb|None|):
          \begin{itemize}
             \item \verb|"raise"|: выдает ошибку
 
@@ -1308,7 +1331,7 @@ LU-разложение \cite[с. 259]{book:levitin}. Для получения
          дополнительную информацию о работе алгоритма.
    \item \(mesg\) --- str
 
-         Cтроковое сообщение с информацией о решении.
+         Строковое сообщение с информацией о решении.
    \item \(ier\) --- int
 
          Целочисленный флаг. Если он равен 1, 2, 3 или 4, решение
@@ -1338,7 +1361,7 @@ LU-разложение \cite[с. 259]{book:levitin}. Для получения
          будет рассматриваться как одномерный массив с одним элементом.
    \item \(jac\) --- \{\verb|"2-point"|, \verb|"3-point"|, \verb|"cs"|, callable\}, необязательный
 
-         Метод вычисления матрицы Якоби (матрица размером \(m \times  n\), где элемент \((i, j)\) является частной производной f[i] по x[j]). Ключевые слова выбирают схему конечных разностей для числовой оценки. Схема \verb|"3-point"| более точная, но требует в два раза больше операций, чем \verb|"2-point"| (по умолчанию). Схема \verb|"cs"| использует комплексные шаги и, хотя потенциально является наиболее точной, она применима только тогда, когда \(fun\) правильно обрабатывает комплекснозначные входные данные и может быть аналитически продолжена на комплексную плоскость. \(method\) \verb|"lm"| всегда использует схему \verb|"2-point"|. Если параметр имеет тип \verb|callable|, он используется как \verb|jac(x, *args, **kwargs)| и должен возвращать хорошее приближение (или точное значение) для якобиана в виде array\_like (к выводу \(jac\) будет применен \verb|np.atleast_2d|), разреженной матрицы (предпочтительна \verb|csr_matrix| для производительности) или \verb|scipy.sparse.linalg.LinearOperator|.
+         Метод вычисления матрицы Якоби (матрица размером \(m \times  n\), где элемент \((i, j)\) является частной производной \(f[i]\) по \(x[j]\)). Ключевые слова выбирают схему конечных разностей для числовой оценки. Схема \verb|"3-point"| более точная, но требует в два раза больше операций, чем \verb|"2-point"| (по умолчанию). Схема \verb|"cs"| использует комплексные шаги и, хотя потенциально является наиболее точной, она применима только тогда, когда \(fun\) правильно обрабатывает комплекснозначные входные данные и может быть аналитически продолжена на комплексную плоскость. \(method\) \verb|"lm"| всегда использует схему \verb|"2-point"|. Если параметр имеет тип \verb|callable|, он используется как \verb|jac(x, *args, **kwargs)| и должен возвращать хорошее приближение (или точное значение) для якобиана в виде array\_like (к выводу \(jac\) будет применен \verb|np.atleast_2d|), разреженной матрицы (предпочтительна \verb|csr_matrix| для производительности) или \verb|scipy.sparse.linalg.LinearOperator|.
    \item \(bounds\) --- tuple из 2 array\_like / \verb|Bounds|, необязательный
 
          Есть два способа указать границы:
@@ -1408,44 +1431,497 @@ LU-разложение \cite[с. 259]{book:levitin}. Для получения
          \begin{itemize}
             \item Для  \verb|"trf"| и  \verb|"dogbox"|: \(100 \cdot n\).
 
-            \item Для \verb|"lm"|: \(100 \cdot n\), если \(jac\) является вызываемым, и \(100 \cdot n \cdot (n + 1)\) в противном случае (поскольку \verb|"lm"| учитывает вызовы функций в оценке Якобиана).
+            \item Для \verb|"lm"|: \(100 \cdot n\), если \(jac\) является вызываемым, и \(100 \cdot n \cdot (n + 1)\) в противном случае (поскольку \verb|"lm"| учитывает вызовы функций в оценке якобиана).
          \end{itemize}
    \item \(diff\_step\) --- \verb|None| / array\_like, необязательный
 
+         Определяет относительный размер шага для аппроксимации якобиана конечной разностью. Фактический шаг вычисляется как \verb|x * diff_step|. Если \verb|None| (по умолчанию), то \(diff_step\) принимается за обычную "оптимальную" степень машинного эпсилона для используемой схемы конечных разностей.
+   \item \(tr\_solver\) --- {\verb|None|, \verb|"exact"|, \verb|"lsmr"|}, необязательный
 
-   \item 
+         Метод решения подзадач доверительной области, применим только для методов \verb|"trf"| и \verb|"dogbox"|.
+         \begin{itemize}
+            \item \verb|"exact"| подходит для не очень больших задач с плотными матрицами Якоби. Вычислительная сложность на итерацию сравнима с сингулярным разложением матрицы якобиана.
+         \end{itemize}
+   \item \(tr\_options\) --- dict, необязательный
+
+         Параметры ключевых слов передаются в решатель доверительной
+         области.
+         \begin{itemize}
+            \item \verb|tr_solver="exact"|: \(tr\_options\) игнорируются.
+            \item \verb|tr_solver="lsmr"|: опции для \textbf{scipy.sparse.linalg.lsmr}.
+                  Кроме того, \(method\)=\verb|"trf"| поддерживает опцию \verb|'regularize'| (bool, по умолчанию --- \verb|True|), которая добавляет член регуляризации к нормальному уравнению, что улучшает сходимость, если якобиан имеет недостаточный ранг.
+         \end{itemize}
+   \item \(jac\_sparsity\) --- {\verb|None|, array\_like, разреженная матрица}, необязательный
+
+         Определяет структуру разреженности матрицы якобиана для
+         конечно-разностной оценки, ее форма должна быть (m, n).
+         Если якобиан имеет лишь несколько ненулевых элементов в
+         каждой строке, обеспечение разреженной структуры значительно
+         ускорит вычисления. Нулевая запись означает, что
+         соответствующий элемент якобиана тождественно равен нулю.
+         Если предусмотрено, принудительно используется решатель
+         доверительной области \verb|"lsmr"|. Если \verb|None|
+         (по умолчанию), будет использоваться плотная разность.
+         Не имеет эффекта при \(method\)=\verb|"lm"|.
+   \item \(verbose\) --- {0, 1, 2}, необязательный
+
+         Уровень детализации алгоритма:
+         \begin{itemize}
+            \item 0 (по умолчанию): бесшумная работа функции.
+            \item 1: отображение отчета о завершении.
+            \item 2: отображение хода выполнения во время итераций
+                  (не поддерживается методом \verb|"lm"|).
+         \end{itemize}
+   \item \(args, kwargs\) --- tuple / dict, необязательный
+
+         Дополнительные аргументы передаваемые в \(fun\) и \(jac\).
+         Оба пусты по умолчанию. Сигнатура вызова для \(fun\) ---
+         \(fun(x, *args, **kwargs)\), та же самая для \(jac\).
 \end{enumerate}
-\section{Численное интегрирование}
-\subsection{Метод прямоугольников}
-\subsubsection{Описание метода}
-\subsubsection{Реализации метода в библиотеках numpy, scipy}
+Функция возвращает результат в виде экземпляра класса \verb|OptimizeResult|, у которого определены поля \(x, cost, fun, jac, grad, optimality, active\_mask, nfev, njev,status, \) \(message, success\).
 
+
+\section{Численное интегрирование}
+Задача численного интегрирования состоит в вычислении определенного
+интеграла:
+\begin{equation}
+   I = \int_{a}^{b} f(x) \, dx
+   \label{formula:nm-integral1}
+\end{equation}
+
+Для решения данной задачи выберем на отрезке \([a;b]\) \(n\)
+различных узлов:
+\(a = x_0 < x_1 < x_2 < ... < x_{n-1} < x_n = b\), затем по выбранным
+узлам интерполируем \(f(x)\). Если в роли интерполяционной функции
+выбрать полином \(P_m(x)\), то интеграл (\ref{formula:nm-integral1})
+можно вычислить по формуле
+\begin{equation*}
+   I \approx \int_{a}^{b} P_m(x) \, dx \quad.
+\end{equation*}
+
+Данную формулу также называют квадратурной формулой интерполяционного
+типа.
+
+К наиболее распространенным методам относят метод прямоугольников,
+метод трапеций и метод парабол (Симпсона). Первый метод использует
+использует полиномы 0 степени, второй --- 1 степени, третий ---
+2 степени.
+
+В данной работе будет рассмотрены методы парабол и трапеций.
+
+Метод парабол будет точнее других в большинстве случаев
+(он позволяет находить точное решение для любых
+f\(x\), если \(f^{(4)}(x) = 0, x \in [a;b] \), в соответствии с
+формулой Джузеппе Пеано). Из его недостатков можно отметить низкую
+точность на пилообразных функциях (т.е. значение которой резко
+возрастают на отрезках малой длины).
+
+При этом, если количество точек, по которым строится \(P_m(x)\), четно, то метод трапеций может оказаться удобнее, тем самым прекрасно дополняя
+метод парабол.
 \subsection{Метод трапеций}
 \subsubsection{Описание метода}
+На отрезке \([a;b]\) выбирается \(n\) узлов, на отрезках \([x_i;x_{i+1}]\) строятся интерполяционные многочлены 1 степени \(P_1(x)\).
+Если используется метод Лагранжа, то \(P_1(x)\) определяется так:
+\begin{equation*}
+   P_1(x) = f(x_i) \frac{x-x_{i+1}}{x_i-x_{i+1}} + f(x_{i+1}) \frac{x-x_i}{x_{i+1}-x_i}
+\end{equation*}
+При этом
+\begin{equation*}
+   \int_{x_i}^{x_{i+1}} P_1(x) \, dx = \frac{1}{2} \left(f(x_i)+f(x_{i+1})\right)(x_{i+1} - x_i)
+\end{equation*}
+В случае равноотстоящих узлов
+\(x_0, x_1 = x_0 + h, ..., x_n = x_0 + nh\) значение интеграла будет таким:
+\begin{equation*}
+   I \approx \frac{h}{2} \sum_{i=0}^{n-1} (f(x_i) + f(x_{i+1}))
+\end{equation*}
 \subsubsection{Реализации метода в библиотеках numpy, scipy}
+В библиотеке scipy найдена функция \textbf{trapezoid} в модуле
+\textbf{scipy.integrate}. Она имеет следующие параметры:
+\begin{enumerate}
+   \item \(y\) --- array\_like
 
+         Массив входных данных по которым будет вычисляться интеграл
+   \item \(x\) --- array\_like, необязательный
+
+         Массив значений \(x\), соответствующих \(y\). Если \verb|None|
+         (по умолчанию), то в роли \(x\) будет создан массив равноотстоящих
+         значений (\(h = dx\))
+   \item \(dx\) --- scalar, необязательный
+
+         Расстояние между соседними значениями \(x\). Значение используется
+         когда \(x\)=\verb|None|. По умолчанию --- 1.
+   \item \(axis\) --- int, необязательный
+
+         Ось, относительно которой производить вычисление интеграла. По
+         умолчанию --- \(-1\).
+\end{enumerate}
+
+Данная функция возвращает число, если \(y\) одномерный и массив
+размерности \((n-1)\) для \(n\)-мерного \(y\).
 \subsection{Метод парабол (Симпсона)}
 \subsubsection{Описание метода}
+На отрезке \([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
 
-\section{Численное решение обыкновенных дифференциальных уравнений}
+         Массив, который необходимо интегрировать.
+   \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}
 
-\chapter{Экспериментальное исследование возможностей библиотек}
+Функция \textbf{solve\_ivp} может найти решение задачи Коши для системы
+ОДУ, и включает в себя реализации нескольких методов. Далее будут
+описаны только те значения параметров, которые необходимы для решения
+задачи исследуемым методом (Рунге-Кутта). Также стоит учесть, что при
+описании параметров, вместо \(x\) будет использоваться \(t\), как и
+принято в зарубежных источниках.
+Данная функция имеет следующие параметры:
+\begin{enumerate}
+   \item \(fun\) --- callable
 
-\chapter*{Заключение}
+         Правая часть системы: производная \(\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}

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}

sources.tex

@@ -4,25 +4,25 @@
    \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:bhatia} R. Bhatia Positive definite matrices -- Princeton Series in Applied Mathematics -- 2007.
    \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}