mathmod - библиотека программ для численного моделирования

Содержание

• mathmod - библиотека программ для численного моделирования
  • Содержание
    • Теория погрешностей
  • Решение нелинейных уравнений
    • Постановка задачи
    • Основные вопросы
    • 1. Метод Ньютона
    • 2. Упрощённый метод Ньютона
    • 3. Метод секущих
    • 4. Метод ложного положения
    • 5. Метод бисекций
    • 6. Метод простой итерации
  • Решение систем линейных алгебраических уравнений (СЛАУ)
    • Постановка задачи
    • Основные вопросы
    • Прямые методы
      • 1. Метод Гаусса (схема единственного деления)
      • 2. Метод Гаусса (схема частичного выбора)
      • 3. Метод Холецкого (LLT-разложение)
      • 4. Метод LU-разложения
      • 5. Метод прогонки
    • Итерационные методы
      • 1. Метод Якоби
      • 2. Метод Гаусса-Зейделя
      • 3. Метод релаксации
    • Приближение функций
      • Постановка задачи
      • Постановка задачи интерполяции
      • Полиномиальная интерполяция. Многочлен Лагранжа
      • Интерполяционный многочлен Ньютона с конечными разностями
      • Метод наименьших квадратов
  • Установка

---

Теория погрешностей

Пусть $a$ - точное решение. $a^{\ast }$ - приближенное значение. Тогда абсолютная погрешность:

Абсолютная погрешность:

$$\mathrm{\Delta }(a^{\ast })=|a-a^{\ast }|$$

Относительная погрешность:

$$\delta (a^{\ast })=\frac{|a-a^{\ast }|}{|a|}=\frac{\mathrm{\Delta }(a^{\ast })}{|a|}$$

Верная цифра - значащая цифра числа $a^{\ast }$, абсолютная погрешность которой не превосходит единицы разряда, соответствующей этой цифре.

Округление - замена числа $a$ его другим числом $a^{\ast }$ с меньшим количеством значащих цифр.

При округлении возникает погрешность округления.

Решение нелинейных уравнений

Постановка задачи

Задача отыскания корней нелинейного уравнения с одним неизвестным вида:

$$f(x)=0$$

Корень уравнения - значение $\bar{x}$, при котором $f(\bar{x})=0$. Геометрически, $\bar{x}$ является абсциссой точки пересечения графика $f(x)$ с осью $Ox$.

Для заданной точности $ϵ$ требуется найти приближенное значение корня ${\bar{x}}^{\ast }$ такое, что:

$$|\bar{x}-{\bar{x}}^{\ast }|\le ϵ$$

Пусть функция $f(x)$ дифференцируема $m>1$ раз в точке $\bar{x}$, тогда:

Простой корень - корень уравнения $\bar{x}$ называется простым, если $f^{\prime }(\bar{x})\ne 0$.

Кратный корень - корень уравнения $\bar{x}$ называется простым, если $f^{\prime }(\bar{x})=0$. Целое число $m$ назовем кратностью корня $\bar{x}$, если $f^{(k)}(\bar{x})=0$, для $k=1,2,...,m-1$ и $f^{(m)}(\bar{x})\ne 0$.

Этапы решения уравнения:

1. Локализация корня;
2. Вычисление корня с заданной точностью.

Отрезок локализации - отрезок $[a,b]$, который содержит 1 корень уравнения.

Основные вопросы

Сходящийся итерационный процесс - метод отыскания корня $\bar{x}$, заключающийся в построении последовательности приближений к нему $x^{(0)},x^{(1)},...,,x^{(k)},...$ такой, что:

$$\underset{k\to \mathrm{\infty }}{lim}{x}^{(k)}=\bar{x}$$

Порядок сходимости:

Пусть в некоторой малой окрестности корня $\bar{x}$ уравнения $f(x)=0$ итерационная поледовательность удовлетворяет неравенству:

$$|\bar{x}-x^{(k)}|\le C|\bar{x}-x^{(k-1)}{|}^p$$

, где $С>0$ и $p\ge 1$ - постоянные. Тогда $p$ называется порядком сходимости.

Порядок сходимости - скорость уменьшения погрешности между последовательными приближениями решения.

Одношаговый итерационный метод - метод, у которого очередное приближение $x^{(k+1)}$ находится только через одно предыдущее $x^{(k)}$. Для его работы нужно знать только одно начальное приближение $x^{(0)}$. (Пример - метод Ньютона)

Многошаговый итерационный метод ($l$ - шаговый) - метод, у которого очередное приближение $x^{(k+1)}$ находится $l$ предыдущих $x^{(k)},x^{(k-1)},...,x^{(k-l+1)}$. Для него следует задать $l$ начальных приближений $x^{(0)},x^{(1)},...,x^{(l-1)}$. (Пример - метод секщих)

Интервал неопределенности - окрестность $(\bar{x}-ϵ,\bar{x}+ϵ)$ внутри которой любую точку можно принять за приближение к корню.

1. Метод Ньютона

• Быстрый итерационный метод для нахождения корня уравнения $f(x)=0$.
• Требует предоставления функции $f(x)$ и её производной $f^{\prime }(x)$.
• Функция: newton(f, df, x, epsilon=1e-6)
• Описание параметров:
  • f - Функция, корень которой нужно найти;
  • df - Производная функции;
  • x - Начальное приближение корня;
  • epsilon - Заданная точность (по умолчанию ${10}^{-6}$).

from mathmod.nonlinear_equations import newton

x, iteration = newton(f, df, x, epsilon=1e-6)

Расчетная формула:

$$x^{(n+1)}=x^{(n)}-\frac{f(x^{(n)})}{f^{\prime }(x^{(n)})}$$

Сходимость метода: квадратичная (при выборе начального приближения из достаточно малой окрестности).

Критерий окончания:

$$|x^{(n+1)}-x^{(n)}|<\epsilon$$

2. Упрощённый метод Ньютона

• Упрощённая версия метода Ньютона, где производная функции вычисляется только один раз.
• Функция: simplified_newton(f, df, x0, epsilon=1e-6)
  
• Описание параметров:
  • f - Функция, корень которой нужно найти;
  • df - Производная функции;
  • x - Начальное приближение корня;
  • epsilon - Заданная точность (по умолчанию ${10}^{-6}$).

from mathmod.nonlinear_equations import simplified_newton

x, iteration = simplified_newton(f, df, x0, epsilon=1e-6)

Расчетная формула:

$$x^{(n+1)}=x^{(n)}-\frac{f(x^{(n)})}{f^{\prime }(x^{(0)})}$$

Сходимость метода: скорость сходимости тем выше, чем ближе начальное приближение $x^{(0)}$ к решению $\bar{x}$.

Критерий окончания:

$$|x^{(n+1)}-x^{(n)}|<\epsilon$$

3. Метод секущих

• Не требует аналитической производной функции.
• Использует приближённую производную.
• Функция: secant(f, x_minus_1, x_n, epsilon=1e-6)
• Описание параметров:
  • f - Функция, корень которой нужно найти;
  • x_minus_1 - Первой начальное приближение;
  • x - Второе начальное приближение;
  • epsilon - Заданная точность (по умолчанию ${10}^{-6}$).

from mathmod.nonlinear_equations import secant

x, iteration = secant(f, x_minus_1, x_n, epsilon=1e-6)

Расчетная формула:

$$x^{(n+1)}=x^{(n)}-\frac{x^{(n-1)}-x^{(n)}}{f(x^{(n-1)})-f(x^{(n)})}f(x^{(n)})$$

Сходимость метода: $p=\frac{\sqrt{5}+1}{2}\approx 1.618$ - сверхлинейная, если вычисляется простой корень. При неудачном выборе приближения, метод расходится. Поэтому требуется выбор двух близких к $\bar{x}$ начальных приближений $x^{(0)}$ и $x^{(1)}$.

Критерий окончания:

$$|x^{(n+1)}-x^{(n)}|<\epsilon$$

4. Метод ложного положения

• Функция: false_position(f, a, b, epsilon=1e-6)
• Описание параметров:
• f: Функция, корень которой нужно найти;
• a: Левый конец начального отрезка локализации корня;
• b: Правый конец начального отрезка локализации корня;
• epsilon: Заданная точность (по умолчанию ${10}^{-6}$).

from mathmod.nonlinear_equations import false_position

x, iteration = false_position(f, a, b, epsilon=1e-6)

Расчетная формула:

$$x^{(n+1)}=x^{(n)}-\frac{c-x^{(n)}}{f(c)-f(x^{(n)})}f(x^{(n)})$$

, где $c$ - некторая точнка из окрестности корня.

Сходимость метода: линейная. Для достижения заданной точности требуется тем меньше итераций, чем ближе к корню лежит точка $c$.

Критерий окончания:

$$|x^{(n+1)}-x^{(n)}|<\epsilon$$

5. Метод бисекций

• Работает на основе деления отрезка, учитывая значения функции.
• Требует, чтобы начальный отрезок $[a,b]$ удовлетворял условию $f(a)\ast f(b)<0$. Функция: bisection(f, a, b, epsilon)

Описание параметров:

• f - Функция, корень которой нужно найти;
• a - Левый конец начального отрезка локализации корня;
• b - Правый конец начального отрезка локализации корня;
• epsilon - Заданная точность (по умолчанию ${10}^{-6}$).

from mathmod.nonlinear_equations import bisection

x, iteration = bisection(f, a, b, epsilon=1e-6)

Расчетная формула:

$$c=\frac{a+b}{2}$$

Сходимость метода: со скоростью геометрической прогрессии, знаменатель которой $q=\frac{1}{2}$.

Критерий окончания:

$$b^{(n)}-a^{(n)}<2\epsilon$$

Тогда $x^{(n)}=\frac{a^{(n)}+b^{(n)}}{2}$ является искомым приближением к корню с точностью $ϵ$.

6. Метод простой итерации

Функция: simple_iteration_method(phi, x0, epsilon=1e-6)

Описание параметров:

• phi - Итерационная функция $\varphi (x)$, преобразующая исходное уравнение $f(x)=0$ к виду $x=\varphi (x)$;
• x0 - Начальное приближение корня;
• epsilon - Заданная точность (по умолчанию ${10}^{-6}$).

from mathmod.nonlinear_equations import simple_iteration_method

x, iteration = simple_iteration_method(phi, x0, epsilon=1e-6)

Расчетная формула:

$$x^{(n+1)}=\varphi (x^{(n)})$$

Теорема о сходимости:

Если в окрестности корня функция $\varphi (x)$ непрерывно дифференцируема и удовлетворяет условию:

$$\underset{x\in [a,b]}{max}|\varphi (x)|\le q$$

где $0\le q<1$ - постоянная

Тогда независимо от выбора начального приближения $x^{(0)}$ из указанной окрестности корня итерационная последовательность не выходит из этой окрестности, метод сходится со скоростью геометрической прогресии и справедлива следующая оценка погрешности (априорная оценка):

Априорная оценка - показывет, что итерационный метод сходится

$$|x^{(n)}-\bar{x}|\le q^n|x^{(0)}-\bar{x}|$$

Чем меньше $q$, тем выше скорость сходимости.

Апостериорная оценка - критерий окончания итерационного процесса

$$|x^{(n)}-x^{(n-1)}|\le \frac{1-q}{q}ϵ$$

Если это условие выполнено, то можно считать, что $x^{(n)}$ является приближением к $\bar{x}$ с точностью $ϵ$.

Решение систем линейных алгебраических уравнений (СЛАУ)

Постановка задачи

Система уравнений в общем виде:

$$\begin{array}{rl}{a}_{11}{x}_1+a_{12}{x}_2+a_{13}{x}_3+\cdots +a_{1m}{x}_m& =b_1,\\ a_{21}{x}_1+a_{22}{x}_2+a_{23}{x}_3+\cdots +a_{2m}{x}_m& =b_2,\\ a_{31}{x}_1+a_{32}{x}_2+a_{33}{x}_3+\cdots +a_{3m}{x}_m& =b_3,\\ & \dots \\ a_{m1}{x}_1+a_{m2}{x}_2+a_{m3}{x}_3+\cdots +a_{mm}{x}_m& =b_m.\end{array}$$

В матричной форме эта система принимает вид:

$$A\mathbf{x}=\mathbf{b}$$

, где

$$A=\left[\begin{array}{ccccc}{a}_{11}& a_{12}& a_{13}& \dots & a_{1m}\\ a_{21}& a_{22}& a_{23}& \dots & a_{2m}\\ a_{31}& a_{32}& a_{33}& \dots & a_{3m}\\ \dots & \dots & \dots & \dots & \dots \\ a_{m1}& a_{m2}& a_{m3}& \dots & a_{mm}\end{array}\right],\mathbf{x}=\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ ⋮\\ x_m\end{array}\right],\mathbf{b}=\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\\ ⋮\\ b_m\end{array}\right]$$

$m$ - порядок матрицы;

$A$ - невырожденная матрица ($\mathrm{\Delta }A\ne 0$, т.е. $\mathrm{\exists }$ единственное решение);

$x$ - вектор неизвестных;

$b$ - вектор свободных членов;

Пусть $\bar{x}$ - точное решение, $x^{\ast }$ - приближенное решение.

$$\bar{\mathbf{x}}=\left[\begin{array}{c}{\bar{x}}_1\\ {\bar{x}}_2\\ {\bar{x}}_3\\ ⋮\\ {\bar{x}}_m\end{array}\right]{\mathbf{x}}^{\ast }=\left[\begin{array}{c}{x}_1^{\ast }\\ x_2^{\ast }\\ x_3^{\ast }\\ ⋮\\ x_m^{\ast }\end{array}\right]$$

Тогда вектор $ϵ=\bar{x}-x^{\ast }$ называется вектором погрешности.

Задача:

Найти решение системы $Ax=b$ с точностью $ϵ$. Это означает, что нужно найти вектор $x^{\ast }$ такой, что $||\bar{x}-x^{\ast }||\le ϵ$, где $||\cdot ||$ - одна из норм (единичная, евклидова, бесконечности).

Основные вопросы

Прямой метод - метод, который позволяет получить решение после выполнения конечного числа элементарных операций;

Итерационный метод - метод, который строит последовательность приближений к решению;

Норма - будем говорить, что в $R^m$ введена норма, если каждому вектору $x\in R^m$ сопоставлено вещественное число, обозначаемое $||x||$.

1. Нормы векторов:

   

   Свойства норм векторов:

   

2. Нормы матриц:

   

   Свойства норм матриц:

   

3. Абсолютна и относительная погрешность векторов и матриц:

• Погрешность векторов:

  $\mathrm{\Delta }{x}^{\ast }=||\bar{x}-x^{\ast }||$ - абсолютная погрешность;

  $\delta x^{\ast }=\frac{\mathrm{\Delta }{x}^{\ast }}{||x^{\ast }||}$ - относительная погрешность;

4. Вектор невязки - вектор, показывающий насколько найденное решение СЛАУ отклоняется от точного решения.

$$r=b-A{x}^{\ast }$$

Число обусловленности - коэффициент возможного возрастания относительной погрешности решения, вызванное погрешностью задания правой части.

Пусть $\bar{x}$ - точное решение системы $A\bar{x}=b$, а $x^{\ast }$ - решение системы. Тогда верны следующие оценки:

$$\delta (x^{\ast })\le {\nu }_{\delta }(\delta (A^{\ast })+\delta (b^{\ast }))$$

, где

$${\nu }_{\delta }=||A||\ast ||A^{-1}||\delta (A^{\ast })=\frac{||A-A^{-1}||}{||A||}|$$

$\nu (A)=cond(A)=||A^{-1}||\ast ||A||$ - стандартное число обусловленности.

Матрица плохо обусловлена, если $cond(A)>>1$. Следовательно, тогда существует решение, обладающее черезвычайно высокой чувствительностью к малым погрешностям входного данного $b$.

Прямые методы

Прямой метод - метод, который позволяет получить решение после выполнения конечного числа элементарных операций.

1. Метод Гаусса (схема единственного деления)

• Функция: gauss_single_division(A, b)

• A - Матрица левой части;

• b - Вектор правой части;

Трудоемкость метода - $\frac{2}{3}{m}^3$

• Прямой ход - матрица $A$ преобразуется к треугольному виду ($m-1$ - шагов).
• Обратный ход - вычисляются значения неизвестных, начиная с последнего уравнения ($m^2$ - шагов).

Условие применимости - схема единственного деления не может быть реализована, если один из главных элементов равен нулю.

Описание метода:

from mathmod.linear_systems import gauss_single_division

x = gauss_single_division(A, b)

---

Пример:

$$\begin{array}{rl}2x_1+x_2-x_3& =1,\\ 4x_1+3x_2-x_3& =7,\\ 8x_1+7x_2+3x_3& =25.\end{array}$$

В матричной форме система записывается так:

$$\left[\begin{array}{cccc}2& 1& -1& 1\\ 4& 3& -1& 7\\ 8& 7& 3& 25\end{array}\right]$$

Прямой ход

Шаг 1: Приведение первого столбца

Делим первую строку на ведущий элемент $a_{11}=2$:

$$\left[\begin{array}{cccc}1& 0.5& -0.5& 0.5\\ 4& 3& -1& 7\\ 8& 7& 3& 25\end{array}\right]$$

Обнуляем элементы ниже ведущего элемента в первом столбце, используя формулу:

$$a_{ij}\leftarrow a_{ij}-{\mu }_{ik}\cdot a_{kj},b_i\leftarrow b_i-{\mu }_{ik}\cdot b_k,$$

, где ${\mu }_{ik}=\frac{a_{ik}}{a_{kk}}$

Для второй строки:

$${\mu }_{21}=4,a_{2j}\leftarrow a_{2j}-4\cdot a_{1j}$$

Для третьей строки:

$${\mu }_{31}=8,a_{3j}\leftarrow a_{3j}-8\cdot a_{1j}$$

Получаем:

$$\left[\begin{array}{cccc}1& 0.5& -0.5& 0.5\\ 0& 1& 1& 5\\ 0& 3& 7& 21\end{array}\right]$$

Шаг 2: Приведение второго столбца

Делим вторую строку на ведущий элемент $a_{22}=1$ (он уже равен $1$):

$$\left[\begin{array}{cccc}1& 0.5& -0.5& 0.5\\ 0& 1& 1& 5\\ 0& 3& 7& 21\end{array}\right]$$

Обнуляем элементы ниже ведущего элемента во втором столбце:

$${\mu }_{32}=3,a_{3j}\leftarrow a_{3j}-3\cdot a_{2j}$$

Получаем:

$$\left[\begin{array}{cccc}1& 0.5& -0.5& 0.5\\ 0& 1& 1& 5\\ 0& 0& 4& 6\end{array}\right]$$

Обратный ход

Решаем треугольную систему методом подстановки.

Шаг 1: Найдем $x_3$:

$$x_3=\frac{6}{4}=1.5$$

Шаг 2: Найдем $x_2$:

$$x_2=5-1\cdot x_3=5-1\cdot 1.5=3.5$$

Шаг 3: Найдем $x_1$:

$$x_1=0.5-0.5\cdot x_2-0.5\cdot x_3=0.5-0.5\cdot 3.5+0.5\cdot 1.5=-0.5$$

Решение системы:

$$x_1=-0.5,x_2=3.5,x_3=1.5$$

---

2. Метод Гаусса (схема частичного выбора)

• Функция: gauss_partial_pivot(a, b)

Трудоемкость метода - $\frac{2}{3}{m}^3$

Описание метода:

Отличие от схемы единственного деления заключается в том, что на $k$-м шаге исключение в качестве главного элемента выбирают максимальный по модулю коэффициент $a_{i_{k}k}$.

Вычислительная устойчивость:

Гарантия ограниченности роста элементов матрицы делает схему частиного выбора вычислительно устойчивой. Становится справедлива оценка погрешности:

$$\delta (x^{\ast })\lesssim f(m)con{d}_E(A){ϵ}_M$$

, где:

$x^{\ast }$ - вычисленное ЭВМ решение системы. $\delta (x^{\ast })=\frac{||x-x^{\ast }|{|}_2}{||x|{|}_2}$ - относительная погрешность;

$con{d}_E(A)=||A|{|}_E||A^{-1}|{|}_E$ - числро обусловленности;

${ϵ}_M$ - машинный эпсилон;

$f(m)=C(m)\varphi (m)$, где $C(m)$ - некоторая медленно растущая функция, зависящая от порядка стистемы;

$\varphi (m)$ - коэффициент роста.

Далее исключение неизвестного $x_k$ производят, как в схеме единственного деления.

from mathmod.linear_systems import gauss_partial_pivot

x = gauss_partial_pivot(A,b)

---

Пример:

Рассмотрим систему:

$$A=\left[\begin{array}{ccc}2& -9& 5\\ 0& 3.5& -10\\ 0& 0.0001& 3\end{array}\right]b=\left[\begin{array}{c}-4\\ -6.5\\ 3.0001\end{array}\right]$$

---

3. Метод Холецкого (LLT-разложение)

• Для решения СЛАУ с симметрично положительно определённой матрицей.
• Функция: cholecky(A, b)

Трудоемкость метода - $\frac{1}{3}{m}^3$

Условия применимости - требуется, чтобы диагональные элементы $l_{ii}$ матрицы $L$ были положительными.

Достоинства метода:

• Гарантированная устойчивость;
• Требует вдвое меньше вычислительных затрат по сравнению с методом Гаусса;
• Позволяет экономично использовать память ЭВМ при записи исходных данных и результато вычислений за счет симметричности матрицы A.

Описание метода:

$A$ - симетрично положительно определенная матроица ($A=A^T$ и $\mathrm{\forall }x\ne 0$ скалярное произведение $(Ax,x)>0$ ).

$L$ - нижнетреугольная матрица. $L^T$ - транспонированная.

$$A=L{L}^T$$

, где

$$L{L}^T=\left[\begin{array}{ccccc}{l}_{11}& 0& 0& \dots & 0\\ l_{21}& l_{22}& 0& \dots & 0\\ l_{31}& l_{32}& l_{33}& \dots & 0\\ \dots & \dots & \dots & \dots & \dots \\ l_{m1}& l_{m2}& l_{m3}& \dots & l_{mm}\end{array}\right]\cdot \left[\begin{array}{ccccc}{l}_{11}& l_{21}& l_{31}& \dots & l_{m1}\\ 0& l_{22}& l_{23}& \dots & l_{m2}\\ 0& 0& l_{33}& \dots & l_{m3}\\ \dots & \dots & \dots & \dots & \dots \\ 0& 0& 0& \dots & l_{mm}\end{array}\right]=\left[\begin{array}{ccccc}{a}_{11}& a_{12}& a_{13}& \dots & a_{1m}\\ a_{21}& a_{22}& a_{23}& \dots & a_{2m}\\ a_{31}& a_{32}& a_{33}& \dots & a_{3m}\\ \dots & \dots & \dots & \dots & \dots \\ a_{m1}& a_{m2}& a_{m3}& \dots & a_{mm}\end{array}\right]$$

Отсюда:

$$l_{k1}^2+l_{k2}^2+\cdots +l_{kk}^2=a_{kk}$$

$$l_{kk}=\sqrt{a_{kk}-\sum _{j=1}^{k-1}{l}_{k,j}^2}k=2\dots m$$

$$l_{ik}=\frac{a_{i,k}-\sum _{j-1}^{k-1}{l}_{i,j}\cdot l_{l,k}}{l_{k,k}}i=k+1,\dots m$$

Если разложение получено, то решение системы:

$$Ly=b{L}^{T}x=y$$

from mathmod.linear_systems import cholecky

x = cholecky(A, b)

---

Пример:

Рассмотрим систему:

$$A=\left[\begin{array}{ccc}6.25& -1& 0.5\\ -1& 5& 2.12\\ 0.5& 2.12& 3.6\end{array}\right]b=\left[\begin{array}{c}7.5\\ -8.68\\ -0.24\end{array}\right]$$

$l_{11}=\sqrt{(}{a}_{11})=\sqrt{6.25}=2.5l_{21}=\frac{a_{21}}{l_{11}}=\frac{-1}{2.5}=-0.4$

$l_{31}=\frac{a_{31}}{l_{11}}=\frac{0.5}{2.5}=0.2l_{22}=\sqrt{a_{22}-l_{21}^2}=\sqrt{5-0.16}=2.2$

$l_{32}=a_{32}-l_{31}{l}_{21}=(2.12-\frac{0.2\cdot (-0.4)}{2.2})=1$

$l_{33}=\sqrt{a_{22}-l_{31}^2-l_{32}^2}=\sqrt{3.6-{0.2}^2-1^2}=1.6$

Матрица $L$:

$$L=\left[\begin{array}{ccc}2.25& 0& 0\\ -0.4& 2.2& 0\\ 0.2& 1& 1.6\end{array}\right]$$

Решение состоит из 2-х шагов:

1. Решаем $Ly=b$ для $y$ методом прямой подстановки.
2. Решаем $L^{T}x=y$ для $x$ методом обратной подстановки.

Решение:

1. Прямой ход для $y$:

$$Ly=b$$

$$\left[\begin{array}{ccc}2.5& 0& 0\\ -0.4& 2.2& 0\\ 0.2& 1& 1.6\end{array}\right]\cdot \left[\begin{array}{c}{y}_1\\ y_2\\ y_3\end{array}\right]=\left[\begin{array}{c}7.5\\ -8.68\\ -0.24\end{array}\right]$$

• Рассчитываем $y$:

$$\begin{array}{rl}2.5\cdot y_1& =7.5,\\ -0.4\cdot y_1+2.2\cdot y_2& =-8.68,\\ 0.2\cdot y_1+y_2+1.6\cdot y_3& =-0.24\end{array}$$

• Решая, получаем:

$$y_1=3,y_2=-3.4,y_3=1.6$$

2. Обратный ход для $x$:

$$L^{T}x=y$$

$$\left[\begin{array}{ccc}2.5& -0.4& 0.2\\ 0& 2.2& 1\\ 0& 0& 1.6\end{array}\right]\cdot \left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]=\left[\begin{array}{c}3\\ -3.4\\ 1.6\end{array}\right]$$

• Рассчитывем $x$:

$$\begin{array}{rl}2.5\cdot x_1-0.4\cdot x_2+0.2\cdot x_3& =3,\\ 2.2\cdot x_2+x_3& =-3.4,\\ 1.6\cdot x_3& =1.6.\end{array}$$

• Решая, получаем:

$$x_1=0.8,x_2=-2,x_3=1$$

---

4. Метод LU-разложения

• Функция: lu(A, b)

Трудоемкость метода - $\frac{2}{3}{m}^3+m^2$

Теорема о возможности применения $LU$ - разложения - если все главные миноры матрицы $A$ отличны от нуля, то существуют единственная нижняя треугольная матрица $L$ и верхняя треугольная матрица $U$ такие, что:

$$A=LU$$

, где

$$L=\left[\begin{array}{ccccc}1& 0& 0& \dots & 0\\ {\mu }_{21}& 1& 0& \dots & 0\\ {\mu }_{31}& {\mu }_{32}& 1& \dots & 0\\ \dots & \dots & \dots & \dots & \dots \\ {\mu }_{m1}& {\mu }_{m2}& {\mu }_{m3}& \dots & 1\end{array}\right]$$

Описание метода:

from mathmod.linear_systems import lu_solve

x = lu_solve(A, b)

---

Пример:

Рассмотрим систмему:

$$A=\left[\begin{array}{ccc}2& -1& -2\\ -4& 6& 3\\ -4& -2& 8\end{array}\right]b=\left[\begin{array}{c}-5\\ 6\\ 8\end{array}\right]$$

Мы хотим представить её в виде произведения:

$$A=LU$$

, где:

$$L=\left[\begin{array}{ccc}1& 0& 0\\ {\mu }_{21}& 1& 0\\ {\mu }_{31}& {\mu }_{32}& 1\end{array}\right],U=\left[\begin{array}{ccc}{u}_{11}& u_{12}& u_{13}\\ 0& u_{22}& u_{23}\\ 0& 0& u_{33}\end{array}\right].$$

Шаг 1: Прямой ход (разложение)

1. Выбираем первый элемент матрицы $A$ как ведущий:

$$u_{11}=2$$

Элементы верхней матрицы $U$:

$$u_{12}=-1,u_{13}=-2.$$

2. Вычисляем коэффициенты для матрицы $L$:

$${\mu }_{21}=\frac{a_{21}}{u_{11}}=\frac{-4}{2}=-2,{\mu }_{31}=\frac{a_{31}}{u_{11}}=\frac{-4}{2}=-2.$$

3. Обновляем элементы второй строки:

$$u_{22}=a_{22}-{\mu }_{21}\cdot u_{12}=6-(-2)\cdot (-1)=4,$$

$$u_{23}=a_{23}-{\mu }_{21}\cdot u_{13}=3-(-2)\cdot (-2)=-1.$$

4. Обновляем элементы третьей строки:

$$u_{32}=a_{32}-{\mu }_{21}\cdot u_{12}=-2-(-2)\cdot (-1)=-4$$

$$u_{33}=a_{33}-{\mu }_{21}\cdot u_{13}=8-(-2)\cdot (-2)=4$$

5. Промежуточные результаты:

$$L=\left[\begin{array}{ccc}1& 0& 0\\ -2& 1& 0\\ -2& {\mu }_{32}& 1\end{array}\right],U=\left[\begin{array}{ccc}2& -1& -2\\ 0& 4& -1\\ 0& -4& 4\end{array}\right],$$

6. Вычисляем ${\mu }_{32}$:

$${\mu }_{32}=\frac{u_{32}}{u_{22}}=-1$$

7. Обновляем элементы третьей строки:

$$u_{33}=u_{33}-{\mu }_{32}\cdot u_{23}=4-(-1)\cdot (-1)=3$$

8. Итоговые матрицы

$$L=\left[\begin{array}{ccc}1& 0& 0\\ -2& 1& 0\\ -2& -1& 1\end{array}\right],U=\left[\begin{array}{ccc}2& -1& -2\\ 0& 4& -1\\ 0& 0& 3\end{array}\right].$$

Шаг 2: решение системы

Для решения системы $Ax=b$, где:

$$A=LU,$$

Решение состоит из 2-х шагов:

1. Решаем $Ly=b$ для $y$ методом прямой подстановки.
2. Решаем $Ux=y$ для $x$ методом обратной подстановки.

Решение:

1. Прямой ход для $y$:

$$Ly=b$$

$$\left[\begin{array}{ccc}1& 0& 0\\ -2& 1& 0\\ -2& -1& 1\end{array}\right]\cdot \left[\begin{array}{c}{y}_1\\ y_2\\ y_3\end{array}\right]=\left[\begin{array}{c}-5\\ 6\\ 8\end{array}\right]$$

• Рассчитывем $y$:

$$\begin{array}{rl}{y}_1& =-5,\\ -2\cdot y_1+y_2& =6,\\ -2\cdot y_1-1\cdot y_2+y_3& =8.\end{array}$$

• Решая, получаем:

$$y_1=-5,y_2=-4,y_3=-6$$

2. Обратный ход для $x$:

$$Ux=y$$

$$\left[\begin{array}{ccc}2& -1& -2\\ 0& 4& -1\\ 0& 0& 3\end{array}\right]\cdot \left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]=\left[\begin{array}{c}-5\\ -4\\ -6\end{array}\right]$$

• Рассчитывем $x$:

$$\begin{array}{rl}2\cdot x_1-x_2-2\cdot x_3& =-5,\\ 4\cdot x_2-x_3& =-4,\\ 3\cdot x_3& =-6.\end{array}$$

• Решая, получаем:

$$x_1=-5.25,x_2=-1.5,x_3=-2$$

---

5. Метод прогонки

• Функция: three_diag(A,b)

Трудоемкость метода - $8m$

Условие применимости - коэффициенты системы удовлетворяют условиям диагонального преобладания:

$$|b_k|\ge |a_k|+|c_k||b_k|>|a_k|$$

Тогда:

$${\gamma }_i=b_i+a_i{\alpha }_{i-1}\ne 0|{\alpha }_i|\le 1\mathrm{\forall }i=1,2,\dots m$$

Описание метода:

$A$ - терхдиагональная матрица:

$$A=\left[\begin{array}{cccccc}{b}_1& c_1& 0& \dots & \dots & 0\\ a_2& b_2& c_2& \dots & \dots & 0\\ 0& a_3& b_3& c_3& \dots & 0\\ \dots & \dots & \dots & \dots & \dots & \dots \\ \dots & \dots & \dots & a_{m-1}& b_{m-1}& c_{m-1}\\ 0& 0& 0& \dots & a_m& b_m\end{array}\right]b=\left[\begin{array}{c}{d}_1\\ d_2\\ d_3\\ ⋮\\ b_{m-1}\\ b_m\end{array}\right]$$

• Прямой ход (прямая прогонка) - вычисление прогоночных коэффициентов

$${\alpha }_i=-\frac{c_i}{{\gamma }_i}{\beta }_i=\frac{d_i-{\alpha }_i{\beta }_{i-1}}{{\gamma }_i}{\gamma }_i=b_i+a_i{\alpha }_{i-1}$$

• Обратная прогонка (обратная прогонка) - вычисление значения незвестных. Сначала $x_m={\beta }_m$. Затем значения осталных неизветных по формуле:

$$x_i={\alpha }_i{x}_{i+1}+{\beta }_{i}i=m-1,m-2,\dots ,1$$

from mathmod.linear_systems import three_diag

x = three_diag(A,b)

---

Пример:

Рассмотрим систмему:

$$A=\left[\begin{array}{cccc}5& -1& 0& 0\\ 2& 4.6& -1& 0\\ 0& 2& 3.6& -0.8\\ 0& 0& 3& 4.4\end{array}\right]b=\left[\begin{array}{c}2\\ 3.3\\ 2.6\\ 7.2\end{array}\right]$$

Прямой ход

$${\gamma }_1=b_1=5,{\alpha }_1=-c_1/{\gamma }_1=0.2,{\beta }_1=d_1/{\gamma }_1=2.0/5=0.4,$$

$${\gamma }_2=b_2+a_2{\alpha }_1=4.6+2.0\cdot 0.2=5,{\alpha }_2=-c_2/{\gamma }_2=1/5=0.2,$$

$${\beta }_2=(d_2-a_2{\beta }_1)/{\gamma }_2=(3.3-2.0\cdot 0.4)/5=0.5,$$

$${\gamma }_3=b_3+a_3{\alpha }_2=3.6+2.0\cdot 0.2=4,{\alpha }_3=-c_3/{\gamma }_3=0.8/4=0.2,$$

$${\beta }_3=(d_3-a_3{\beta }_2)/{\gamma }_3=(2.6-2.0\cdot 0.5)/4=0.4,$$

$${\gamma }_4=b_4+a_4{\alpha }_3=4.4+3.0\cdot 0.2=5,{\beta }_4=(d_4-a_4{\beta }_3)/{\gamma }_4=(7.2-3.0\cdot 0.4)/5=1.2.$$

Обратный ход

$$x_4={\beta }_4=1.2,$$

$$x_3={\alpha }_3{z}_4+{\beta }_3=0.2\cdot 1.2+0.4=0.64,$$

$$x_2={\alpha }_2{z}_3+{\beta }_2=0.2\cdot 0.64+0.5=0.628,$$

$$x_1={\alpha }_1{z}_2+{\beta }_1=0.2\cdot 0.628+0.4=0.5256.$$

Итак, получаем решение:

$$x_1=0.5256,x_2=0.628,x_3=0.64,x_4=1.2.$$

---

Итерационные методы

Итерационный метод - метод, который строит последовательность приближений к решению;

1. Метод Якоби

• Функция: jacobi(A, b, epsilon=1e-6, norma=1)

  • A - Матрица коэффициентов (n x n);
  • b - Вектор правой части;
  • epsilon - Заданная точность (по умолчанию ${10}^{-6}$);
  • norma - Норма, по которой считается критерий окончания (например, 1, 2, np.inf).

Теорема о сходимости:

Пусть выполнено условие:

$$||B||<1$$

Тогда решение системы $\bar{x}$ существует и единственно при произволном приближении $x^{(0)}$ МПИ сходится и справедлива оценка погрешности (априорная оценка):

$$||x^{(n)}-\bar{x}||\le ||B|{|}^n||x^{(0)}-\bar{x}||$$

Апостериорная оценка:

$$||x^{(n)}-\bar{x}||\le \frac{||B||}{1-||B||}||x^{(n)}-x^{(n-1)}||$$

Критерий окончания:

$$||x^{(n)}-x^{(n-1)}||\le {ϵ}_1$$

, где:

$${ϵ}_1=\frac{||B||}{1-||B||}ϵ$$

Более простой критерий окончания:

$$||x^{(n)}-x^{(n-1)}||\le ϵ$$

Описание метода:

$$\begin{array}{r}{x}_1^{k+1}=b_{11}{x}_1^k+b_{12}{x}_2^k+b_{13}{x}_3^k+\cdots +b_{1m}{x}_m^k+c_1,\\ x_2^{k+1}=b_{21}{x}_1^k+b_{22}{x}_2^k+b_{23}{x}_3^k+\cdots +b_{2m}{x}_m^k+c_2,\\ x_3^{k+1}=b_{31}{x}_1^k+b_{32}{x}_2^k+b_{33}{x}_3^k+\cdots +b_{3m}{x}_m^k+c_3,\\ \dots \\ x_m^{k+1}=b_{m1}{x}_1^k+b_{m2}{x}_2^k+b_{m3}{x}_3^k+\cdots +b_{mm}{x}_m^k+c_m,\end{array}$$

from mathmod.linear_systems import jacobi

x, iteration_count = jacobi(A, b, epsilon=1e-6, norma=1)

---

Пример:

---

2. Метод Гаусса-Зейделя

• Итерационный метод для решения СЛАУ с диагонально преобладающей матрицей.

• Функция: gauss_zeydel(A, b, epsilon=1e-6, norma=1)

  • A - Матрица коэффициентов (n x n);
  • b - Вектор правой части;
  • epsilon - Заданная точность (по умолчанию ${10}^{-6}$);
  • norma - Норма, по которой считается критерий окончания (например, 1, 2, np.inf).

3. Метод релаксации

• Функция: relaxation_method(A, b, epsilon=1e-6, omega=1, norma=1)

  • A - Матрица коэффициентов (n x n);
  • b - Вектор правой части;
  • epsilon - Заданная точность (по умолчанию ${10}^{-6}$);
  • omega - Параметр релаксации (по умолчанию 1.0 — метод Зейделя);
  • norma - Норма, по которой считается критерий окончания (например, 1, 2, np.inf).

from mathmod.linear_systems import relaxation_method

x, iteration_count = relaxation_method(A, b, epsilon=1e-6, omega=1, norma=1)

---

Приближение функций

Постановка задачи

Известны значения некоторой функции $f(x)$ только на множестве дискретных точек $x_0,x_1,\dots ,x_n$, но само аналитическое выражение для функции неизвестно. Заменим функцию $f(x)$ некоторой известной и достаточно легко вычисляемой функцией $\mathrm{\Phi }(x)$ такой, что $\mathrm{\Phi }(x)\approx f(x)$. Подобный процесс замены неизвестной функции некоторой близкой функцией называется аппроксимацией, а функция $\mathrm{\Phi }(x)$ называется аппроксимирующей функцией.

Для аппроксимации функций широко используются классы функций вида:

$${\mathrm{\Phi }}_m(x)=a_0{\varphi }_0(x)+a_1{\varphi }_1(x)+\dots +a_m{\varphi }_m(x),$$

являющиеся линейными комбинациями фиксированного набора базисных функций ${\varphi }_0(x),{\varphi }_1(x),\dots ,{\varphi }_m(x)$.

Функцию ${\varphi }_m(x)$ называют обобщенным многочленом по системе функций ${\varphi }_0(x),{\varphi }_1(x),\dots ,{\varphi }_m(x)$

Число $m$ - степенью многочлена.

Существуют два основных подхода в аппроксимации функций:

1. Пусть точки $f(x_i),i=0,1,\dots ,n$ получены в результате достаточно точных измерений или вычислений, т.е. есть основания считать их лишенными ошибок. Тогда следует выбирать аппроксимирующую функцию $\varphi (x)$ такой, чтобы она совпадала со значениями исходной функции в заданных точках. Геометрически это означает, что кривая $\varphi (x)$ проходит через точки $(x_i,f(x_i))$ плоскости. Такой метод приближения называется интерполяцией.

2. Если точки $f(x_i),i=0,1,\dots ,n$ содержат ошибки (данные экспериментов, статистические данные и т.п.), то функция $\varphi (x)$ выбирается из условия минимума некоторого функционала, обеспечивающего сглаживание ошибок. Такой прием называется аппроксимацией функции «в среднем». Геометрически это будет означать, что кривая $\varphi (x)$ будет занимать некоторое «среднее» положение, не обязательно совпадая с исходными точками $(x_i,f(x_i))$ плоскости.

Постановка задачи интерполяции

Пусть в точках $x_0,x_1,\dots ,x_n$, расположенных на отрезке $[a,b]$ и попарно различных.

Тогда задача итерполяции состоит в построении функции $g(x)$, удовлетворяющей условию:

$$g(x)=y_i(i=0,1,\dots ,n)$$

Интерполяция - способ приближения функции $f(x)$ путем построения функии $g(x)$, график которой проходит через точки $(x_i,y_i)$.

Экстраполяция - способ приближения функции $f(x)$ в точке $x<x_{min}$ или $x>x_{max}$.

$[x_{min},x_{max}]$ - минимальный и максимальный из узлов интерполяции.

Теорема о существовании и единственности интерполяционного многочлена:

Если $m=n$, то решение задачи интерполяции обобщенным многочленом ($\Phi_m(x) = a_0\phi_0(x) + a_1\phi_1(x) + \ldots + a_m\phi_m(x)$) существует и единственно при любом наборе данных $y_0,y_1,\dots ,y_n$, тогда и только тогда, когда системы функций ${\varphi }_0(x),{\varphi }_1(x),\dots ,{\varphi }_n(x)$ являются линейно независимыми в точках $x_0,x_0,\dots ,x_n$.

Полиномиальная интерполяция. Многочлен Лагранжа

$$L_n(x)=\sum _{i=0}^n{y}_i{l}_{nj}(x)$$

, где:

$$l_{ij}(x)=\prod _{\genfrac{}{}{0ex}{}{k=1}{k\ne j}}^n\frac{(x-x_k)}{(x_j-x_i)}=\frac{(x-x_0)(x-x_1)\dots (x-x_{j-1})(x-x_{j+1})\dots (x-x_n)}{(x_j-x_0)(x_j-x_1)\dots (x_j-x_{j-1})(x_j-x_{j+1})\dots (x_j-x_n)}$$

Оценка погрешности:

$$\underset{[a,b]}{max}|f(x)-P_n(x)|\le \frac{M_{n+1}}{(n+1)!}\underset{[a,b]}{max}|{\omega }_{n+1}(x)|$$

${\omega }_{n+1}(x)=(x-x_0)(x-x_1)\dots (x-x_n)$

$M_{n+1}=\underset{[a,b]}{max}|f^{n+1}(x)|$

или

$$\underset{[x_0,x_n]}{max}|f(x)-P_n(x)|\le \frac{M_{n+1}}{4(n+1)}{h}_{max}^{n+1}$$

$h_{max}=\underset{1\le i\le n}{max}{h}_i$

Эта формула позволяет утвержать, что для достаточно гладкой функции $f$ при фиксированной степени интерполяционного многочлена погрешность интерполяции на отрезке $[x_0,x_n]$ при $h_{max}\to 0$ стрепится к нулю не медленее, чем некоторая величина, пропорциональная $h_{max}^{n+1}$.

Итерполяция многочленом $n$ имеет $(n+1)$-й порядок точности относительно $h_{max}$.

---

Пример:

Пусть даны точки:

• $x_0=0$, $y_0=1$
• $x_1=1$, $y_1=2$
• $x_2=2$, $y_2=4$

Построим интерполяционный многочлен Лагранжа:

$$L_2(x)=y_0\cdot \frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}+y_1\cdot \frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}+y_2\cdot \frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}$$

Подставляя значения:

$$L_2(x)=1\cdot \frac{(x-1)(x-2)}{(0-1)(0-2)}+2\cdot \frac{(x-0)(x-2)}{(1-0)(1-2)}+4\cdot \frac{(x-0)(x-1)}{(2-0)(2-1)}$$

$$L_2(x)=1\cdot \frac{(x-1)(x-2)}{2}-2\cdot \frac{x(x-2)}{1}+2\cdot \frac{x(x-1)}{2}$$

---

Интерполяционный многочлен Ньютона с конечными разностями

Пусть интерполируемая функция задана таблицей с постоянными шагом $h$, т.е. $x_i=x_0+ih,i=0,1\dots n$. Тогда введя безразмерную величину $t=\frac{x-x_0}{h}$, можно записать многочлен Ньютона:

$$P_n(x)=P_n(x+ht)=y_0+\frac{\mathrm{\Delta }{y}_0}{1!}t+\frac{{\mathrm{\Delta }}^2{y}_0}{2!}t(t-1)+\frac{{\mathrm{\Delta }}^3{y}_0}{3!}t(t-1)(t-2)+\dots +\frac{{\mathrm{\Delta }}^n{y}_0}{n!}t(t-1)(t-2)\dots (t-n+1)$$

---

Метод наименьших квадратов

Постановка задачи

Пусть даны точки $x_0,x_1,\dots ,x_n$, и известны значения исходной фукнции $f_i=f(x_i),i=0,1,\dots n$. Требуется найти многочлен $P_m$ заданной степени $m(m=n)$ такой, чтобы величина среднеквадратического отклонения:

$$\sigma (P_m,f)=\sqrt{\frac{1}{1+n}\sum _{i=0}^n(P_m(x_i)-f_i{)}^2}=\sqrt{\frac{1}{1+n}\sum _{i=0}^n(\sum _{j=0}^n{a}_j{x}_i^j-f_i{)}^2}$$

была минимальной

Нормальная система:

$$\sum _{j=0}^m{a}_j\sum _{i=0}^n{x}_i^{k+j}=\sum _{i=0}^n{f}_i{x}_i^{k}k=0,1,\dots m$$

$$s_k=\sum _{i=0}^n{x}_i^k{b}_k=\sum _{i=0}^n{f}_i{x}_i^k$$

$$\begin{array}{r}(n+1)a_0+(\sum _{i=0}^n{x}_i)a_1+(\sum _{i=0}^n{x}_i^2)a_2=\sum _{i=0}^n{y}_i,\\ (\sum _{i=0}^n{x}_i)a_0+(\sum _{i=0}^n{x}_i^2)a_1+(\sum _{i=0}^n{x}_i^3)a_2=\sum _{i=0}^n{y}_i{x}_i,\\ (\sum _{i=0}^n{x}_i^2)a_0+(\sum _{i=0}^n{x}_i^3)a_1+(\sum _{i=0}^n{x}_i^4)a_2=\sum _{i=0}^n{y}_i{x}_i^2,\end{array}$$

---

Установка

Для использования библиотеки склонируйте репозиторий и установите необходимые зависимости:

git clone https://github.com/BaranovSerV/mathmod.git
cd mathmod
pip install -r requirements.txt