Содержание

§ 5. Интерполяция на регулярной сетке

Виды сеток

Регулярной сеткой называется сетка, каждый элемент которой является отрезком (в случае одного измерения), параллелограммом (в случае двух измерений) или параллелепипедом (в случае трех).

Нерегулярной треугольной сеткой называется сетка, каждый элемент которой является треугольником (в случае двух измерений) или тетраэдром (в случае трех измерений).

Нерегулярной сеткой называется сетка, в которой нет какой-либо структуры. Иногда такие сетки называют облаками точек.

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

Математическая постановка задачи

Интерполяция — это поиск значений функции в промежуточных точках по известным значениям функции в окружающих точках. Математически это записывается как

yi=f(xi),i=1...n,yi=F(xi),i=1...n.y_i = f(x_i), \qquad i=1...n, \\ y_i = F(x_i), \qquad i=1...n.

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

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

Понятия интерполяции и сглаживания не являются математически строгими. Соответствующие методы схожи, однако находятся на противоположных концах спектра: интерполянт всегда проходит через точки, в которых задана функция, а сглаживающая функция проходит по некоторой усредненной с учетом окружающих известных точек кривой. Общим термином для этих двух методов в одномерном пространстве является приближение с помощью кривых. Если же мы хотим найти значения функции в точках, которые не окружены точками, в которых значения функции известны, то такие методы называются экстраполяцией.

Линейная интерполяция

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

f(x)=f(x0)+f(x1)f(x0)x1x0(xx0).f(x)=f(x0)(1xx0x1x0)+f(x1)xx0x1x0.f(x) = f(x_0) + \frac{f(x_1)-f(x_0)}{x_1-x_0} (x-x_0). \\ f(x) = f(x_0)\left(1 - \frac{x-x_0}{x_1-x_0}\right) + f(x_1)\frac{x-x_0}{x_1-x_0}.

В двухмерном случае формула запишется как

i=01j=01f(xi,yj)k=12{tk,если (k=1i=0)(k=2j=0),1tk,во всех других случаях,t1=xx0x1x0,t2=yy0y1y0.\sum\limits_{i=0}^{1} \sum\limits_{j=0}^{1} f(x_i,y_j) \prod\limits_{k=1}^{2} \begin{cases} t_k, & \text{если }(k=1 \wedge i=0) \vee (k=2 \wedge j=0), \\ 1-t_k, & \text{во всех других случаях}, \end{cases} \\ t_1 = \frac{x-x_0}{x_1-x_0}, \\ t_2 = \frac{y-y_0}{y_1-y_0}.

В nn-мерном случае формула запишется как

i=01f(xi)k=1n{tk,ik=0,1tk,ik0,t=xx0x1x0.\sum\limits_{\vec{i}=\vec{0}}^{\vec{1}} f(\vec{x}_{\vec{i}}) \prod\limits_{k=1}^{n} \begin{cases} \vec{t}_k, & \vec{i}_k=0, \\ 1-\vec{t}_k, & \vec{i}_k \neq 0, \end{cases} \\ \vec{t}= \frac{\vec{x}-\vec{x}_{\vec{0}}}{\vec{x}_{\vec{1}}-\vec{x}_{\vec{0}}}.

Здесь векторные индексы означают выбор дискретных точек, а скалярные индексы означают выбор компоненты вектора, 0\vec{0} — вектор нулей, 1\vec{1} — вектор единиц.

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

Интерполяция полиномами Лагранжа

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

yi=f(xi),i=1...n,pn(x)=a0+a1x+a2x2++anxn,pn(xi)=yi,[1x0x02x0n1x1x12x1n1xnxn2xnn][a0a1an]=[y0y1yn]y_i = f(x_i), \qquad i=1...n, \\ p_n(x) = a_0 + a_1 x + a_2 x^2 + \ldots + a_n x^n, \\ p_n(x_i) = y_i, \\ \left[\begin{array}{lllll} 1 & x_0 & x_0^2 & \ldots & x_0^n \\ 1 & x_1 & x_1^2 & \ldots & x_1^n \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & x_n & x_n^2 & \ldots & x_n^n \\ \end{array} \right] \left[\begin{array}{l}a_0\\a_1\\\vdots\\a_n\end{array} \right] = \left[\begin{array}{l}y_0\\y_1\\\vdots\\y_n\end{array} \right]

Решая эту систему уравнений, получаем следующую формулу для полинома:

p(x)=i=0n(0jn,jixxjxixj)yi.p(x) = \sum\limits_{i=0}^{n} \left(\prod\limits_{0 \leq j \leq n,j \neq i} \frac{x-x_j}{x_i-x_j} \right) y_i.

Эта форма записи полинома называется полиномом Лагранжа. При n=1n=1 эта формула сводится к формуле линейной интерполяции. Таким образом, линейная интерполяция является частным случаем полиномиальной, когда степень полинома равна единице.

Интерполяция полиномами Ньютона

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

nj(x)=i=0j(xxi),j=1...n,n_j(x) = \prod\limits_{i=0}^j \left(x-x_i\right), \qquad j=1...n,f(xj)=i=0jaini(xj),j=1...n.f(x_j) = \sum\limits_{i=0}^j a_i n_i(x_j), \qquad j=1...n.

Это приводит к системе уравнений

[101x1x01x2x0(x2x0)(x2x1)1xkx0j=0n1(xkxj)][a0a1an]=[y0y1yn]\left[\begin{array}{lllll} 1 & & \ldots & & 0 \\ 1 & x_1-x_0 & & & \\ 1 & x_2-x_0 & (x_2-x_0)(x_2-x_1) & & \vdots \\ \vdots & \vdots & & \ddots & \\ 1 & x_k-x_0 & \ldots & & \prod\limits_{j=0}^{n-1}(x_k-x_j) \\ \end{array}\right] \left[\begin{array}{l}a_0\\a_1\\\vdots\\a_n\end{array} \right] = \left[\begin{array}{l}y_0\\y_1\\\vdots\\y_n\end{array} \right]

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

Интерполяция сплайнами

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

Кубический сплайн (общий случай)

Запишем систему уравнений для одномерного кубического сплайна (сплайна третьего порядка). Нам известно значение функции в n+1n+1 точках f(x1),,f(xn+1)f(x_1), \ldots, f(x_{n+1}). Каждый интервал интерполируется полиномом третьей степени pi(t)=ai+bit+cit2+dit3p_i(t) = a_i + b_i t + c_i t^2 + d_i t^3, где t[0,1]t \in [0,1] — параметр, а i=1,,ni=1,\ldots,n. На концах интервалов полином запишется как

pi(0)=f(xi)=ai,pi(1)=f(xi+1)=ai+bi+ci+di,p_i(0) = f(x_i) = a_i, \\ p_i(1) = f(x_{i+1}) = a_i+b_i+c_i+d_i,

а его производная запишется как

pi(0)=f(xi)=bi,pi(1)=f(xi+1)=bi+2ci+3di.p_i'(0) = f'(x_i) = b_i, \\ p_i'(1) = f'(x_{i+1}) = b_i + 2 c_i + 3 d_i.

Решая последние четыре уравнения относительно коэффициентов, получим

ai=f(xi),bi=f(xi),ci=2(f(xi+1)f(xi))2f(xi)f(xi+1),di=2(f(xi)f(xi+1))+f(xi)+f(xi+1).a_i = f(x_i), \\ b_i = f'(x_i), \\ c_i = 2(f(x_{i+1}) - f(x_i)) - 2f'(x_i) - f'(x_{i+1}), \\ d_i = 2(f(x_i) - f(x_{i+1})) + f'(x_i) + f'(x_{i+1}).

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

Приравняем значения функции, а также ее первые и вторые производные на граничных точках соседних интервалов:

pi1(1)=pi(0),pi1(1)=pi(0),pi1(1)=pi(0).p_{i-1}(1) = p_i(0), \\ p_{i-1}'(1) = p_i'(0), \\ p_{i-1}''(1) = p_i''(0).

На левой и правой границах области определения функции обнулим вторые производные:

p1(0)=f(x1),pn(1)=f(xn+1),p1(0)=pn(1)=0.p_1(0) = f(x_1), p_n(1) = f(x_{n+1}), \\ p_1''(0) = p_n''(1) = 0.

Из полученных уравнений составляем систему

[2114114114114][f(x1)f(xn+1)]=[3(f(x2)f(x1))3(f(x3)f(x1))3(f(x4)f(x2))3(f(xn1)f(xn3))3(f(xn)f(xn2))3(f(xn)f(xn1))].\left[\begin{array}{llllll} 2 & 1 \\ 1 & 4 & 1 \\ & 1 & 4 & 1 \\ & & \ddots & \ddots & \ddots \\ & & & 1 & 4 & 1 \\ & & & & 1 & 4 \\ \end{array}\right] \left[\begin{array}{l}f'(x_1)\\\vdots\\f'(x_{n+1})\end{array}\right] = \left[\begin{array}{l} 3(f(x_2)-f(x_1)) \\ 3(f(x_3)-f(x_1)) \\ 3(f(x_4)-f(x_2)) \\ \vdots \\ 3(f(x_{n-1})-f(x_{n-3})) \\ 3(f(x_n)-f(x_{n-2})) \\ 3(f(x_n)-f(x_{n-1})) \\ \end{array}\right].

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

Кубический сплайн Эрмита

Эти сплайны не используют полином Эрмита, а только названы в честь этого ученого. Здесь производные аппроксимируются конечными разностями (подробный разговор о которых пойдет в последующих разделах, а здесь представлены готовые формулы). Полином для каждого промежутка записывается как

pi(x)=h00(t)f(xi)+h10(t)(xi+1xi)mi+h01(t)f(xi+1)+h11(t)(xi+1xi)mi+1,t=xxixi+1xi,h00=(1+2t)(1t)2,h10=t(1t)2,h01=t2(32t),h11=t2(t1).p_i(x) = h_{00}(t)f(x_i) + h_{10}(t) (x_{i+1}-x_i) m_i + h_{01}(t) f(x_{i+1}) + h_{11}(t) (x_{i+1}-x_i) m_{i+1}, \\ t = \frac{x - x_i}{x_{i+1}-x_i}, \\ h_{00} = (1 + 2t)(1-t)^2, \\ h_{10} = t(1-t)^2, \\ h_{01} = t^2(3-2t), \\ h_{11} = t^2(t-1).

Здесь hh — это базисные функции, mim_i — аппроксимации производных, которые записываются одним из следующих способов.

mi=12(f(xi+1)f(xi)xi+1xi+f(xi)f(xi1)xixi1) или mi=tf(xi+1)f(xi1)xi+1xi1,m_i = \frac{1}{2} \left( \frac{f(x_{i+1}) - f(x_i)}{x_{i+1}-x_i} + \frac{f(x_{i}) - f(x_{i-1})}{x_{i}-x_{i-1}} \right) \text{ или } m_i = t \frac{f(x_{i+1}) - f(x_{i-1})}{x_{i+1}-x_{i-1}},

где t[0,1]t \in [0,1] — параметр. Первая формула является средним значением производных двух соседних интервалов, записанных в виде конечных разностей. Вторую формулу называют сплайнами Катмалла—Рома.

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

Задания

Линейная интерполяция1 балл

Аппроксимируйте функцию sin(x)\sin(x) с помощью линейной интерполяции и с помощью полинома Лагранжа второй степени в точках с шагом 2π/102\pi/10. Уменьшая шаг, найдите такую его величину, при которой точность линейной интерполяции в середине промежутка между двумя точками сопоставима с точностью интерполяции полиномом Лагранжа второй степени с исходным шагом.

Экстремумы дискретно заданной функции2 балла

Найдите все точки минимума и максимума функции sin(x)\sin(x), заданной дискретно в точках с шагом 2π/102\pi/10. (Задайте точки таким образом, чтобы на них не попадал максимум и минимум.) Для этого найдите все точки, которые меньше (больше) соседних. В найденных точках аппроксимируйте функцию полиномом Лагранжа второй степени (параболой) и определите по формуле параболы координату точки экстремума. Сравните полученные значения с реальными.

Точки пересечения сложных фигур3 балла

Найдите точки пересечения прямоугольника и синусоиды (см. картинку справа). Прямоугольник повернут на произвольный угол. Синусоида описывается функцией sin(x)\sin(x), заданной дискретно в точках с шагом 2π/102\pi/10. Вычислите координаты точек пересечения синусоиды и прямоугольника с помощью интерполяции Лагранжа второй степени и метода бисекции.

Видео

2023

  • 00:00 Виды сеток.
  • 12:20 Интерполяция, экстраполяция, сглаживание.
  • 19:13 Линейная интерполяция.
  • 13:10 Полиномы Лагранжа.
  • 28:13 Полиномы Ньютона.
  • 30:49 Кубический сплайн (общий случай).
  • 40:26 Кубический сплайн Эрмита.
  • 49:25 Многомерный случай.
  • 53:45 Задания.
Запись лекции 24.03.2023.

2022

  • 00:00 Виды сеток.
  • 09:10 Интерполяция и экстраполяция.
  • 10:40 Линейная интерполяция.
  • 13:10 Интерполяция и сглаживание.
  • 17:40 Полиномы Лагранжа.
  • 21:20 Сплайны.
  • 44:12 Задания.
Запись лекции 17.03.2022.

Доска

2023

interpolation-regular-grid-2023 thumbnail.
Запись доски 24.03.2023.

2022

interpolation-regular-grid thumbnail.
Запись доски 17.03.2022.