Математика

Численные методы

Обзор вычислительных методов и выбор алгоритма по типу задачи.

Как выбрать подходящий материал

  1. Сначала определите тип исходных данных: уравнение, таблица значений, интеграл, линейная система или задача Коши. От этого зависит семейство алгоритмов.
  2. Затем проверьте условия сходимости и критерий остановки. Малое изменение двух итераций само по себе не доказывает близость к точному ответу.

5 формул

Интерполяция

Восстановление промежуточных значений по данным узлам.

44 формулы

Формулы темы

Составная формула трапеций для вычисления интеграла

Составная формула трапеций приближает определённый интеграл суммой площадей под линейными отрезками между соседними узлами. Концевые значения получают половинный вес.

$T_n=h\left[\frac{f(a)+f(b)}2+\sum_{i=1}^{n-1}f(a+ih)\right],\quad h=\frac{b-a}{n}$

Формула Симпсона: интеграл по параболическим дугам

Составная формула Симпсона объединяет соседние промежутки попарно и заменяет функцию параболой через три узла. Для равномерной сетки число промежутков должно быть чётным.

$S_n=\frac h3\left[f_0+f_n+4\!\sum_{\substack{i=1\\i\text{ нечётное}}}^{n-1}f_i+2\!\sum_{\substack{i=2\\i\text{ чётное}}}^{n-2}f_i\right]$

Составная формула средних прямоугольников

Составная формула средних прямоугольников заменяет интеграл суммой площадей с высотами в серединах равных отрезков. Для гладкой функции погрешность имеет второй порядок по шагу, хотя значения на концах не используются.

$M_n=h\sum_{i=0}^{n-1}f\!\left(a+\left(i+\tfrac12\right)h\right),\quad h=\frac{b-a}{n}$

Квадратура Гаусса—Лежандра с двумя узлами

Двухузловая квадратура Гаусса—Лежандра вычисляет функцию в двух специально выбранных внутренних точках. На конечном отрезке с единичным весом она точно интегрирует любой многочлен степени не выше трёх.

$G_2=\frac{b-a}{2}\left[f\!\left(m-\frac{b-a}{2\sqrt3}\right)+f\!\left(m+\frac{b-a}{2\sqrt3}\right)\right],\quad m=\frac{a+b}{2}$

Метод Ромберга: таблица уточнений интеграла

Метод Ромберга строит таблицу из составных трапеций с последовательным удвоением числа отрезков. Экстраполяция устраняет ведущие чётные степени шага и повышает точность без смены уже вычисленных узлов.

$R_{k,0}=T_{2^k},\quad R_{k,j}=R_{k,j-1}+\frac{R_{k,j-1}-R_{k-1,j-1}}{4^j-1},\quad 1\le j\le k$

Адаптивный Симпсон: локальная оценка ошибки

Адаптивный метод Симпсона сравнивает одну параболическую оценку с двумя оценками на половинах. Разность даёт локальный индикатор ошибки, по которому алгоритм уточняет только те участки, где не достигнут допуск.

$S_1=S(a,b),\ S_2=S(a,m)+S(m,b),\quad \varepsilon_{\mathrm{est}}=\frac{|S_2-S_1|}{15},\quad I\approx S_2+\frac{S_2-S_1}{15}$

Интеграл методом Монте-Карло: среднее и ошибка

Метод Монте-Карло оценивает интеграл как длину отрезка, умноженную на среднее значений функции в случайных равномерных точках. Разброс оценки измеряют стандартной ошибкой, которая при независимой выборке убывает как N⁻¹ᐟ².

$\hat I_N=\frac{b-a}{N}\sum_{i=1}^N f(X_i),\quad X_i\overset{\mathrm{iid}}\sim U(a,b),\quad\widehat{\mathrm{SE}}=(b-a)\frac{s_f}{\sqrt N}$

Линейная интерполяция: значение между двумя узлами

Линейная интерполяция оценивает значение между двумя заданными точками по соединяющей их прямой. Доля положения внутри интервала задаёт веса концов; оценка учитывает наклон, но не кривизну исходного графика.

$P(x)=y_0+\frac{y_1-y_0}{x_1-x_0}(x-x_0)=(1-t)y_0+ty_1,\quad t=\frac{x-x_0}{x_1-x_0}$

Интерполяция Лагранжа: полином через заданные точки

Формула Лагранжа строит полином степени не выше n по n+1 различным узлам. Каждый базисный множитель выбирает значение своего узла и обращает вклады остальных узлов в ноль.

$P_n(x)=\sum_{k=0}^{n}y_kL_k(x),\qquad L_k(x)=\prod_{\substack{0\leq j\leq n\\j\ne k}}\frac{x-x_j}{x_k-x_j}$

Интерполяция Ньютона: таблица разделённых разностей

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

$P_n(x)=f[x_0]+\sum_{k=1}^{n}f[x_0,\ldots,x_k]\prod_{j=0}^{k-1}(x-x_j),\quad f[x_i,\ldots,x_j]=\frac{f[x_{i+1},\ldots,x_j]-f[x_i,\ldots,x_{j-1}]}{x_j-x_i}$

Кубическая интерполяция Эрмита по значениям и наклонам

Кубический полином Эрмита согласует значения и первые производные на двух концах интервала. Четыре заданных условия определяют единственную кубическую или более простую кривую между узлами.

$H(x)=h_{00}(s)y_0+h_{10}(s)\ell m_0+h_{01}(s)y_1+h_{11}(s)\ell m_1,\quad s=\frac{x-x_0}{\ell},\quad \ell=x_1-x_0$

Барицентрическая интерполяция: веса и вычисление полинома

Барицентрическая форма вычисляет полином Лагранжа через заранее рассчитанные веса узлов. Она удобна для множества запросов по одной таблице; при точном совпадении аргумента с узлом нужно сразу вернуть узловое значение.

$P(x)=\frac{\sum_{k=0}^{n}\dfrac{w_ky_k}{x-x_k}}{\sum_{k=0}^{n}\dfrac{w_k}{x-x_k}},\qquad w_k=\frac{1}{\prod_{j\ne k}(x_k-x_j)},\quad x\ne x_k$

Правая разность: производная по двум соседним значениям

Правая разность оценивает производную по значениям функции в текущем и следующем узле. Это удобная односторонняя схема у левой границы таблицы, но её ошибка обычно линейна по шагу.

$f'(x)\approx D_h^+f(x)=\frac{f(x+h)-f(x)}{h},\qquad D_h^+f-f'=O(h)$

Центральная разность для первой производной

Центральная разность оценивает первую производную по двум узлам на одинаковом расстоянии слева и справа. Симметрия устраняет линейную по шагу ошибку и обычно даёт второй порядок точности.

$f'(x)\approx D_hf(x)=\frac{f(x+h)-f(x-h)}{2h},\qquad D_hf-f'=O(h^2)$

Вторая производная по трём равноотстоящим узлам

Трёхточечная центральная схема оценивает вторую производную по отклонению центрального значения от среднего соседей. Она измеряет изменение наклона, а не сам наклон, и имеет второй порядок по шагу.

$f''(x)\approx Q_hf(x)=\frac{f(x+h)-2f(x)+f(x-h)}{h^2},\qquad Q_hf-f''=O(h^2)$

Пятиточечная формула первой производной

Пятиточечная центральная формула оценивает первую производную с четвёртым порядком по шагу. Для неё нужны четыре значения вокруг центра: вес самого центрального значения равен нулю.

$f'(x)\approx P_hf=\frac{f(x-2h)-8f(x-h)+8f(x+h)-f(x+2h)}{12h},\qquad P_hf-f'=O(h^4)$

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

Экстраполяция Ричардсона объединяет две центральные оценки производной на шагах h и h/2. При гладкой функции она сокращает ведущую ошибку второго порядка и даёт оценку четвёртого порядка.

$R_h=\frac{4D_{h/2}-D_h}{3},\qquad D_h=\frac{f(x+h)-f(x-h)}{2h},\qquad R_h-f'(x)=O(h^4)$

Метод Якоби для системы линейных уравнений

Метод Якоби строит новое приближение к решению Ax=b, одновременно пересчитывая все неизвестные из значений предыдущего шага. Он удобен для разреженных систем и хорошо показывает, зачем проверять сходимость, а не только выполнять заданное число итераций.

$x_i^{(k+1)}=\frac{b_i-\sum_{j\ne i}a_{ij}x_j^{(k)}}{a_{ii}},\quad i=1,\ldots,n$

Метод Гаусса—Зейделя для системы линейных уравнений

Метод Гаусса—Зейделя решает Ax=b последовательными проходами: новая компонента сразу участвует в вычислении следующих. Это обычно ускоряет простую итерацию, но делает результат одного прохода зависимым от порядка переменных.

$x_i^{(k+1)}=\frac{b_i-\sum_{j<i}a_{ij}x_j^{(k+1)}-\sum_{j>i}a_{ij}x_j^{(k)}}{a_{ii}}$

Метод SOR: последовательная сверхрелаксация

Метод SOR добавляет к проходу Гаусса—Зейделя параметр релаксации ω. При 0<ω<2 он дозирует поправку; ω=1 возвращает обычный метод, а удачный выбор ω>1 способен ускорить сходимость.

$x_i^{(k+1)}=(1-\omega)x_i^{(k)}+\frac{\omega}{a_{ii}}\left(b_i-\sum_{j<i}a_{ij}x_j^{(k+1)}-\sum_{j>i}a_{ij}x_j^{(k)}\right)$

Метод прогонки для трёхдиагональной системы

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

$c_i x_{i-1}+a_i x_i+b_i x_{i+1}=d_i;\quad \alpha_i=-\frac{b_i}{a_i+c_i\alpha_{i-1}},\;\beta_i=\frac{d_i-c_i\beta_{i-1}}{a_i+c_i\alpha_{i-1}}$

Метод Ричардсона: стационарная итерация для Ax=b

Метод Ричардсона прибавляет к текущему приближению невязку, умноженную на параметр τ. Для симметричной положительно определённой матрицы разумный τ превращает эту простую формулу в сходящуюся стационарную итерацию.

$x^{(k+1)}=x^{(k)}+\tau(b-Ax^{(k)}),\qquad r^{(k)}=b-Ax^{(k)}$

Явный метод Эйлера: решение задачи Коши по шагам

Явный метод Эйлера приближённо решает уравнение y′ = f(t,y), продолжая решение по наклону в начале каждого шага. При гладком решении глобальная ошибка имеет первый порядок по шагу.

$y_{n+1}=y_n+h f(t_n,y_n),\qquad t_{n+1}=t_n+h$

Неявный метод Эйлера: шаг и устойчивость затухания

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

$y_{n+1}=y_n+h f(t_{n+1},y_{n+1})$

Явный метод средней точки для задачи Коши

Явный метод средней точки делает пробный полушаг по Эйлеру, оценивает производную в середине и использует её на полном шаге. Это двухстадийная схема Рунге — Кутты второго порядка.

$k_1=f(t_n,y_n),\quad k_2=f(t_n+h/2,y_n+hk_1/2),\quad y_{n+1}=y_n+hk_2$

Метод Хойна: явное усреднение двух наклонов

В этой двухстадийной схеме сначала делают прогноз Эйлера на полный шаг, затем усредняют начальный и прогнозный конечный наклоны. Метод имеет второй порядок, но названия его вариантов различаются в литературе.

$k_1=f(t_n,y_n),\quad k_2=f(t_n+h,y_n+hk_1),\quad y_{n+1}=y_n+\frac h2(k_1+k_2)$

Классический метод Рунге — Кутты четвёртого порядка

Классический RK4 вычисляет четыре наклона за шаг и объединяет их с весами 1:2:2:1. Для гладкой задачи Коши он даёт четвёртый порядок общей ошибки без явного вычисления высших производных.

$\begin{aligned}k_1&=f(t_n,y_n),&k_2&=f(t_n+h/2,y_n+hk_1/2),\\k_3&=f(t_n+h/2,y_n+hk_2/2),&k_4&=f(t_n+h,y_n+hk_3),\\y_{n+1}&=y_n+\frac h6(k_1+2k_2+2k_3+k_4).\end{aligned}$

Метод Тейлора второго порядка для задачи Коши

Метод Тейлора второго порядка добавляет к шагу Эйлера член с полной производной правой части. Для скалярного уравнения y′=f(t,y) вторая производная решения равна fₜ+fᵧf, вычисленной в начале шага.

$y_{n+1}=y_n+h f_n+\frac{h^2}{2}\left(f_t+f_y f\right)_n$

Неявная средняя точка: шаг и сохранение структуры

Неявная средняя точка вычисляет производную в среднем старого и нового состояний. Новый уровень входит в правую часть, поэтому шаг находят решением уравнения. Схема второго порядка обратима по времени при согласованном решении.

$y_{n+1}=y_n+h f\!\left(t_n+\frac h2,\frac{y_n+y_{n+1}}2\right)$

Неявная трапециевидная схема Кранка — Николсон для ОДУ

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

$y_{n+1}=y_n+\frac h2[f(t_n,y_n)+f(t_{n+1},y_{n+1})]$

Метод золотого сечения: сужение интервала минимума

Золотое сечение находит минимум унимодальной функции без производных. Два внутренних значения позволяют отбросить часть интервала, а удачное отношение длин сохраняет одно вычисление для следующего шага.

$\tau=\frac{\sqrt5-1}{2},\quad u=a+(1-\tau)(b-a),\quad v=a+\tau(b-a),\quad L_n=\tau^nL_0$

Тернарный поиск минимума непрерывной унимодальной функции

Тернарный поиск сравнивает функцию в двух точках, делящих интервал на трети. При унимодальности одно сравнение позволяет удалить крайнюю треть и сохранить интервал, содержащий минимум.

$u=a+\frac{b-a}{3},\quad v=b-\frac{b-a}{3},\quad [a,b]\leftarrow\begin{cases}\left[u,b\right],&f(u)>f(v),\\\left[a,v\right],&f(u)\le f(v),\end{cases}\quad L_n=(2/3)^nL_0$

Параболическая интерполяция для поиска минимума функции

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

$A=\frac{(f_3-f_2)/(x_3-x_2)-(f_2-f_1)/(x_2-x_1)}{x_3-x_1},\quad B=\frac{f_2-f_1}{x_2-x_1}-A(x_1+x_2),\quad x_p=-\frac{B}{2A}$

Градиентный спуск с постоянным шагом: условие убывания

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

$x_{k+1}=x_k-\alpha\nabla f(x_k),\quad 0<\alpha<2/L,\quad f(x_{k+1})\le f(x_k)-\alpha(1-L\alpha/2)\|\nabla f(x_k)\|^2$

Правило Армихо: обратный поиск длины шага спуска

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

$f(x+\alpha p)\le f(x)+c\alpha\nabla f(x)^Tp,\quad \nabla f(x)^Tp<0,\quad \alpha\leftarrow\beta\alpha,\quad 0<c<1,\ 0<\beta<1$

Метод Ньютона для минимума: шаг через матрицу Гессе

Метод Ньютона выбирает шаг по квадратичной модели целевой функции: матрица Гессе учитывает кривизну, а градиент задаёт линейную часть. Решение системы даёт направление, но длину шага иногда приходится уменьшать.

$H_k d_k=-g_k,\qquad x_{k+1}=x_k+\alpha_k d_k$

Покоординатный спуск для выпуклой квадратичной функции

Покоординатный спуск минимизирует функцию по одной переменной при фиксированных остальных. Для положительно определённой квадратичной цели точное обновление получается делением компоненты градиента на диагональ матрицы.

$x_i^{+}=x_i-\frac{(Ax-b)_i}{A_{ii}},\quad f(x)=\tfrac12x^TAx-b^Tx$

Проекционный градиент: минимум с ограничениями на переменные

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

$x_{k+1}=\Pi_C(x_k-\alpha\nabla f(x_k))$

Метод Франка — Вульфа: линейный подшаг на симплексе

Метод Франка — Вульфа заменяет проекцию решением линейной задачи на допустимой области. Новая точка лежит на отрезке между старой точкой и решением подзадачи; на симплексе линейный минимум достигается в одной из его вершин.

$s_k\in\arg\min_{s\in C}\nabla f(x_k)^Ts,\quad x_{k+1}=(1-\gamma_k)x_k+\gamma_k s_k$

Метод тяжёлого шарика: инерционный шаг для квадратичной цели

Метод тяжёлого шарика добавляет к градиентному шагу долю предыдущего перемещения. На положительно определённой квадратичной функции сходимость можно проверять по корням скалярной рекурсии для каждого собственного направления.

$x_{k+1}=x_k-\alpha\nabla f(x_k)+\beta(x_k-x_{k-1})$