Полиномиальная интерполяция

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

Учитывая набор из n + 1 точек данных $(x_{0},y_{0}),\ldots ,(x_{n},y_{n})$, без двух $x_{j}$ то же, полиномиальная функция $p(x)=a_{0}+a_{1}x+\cdots +a_{n}x^{n}$ Говорят, что интерполирует данные, если $p(x_{j})=y_{j}$ для каждого $j\in \{0,1,\dotsc ,n\}$.

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

Приложения

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

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

В численном анализе полиномиальная интерполяция необходима для выполнения субквадратного умножения и возведения в квадрат, например умножения Карацубы и умножения Тума – Кука , где интерполяция через точки на полиноме произведения дает конкретный требуемый продукт. Например, учитывая a = f ( x ) = a ₀ x ⁰ + 1 _​ х ¹ + ··· и б знак равно г ( Икс ) знак равно б ₀ Икс ⁰ + б ₁ х ¹ + ···, произведение ab представляет собой конкретное значение W ( x ) = f ( x ) g ( x ). Можно легко найти точки вдоль W ( x ) при малых значениях x , и интерполяция на основе этих точек даст члены W ( x ) и конкретное произведение ab . Как сформулировано в умножении Карацубы, этот метод существенно быстрее, чем квадратичное умножение, даже для входных данных скромного размера, особенно на параллельном оборудовании.

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

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

Для любого $n+1$ двумерные точки данных $(x_{0},y_{0}),\dotsc ,(x_{n},y_{n})\in \mathbb {R} ^{2}$, где нет двух $x_{j}$ одинаковы, существует единственный полином $p(x)$ степени максимум $n$ который интерполирует эти точки, т.е. $p(x_{0})=y_{0},\ldots ,p(x_{n})=y_{n}$.

Эквивалентно, для фиксированного выбора узлов интерполяции $x_{j}$полиномиальная интерполяция определяет линейную биекцию $L_{n}$ между ( n +1)-кортежами значений действительных чисел $(y_{0},\ldots ,y_{n})\in \mathbb {R} ^{n+1}$ и векторное пространство $P(n)$ действительных многочленов степени не выше n : $$L_{n}:\mathbb {R} ^{n+1}{\stackrel {\sim }{\longrightarrow }}\,P(n).$$

Это разновидность теоремы о неразрешимости . Теорема справедлива и для любого бесконечного поля вместо действительных чисел. $\mathbb {R}$, например рациональные или комплексные числа.

Первое доказательство

Рассмотрим базисные функции Лагранжа $L_{1}(x),\ldots ,L_{n}(x)$ предоставлено: $$L_{j}(x)=\prod _{i\neq j}{\frac {x-x_{i}}{x_{j}-x_{i}}}={\frac {(x-x_{0})\cdots (x-x_{j-1})(x-x_{j+1})\cdots (x-x_{n})}{(x_{j}-x_{0})\cdots (x_{j}-x_{j-1})(x_{j}-x_{j+1})\cdots (x_{j}-x_{n})}}.$$

Обратите внимание, что $L_{j}(x)$ является полиномом степени $n$, и у нас есть $L_{j}(x_{k})=0$ для каждого $j\neq k$, пока $L_{k}(x_{k})=1$. Отсюда следует, что линейная комбинация: $$p(x)=\sum _{j=0}^{n}y_{j}L_{j}(x)$$имеет $p(x_{k})=\sum _{j}y_{j}\,L_{j}(x_{k})=y_{k}$, так $p(x)$ представляет собой интерполяционный полином степени $n$.

Для доказательства единственности предположим, что существует другой интерполяционный полином $q(x)$ степени максимум $n$, так что $p(x_{k})=q(x_{k})$ для всех $k=0,\dotsc ,n$. Затем $p(x)-q(x)$ является полиномом степени не более $n$ который имеет $n+1$ отдельные нули (т. $x_{k}$). Но ненулевой полином степени не выше $n$ может иметь максимум $n$ нули, так $p(x)-q(x)$ должен быть нулевым полиномом, т.е. $p(x)=q(x)$.

Второе доказательство

Запишите интерполяционный полином в виде

$$p(x)=a_{n}x^{n}+a_{n-1}x^{n-1}+\cdots +a_{2}x^{2}+a_{1}x+a_{0}.$$ ( 1 )

Подставив это в интерполяционные уравнения $p(x_{j})=y_{j}$, получим систему линейных уравнений с коэффициентами $a_{j}$, который читается в матрично-векторной форме как следующее умножение : $${\begin{bmatrix}x_{0}^{n}&x_{0}^{n-1}&x_{0}^{n-2}&\ldots &x_{0}&1\\x_{1}^{n}&x_{1}^{n-1}&x_{1}^{n-2}&\ldots &x_{1}&1\\\vdots &\vdots &\vdots &&\vdots &\vdots \\x_{n}^{n}&x_{n}^{n-1}&x_{n}^{n-2}&\ldots &x_{n}&1\end{bmatrix}}{\begin{bmatrix}a_{n}\\a_{n-1}\\\vdots \\a_{0}\end{bmatrix}}={\begin{bmatrix}y_{0}\\y_{1}\\\vdots \\y_{n}\end{bmatrix}}.$$

Интерполянт $p(x)$ соответствует решению $A=(a_{n},\ldots ,a_{0})$ приведенного выше матричного уравнения $X\cdot A=Y$. Матрица X слева — это матрица Вандермонда , определитель которой, как известно, равен $\textstyle \det(X)=\prod _{1\leq i<j\leq n}(x_{j}-x_{i}),$ что не равно нулю, поскольку узлы $x_{j}$ все различны. Это гарантирует, что матрица обратима и уравнение имеет единственное решение. $A=X^{-1}\cdot Y$; то есть, $p(x)$ существует и является единственным.

Следствие

Если $f(x)$ является полиномом степени не более $n$, то интерполяционный полином $f(x)$ в $n+1$ отдельные точки это $f(x)$ сам.

Построение интерполяционного полинома

Красные точки обозначают точки данных ( x _k , y _k ) , а синяя кривая показывает полином интерполяции.

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

Основная статья: Полином Лагранжа

Мы можем сразу записать полином через полиномы Лагранжа как: $${\begin{aligned}p(x)&={\frac {(x-x_{1})(x-x_{2})\cdots (x-x_{n})}{(x_{0}-x_{1})(x_{0}-x_{2})\cdots (x_{0}-x_{n})}}y_{0}\\[4pt]&+{\frac {(x-x_{0})(x-x_{2})\cdots (x-x_{n})}{(x_{1}-x_{0})(x_{1}-x_{2})\cdots (x_{1}-x_{n})}}y_{1}\\[4pt]&+\cdots \\[4pt]&+{\frac {(x-x_{0})(x-x_{1})\cdots (x-x_{n-1})}{(x_{n}-x_{0})(x_{n}-x_{1})\cdots (x_{n}-x_{n-1})}}y_{n}\\[7pt]&=\sum _{i=0}^{n}{\Biggl (}\prod _{\stackrel {\!0\,\leq \,j\,\leq \,n}{j\,\neq \,i}}{\frac {x-x_{j}}{x_{i}-x_{j}}}{\Biggr )}y_{i}=\sum _{i=0}^{n}{\frac {p(x)}{p'(x_{i})(x-x_{i})}}\,y_{i}\end{aligned}}$$Для матричных аргументов эта формула называется формулой Сильвестра , а полиномы Лагранжа с матричным знаком являются ковариантами Фробениуса .

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

См. Также: Полином Ньютона и разделенные разности.

Теорема

Для полинома $p_{n}$ степени меньше или равной n, которая интерполирует $f$ в узлах $x_{i}$ где $i=0,1,2,3,\cdots ,n$. Позволять $p_{n+1}$ — многочлен степени меньше или равной n+1, который интерполирует $f$ в узлах $x_{i}$ где $i=0,1,2,3,\cdots ,n,n+1$. Затем $p_{n+1}$ дается: $$p_{n+1}(x)=p_{n}(x)+a_{n+1}w_{n}(x)$$где ${\textstyle w_{n}(x):=\prod _{i=0}^{n}(x-x_{i})}$ также известный как базис Ньютона и ${\textstyle a_{n+1}:={f(x_{n+1})-p_{n}(x_{n+1}) \over w_{n}(x_{n+1})}}$.

Доказательство:

Это можно показать для случая, когда $i=0,1,2,3,\cdots ,n$: $$p_{n+1}(x_{i})=p_{n}(x_{i})+a_{n+1}\prod _{j=0}^{n}(x_{i}-x_{j})=p_{n}(x_{i})$$и когда $i=n+1$: $$p_{n+1}(x_{n+1})=p_{n}(x_{n+1})+{f(x_{n+1})-p_{n}(x_{n+1}) \over w_{n}(x_{n+1})}w_{n}(x_{n+1})=f(x_{n+1})$$В силу единственности интерполированных полиномов степени меньше $n+1$, ${\textstyle p_{n+1}(x)=p_{n}(x)+a_{n+1}w_{n}(x)}$ – требуемая полиномиальная интерполяция. Таким образом, функцию можно выразить как:

${\textstyle p_{n}(x)=a_{0}+a_{1}(x-x_{0})+a_{2}(x-x_{0})(x-x_{1})+\cdots +a_{n}(x-x_{0})\cdots (x-x_{n-1}).}$

Полиномиальные коэффициенты

Найти $a_{i}$, нам нужно решить нижнюю треугольную матрицу, образованную расстановкой ${\textstyle p_{n}(x_{i})=f(x_{i})=y_{i}}$ из приведенного выше уравнения в матричной форме:

${\begin{bmatrix}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 &\ldots &\prod _{j=0}^{n-1}(x_{n}-x_{j})\end{bmatrix}}{\begin{bmatrix}a_{0}\\\\\vdots \\\\a_{n}\end{bmatrix}}={\begin{bmatrix}y_{0}\\\\\vdots \\\\y_{n}\end{bmatrix}}$

Коэффициенты выводятся как

$a_{j}:=[y_{0},\ldots ,y_{j}]$

где

$[y_{0},\ldots ,y_{j}]$

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

showДоказательство

The first few coefficients can be calculated using the system of equations. The form of n-th coefficient is assumed for proof by mathematical induction. ${\begin{aligned}a_{0}&=y_{0}=[y_{0}]\\a_{1}&={y_{1}-y_{0} \over x_{1}-x_{0}}=[y_{0},y_{1}]\\\vdots \\a_{n}&=[y_{0},\cdots ,y_{n}]\quad {\text{(let)}}\\\end{aligned}}$ Let Q be polynomial interpolation of points $(x_{1},y_{1}),\ldots ,(x_{n},y_{n})$. Adding $(x_{0},y_{0})$ to the polynomial Q: $Q(x)+a'_{n}(x-x_{1})\cdot \ldots \cdot (x-x_{n})=P_{n}(x),$ where ${\textstyle a'_{n}(x_{0}-x_{1})\ldots (x_{0}-x_{n})=y_{0}-Q(x_{0})}$. By uniqueness of the interpolating polynomial of the points $(x_{0},y_{0}),\ldots ,(x_{n},y_{n})$, equating the coefficients of $x^{n-1}$ we get, ${\textstyle a'_{n}=[y_{0},\ldots ,y_{n}]}$. Hence the polynomial can be expressed as:$P_{n}(x)=Q(x)+[y_{0},\ldots ,y_{n}](x-x_{1})\cdot \ldots \cdot (x-x_{n}).$ Adding $(x_{n+1},y_{n+1})$ to the polynomial Q, it has to satisfiy: ${\textstyle [y_{1},\ldots ,y_{n+1}](x_{n+1}-x_{1})\cdot \ldots \cdot (x_{n+1}-x_{n})=y_{n+1}-Q(x_{n+1})}$ where the formula for ${\textstyle a_{n}}$ and interpolating polynomial are used.The ${\textstyle a_{n+1}}$ term for the polynomial ${\textstyle P_{n+1}}$ can be found by calculating:$${\begin{aligned}&[y_{0},\ldots ,y_{n+1}](x_{n+1}-x_{0})\cdot \ldots \cdot (x_{n+1}-x_{n})\\&={\frac {[y_{1},\ldots ,y_{n+1}]-[y_{0},\ldots ,y_{n}]}{x_{n+1}-x_{0}}}(x_{n+1}-x_{0})\cdot \ldots \cdot (x_{n+1}-x_{n})\\&=\left([y_{1},\ldots ,y_{n+1}]-[y_{0},\ldots ,y_{n}]\right)(x_{n+1}-x_{1})\cdot \ldots \cdot (x_{n+1}-x_{n})\\&=[y_{1},\ldots ,y_{n+1}](x_{n+1}-x_{1})\cdot \ldots \cdot (x_{n+1}-x_{n})-[y_{0},\ldots ,y_{n}](x_{n+1}-x_{1})\cdot \ldots \cdot (x_{n+1}-x_{n})\\&=(y_{n+1}-Q(x_{n+1}))-[y_{0},\ldots ,y_{n}](x_{n+1}-x_{1})\cdot \ldots \cdot (x_{n+1}-x_{n})\\&=y_{n+1}-(Q(x_{n+1})+[y_{0},\ldots ,y_{n}](x_{n+1}-x_{1})\cdot \ldots \cdot (x_{n+1}-x_{n}))\\&=y_{n+1}-P(x_{n+1}).\end{aligned}}$$which implies that $a_{n+1}={y_{n+1}-P_{n}(x_{n+1}) \over w_{n}(x_{n+1})}=[y_{0},\ldots ,y_{n+1}]$. Hence it is proved by principle of mathematical induction.

Форвардная формула Ньютона

Полином Ньютона можно выразить в упрощенной форме, если $x_{0},x_{1},\dots ,x_{k}$ расположены последовательно с одинаковым интервалом.

Если $x_{0},x_{1},\dots ,x_{k}$ расположены последовательно и на равном расстоянии от ${x}_{i}={x}_{0}+ih$ для i = 0, 1, ..., k и некоторой переменной x выражается как ${x}={x}_{0}+sh$, то разница $x-x_{i}$ можно записать как $(s-i)h$. Таким образом, полином Ньютона становится

${\begin{aligned}N(x)&=[y_{0}]+[y_{0},y_{1}]sh+\cdots +[y_{0},\ldots ,y_{k}]s(s-1)\cdots (s-k+1){h}^{k}\\&=\sum _{i=0}^{k}s(s-1)\cdots (s-i+1){h}^{i}[y_{0},\ldots ,y_{i}]\\&=\sum _{i=0}^{k}{s \choose i}i!{h}^{i}[y_{0},\ldots ,y_{i}].\end{aligned}}$

Поскольку связь между разделенными разностями и прямыми разницами определяется как: $$[y_{j},y_{j+1},\ldots ,y_{j+n}]={\frac {1}{n!h^{n}}}\Delta ^{(n)}y_{j},$$принимая $y_{i}=f(x_{i})$, если вместо этого было принято представление x в предыдущих разделах $x=x_{j}+sh$, формула прямой интерполяции Ньютона выражается как: $$f(x)\approx N(x)=N(x_{j}+sh)=\sum _{i=0}^{k}{s \choose i}\Delta ^{(i)}f(x_{j})$$что является интерполяцией всех точек после $x_{j}$. Он расширен как: $$f(x_{j}+sh)=f(x_{j})+{\frac {s}{1!}}\Delta f(x_{j})+{\frac {s(s-1)}{2!}}\Delta ^{2}f(x_{j})+{\frac {s(s-1)(s-2)}{3!}}\Delta ^{3}f(x_{j})+{\frac {s(s-1)(s-2)(s-3)}{4!}}\Delta ^{4}f(x_{j})+\cdots$$

Обратная формула Ньютона

Если узлы переупорядочены как ${x}_{k},{x}_{k-1},\dots ,{x}_{0}$, полином Ньютона становится

$N(x)=[y_{k}]+[{y}_{k},{y}_{k-1}](x-{x}_{k})+\cdots +[{y}_{k},\ldots ,{y}_{0}](x-{x}_{k})(x-{x}_{k-1})\cdots (x-{x}_{1}).$

Если ${x}_{k},\;{x}_{k-1},\;\dots ,\;{x}_{0}$ расположены на равном расстоянии от ${x}_{i}={x}_{k}-(k-i)h$ для i = 0, 1, ..., k и ${x}={x}_{k}+sh$, затем,

${\begin{aligned}N(x)&=[{y}_{k}]+[{y}_{k},{y}_{k-1}]sh+\cdots +[{y}_{k},\ldots ,{y}_{0}]s(s+1)\cdots (s+k-1){h}^{k}\\&=\sum _{i=0}^{k}{(-1)}^{i}{-s \choose i}i!{h}^{i}[{y}_{k},\ldots ,{y}_{k-i}].\end{aligned}}$

Поскольку связь между разделенными разностями и обратными разностями задается как: $$[{y}_{j},y_{j-1},\ldots ,{y}_{j-n}]={\frac {1}{n!h^{n}}}\nabla ^{(n)}y_{j},$$принимая $y_{i}=f(x_{i})$, если вместо этого было принято представление x в предыдущих разделах $x=x_{j}+sh$, формула обратной интерполяции Ньютона выражается как: $$f(x)\approx N(x)=N(x_{j}+sh)=\sum _{i=0}^{k}{(-1)}^{i}{-s \choose i}\nabla ^{(i)}f(x_{j}).$$что является интерполяцией всех точек перед $x_{j}$. Он расширен как: $$f(x_{j}+sh)=f(x_{j})+{\frac {s}{1!}}\nabla f(x_{j})+{\frac {s(s+1)}{2!}}\nabla ^{2}f(x_{j})+{\frac {s(s+1)(s+2)}{3!}}\nabla ^{3}f(x_{j})+{\frac {s(s+1)(s+2)(s+3)}{4!}}\nabla ^{4}f(x_{j})+\cdots$$

Ромбическая диаграмма

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

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

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

Коэффициенты выражаются по формуле: $$C(u+k,n)={\frac {(u+k)(u+k-1)\cdots (u+k-n+1)}{n!}}$$

Доказательство эквивалентности

Если путь идет от $\Delta ^{n-1}y_{s}$ к $\Delta ^{n+1}y_{s-1}$, оно может соединяться посредством трех промежуточных шагов: (а) через $\Delta ^{n}y_{s-1}$, (б) через ${\textstyle C(u-s,n)}$ или (c) через $\Delta ^{n}y_{s}$. Доказательство эквивалентности этих трех двухшаговых путей должно доказать, что все (n-шаговые) пути могут быть преобразованы с одинаковым началом и концом, и все они представляют собой одну и ту же формулу.

Путь (а):

$C(u-s,n)\Delta ^{n}y_{s-1}+C(u-s+1,n+1)\Delta ^{n+1}y_{s-1}$

Путь (б):

$C(u-s,n)\Delta ^{n}y_{s}+C(u-s,n+1)\Delta ^{n+1}y_{s-1}$

Путь (с):

$C(u-s,n){\frac {\Delta ^{n}y_{s-1}+\Delta ^{n}y_{s}}{2}}\quad +{\frac {C(u-s+1,n+1)+C(u-s,n+1)}{2}}\Delta ^{n+1}y_{s-1}$

Вычитая вклады от путей a и b:

${\begin{aligned}{\text{Path a - Path b}}=&C(u-s,n)(\Delta ^{n}y_{s-1}-\Delta ^{n}y_{s})+(C(u-s+1,n+1)-C(u-s,n-1))\Delta ^{n+1}y_{s-1}\\=&-C(u-s,n)\Delta ^{n+1}y_{s-1}+C(u-s,n){\frac {(u-s+1)-(u-s-n)}{n+1}}\Delta ^{n+1}y_{s-1}\\=&C(u-s,n)(-\Delta ^{n+1}y_{s-1}+\Delta ^{n+1}y_{s-1})=0\\\end{aligned}}$

Таким образом, вклад пути (а) или пути (б) одинаков. Поскольку путь (c) является средним значением путей (a) и (b), он также вносит в полином идентичную функцию. Таким образом, доказывается эквивалентность путей с одинаковыми начальной и конечной точками. Чтобы проверить, можно ли сдвинуть пути на разные значения в крайнем левом углу, достаточно сделать всего два шага пути: (a) $y_{s+1}$ к $y_{s}$ через $\Delta y_{s}$ или (b) коэффициент между $y_{s+1}$ и $y_{s}$, к $y_{s}$ через $\Delta y_{s}$ или (c) начиная с $y_{s}$.

Путь (а)

$y_{s+1}+C(u-s-1,1)\Delta y_{s}-C(u-s,1)\Delta y_{s}$

Путь (б)

${\frac {y_{s+1}+y_{s}}{2}}+{\frac {C(u-s-1,1)+C(u-s,1)}{2}}\Delta y_{s}-C(u-s,1)\Delta y_{s}$

Путь (с)

$y_{s}$

С $\Delta y_{s}=y_{s+1}-y_{s}$, подстановка в приведенные выше уравнения показывает, что все приведенные выше члены сводятся к $y_{s}$ и, следовательно, эквивалентны. Следовательно, эти пути можно трансформировать, чтобы они начинались с крайнего левого угла и заканчивались в общей точке.

Формула Ньютона

Взяв поперечный отрицательный наклон из $y_{0}$ к $\Delta ^{n}y_{0}$ дает интерполяционную формулу всех $n+1$ последовательно расположенные точки, что эквивалентно формуле прямой интерполяции Ньютона:

${\begin{aligned}y(s)&=y_{0}+C(s,1)\Delta y_{0}+C(s,2)\Delta ^{2}y_{0}+C(s,3)\Delta ^{3}y_{0}+\cdots \\&=y_{0}+s\Delta y_{0}+{\frac {s(s-1)}{2}}\Delta ^{2}y_{0}+{\frac {s(s-1)(s-2)}{3!}}\Delta ^{3}y_{0}+{\frac {s(s-1)(s-2)(s-3)}{4!}}\Delta ^{4}y_{0}+\cdots \end{aligned}}$

тогда как, взяв поперечный положительный наклон от $y_{n}$ к $\nabla ^{n}y_{n}=\Delta ^{n}y_{0}$, дает интерполяционную формулу всех $n+1$ последовательно расположенные точки, что эквивалентно формуле обратной интерполяции Ньютона:

${\begin{aligned}y(u)&=y_{k}+C(u-k,1)\Delta y_{k-1}+C(u-k+1,2)\Delta ^{2}y_{k-2}+C(u-k+2,3)\Delta ^{3}y_{k-3}+\cdots \\&=y_{k}+(u-k)\Delta y_{k-1}+{\frac {(u-k+1)(u-k)}{2}}\Delta ^{2}y_{k-2}+{\frac {(u-k+2)(u-k+1)(u-k)}{3!}}\Delta ^{3}y_{k-3}+\cdots \\y(k+s)&=y_{k}+(s)\nabla y_{k}+{\frac {(s+1)s}{2}}\nabla ^{2}y_{k}+{\frac {(s+2)(s+1)s}{3!}}\nabla ^{3}y_{k}+{\frac {(s+3)(s+2)(s+1)s}{4!}}\nabla ^{4}y_{k}+\cdots \\\end{aligned}}$

где $s=u-k$ — число, соответствующее введенному в интерполяции Ньютона.

Формула Гаусса

Пройдите зигзагообразной линией вправо, начиная с $y_{0}$ с отрицательным наклоном мы получаем формулу Гаусса:

$y(u)=y_{0}+u\Delta y_{0}+{\frac {u(u-1)}{2}}\Delta ^{2}y_{-1}+{\frac {(u+1)u\left(u-1\right)}{3!}}\Delta ^{3}y_{-1}+{\frac {(u+1)u\left(u-1\right)(u-2)}{4!}}\Delta ^{4}y_{-2}+\cdots$

тогда как начиная с $y_{0}$ с положительным наклоном мы получаем обратную формулу Гаусса:

$y(u)=y_{0}+u\Delta y_{-1}+{\frac {(u+1)u}{2}}\Delta ^{2}y_{-1}+{\frac {(u+1)u\left(u-1\right)}{3!}}\Delta ^{3}y_{-2}+{\frac {(u+2)(u+1)u\left(u-1\right)}{4!}}\Delta ^{4}y_{-2}+\cdots$

Формула Стирлинга

Пройдя горизонтальный путь вправо, начиная с $y_{0}$, получаем формулу Стирлинга:

${\begin{aligned}y(u)&=y_{0}+u{\frac {\Delta y_{0}+\Delta y_{-1}}{2}}+{\frac {C(u+1,2)+C(u,2)}{2}}\Delta ^{2}y_{-1}+C(u+1,3){\frac {\Delta ^{3}y_{-2}+\Delta ^{3}y_{-1}}{2}}+\cdots \\&=y_{0}+u{\frac {\Delta y_{0}+\Delta y_{-1}}{2}}+{\frac {u^{2}}{2}}\Delta ^{2}y_{-1}+{\frac {u(u^{2}-1)}{3!}}{\frac {\Delta ^{3}y_{-2}+\Delta ^{3}y_{-1}}{2}}+{\frac {u^{2}(u^{2}-1)}{4!}}\Delta ^{4}y_{-2}+\cdots \end{aligned}}$

Формула Стирлинга представляет собой среднее арифметическое прямой и обратной формул Гаусса.

Формула Бесселя

Пройдя горизонтальный путь вправо, начиная с фактора между $y_{0}$ и $y_{1}$, получаем формулу Стирлинга:

${\begin{aligned}y(u)&=1{\frac {y_{0}+y_{1}}{2}}+{\frac {C(u,1)+C(u-1,1)}{2}}\Delta y_{0}+C(u,2){\frac {\Delta ^{2}y_{-1}+\Delta ^{2}y_{0}}{2}}+\cdots \\&={\frac {y_{0}+y_{1}}{2}}+\left(u-{\frac {1}{2}}\right)\Delta y_{0}+{\frac {u(u-1)}{2}}{\frac {\Delta ^{2}y_{-1}+\Delta ^{2}y_{0}}{2}}+{\frac {\left(u-{\frac {1}{2}}\right)u\left(u-1\right)}{3!}}\Delta ^{3}y_{0}+{\frac {(u+1)u(u-1)(u-2)}{4!}}{\frac {\Delta ^{4}y_{-1}+\Delta ^{4}y_{-2}}{2}}+\cdots \\\end{aligned}}$

Алгоритмы Вандермонда

Матрица Вандермонда во втором доказательстве выше может иметь большое число обусловленности , вызывая большие ошибки при вычислении коэффициентов a _i, если система уравнений решается методом исключения Гаусса .

Поэтому несколько авторов предложили алгоритмы, которые используют структуру матрицы Вандермонда для вычисления численно устойчивых решений за O( n ² ) вместо операций O( n ³ ) требуется методом исключения Гаусса. Эти методы основаны на построении сначала интерполяции Ньютона полинома, а затем на преобразовании его в мономиальную форму .

Алгоритмы, не относящиеся к Вандермонду

Чтобы найти интерполяционный полином p ( x ) в векторном пространстве P ( n ) полиномов степени n , мы можем использовать обычный мономиальный базис для P ( n ) и инвертировать матрицу Вандермонда методом исключения Гаусса, что дает вычислительные затраты O ( н ³ ) операции. Чтобы улучшить этот алгоритм, более удобный базис для P ( n ) может упростить вычисление коэффициентов, которые затем необходимо перевести обратно в термины мономиального базиса .

Один из методов — записать интерполяционный полином в форме Ньютона (т.е. использовать базис Ньютона) и использовать метод разделенных разностей для построения коэффициентов, например алгоритм Невилла . Стоимость равна O( n ² ) операции. Более того, вам нужно выполнить дополнительную работу O( n ) только в том случае, если к набору данных добавляется дополнительная точка, тогда как для других методов вам придется переделать все вычисления.

Другой метод предпочтителен, когда целью является вычисление не коэффициентов ( x p ) , а только одного значения p ( a ) в точке x = a, отсутствующей в исходном наборе данных. Форма Лагранжа вычисляет значение p ( a ) со сложностью O( n ² ).

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

Интерполяции как линейные комбинации значений

Учитывая набор точек данных (позиция, значение) $(x_{0},y_{0}),\ldots ,(x_{j},y_{j}),\ldots ,(x_{n},y_{n})$ где нет двух позиций $x_{j}$ одинаковы, интерполяционный полином $y(x)$ можно рассматривать как линейную комбинацию значений $y_{j}$, используя коэффициенты, которые являются полиномами от $x$ в зависимости от $x_{j}$. Например, интерполяционный полином в форме Лагранжа представляет собой линейную комбинацию $$y(x):=\sum _{j=0}^{k}y_{j}c_{j}(x)$$с каждым коэффициентом $c_{j}(x)$ задается соответствующим базисным полиномом Лагранжа на заданных позициях $x_{j}$: $$c_{j}(x)=L_{j}(x_{0},\ldots ,x_{n};x)=\prod _{0\leq i\leq n \atop i\neq j}{\frac {x-x_{i}}{x_{j}-x_{i}}}={\frac {(x-x_{0})}{(x_{j}-x_{0})}}\cdots {\frac {(x-x_{j-1})}{(x_{j}-x_{j-1})}}{\frac {(x-x_{j+1})}{(x_{j}-x_{j+1})}}\cdots {\frac {(x-x_{n})}{(x_{j}-x_{n})}}.$$

Поскольку коэффициенты зависят только от позиций $x_{j}$, а не значения $y_{j}$, мы можем использовать те же коэффициенты , чтобы найти интерполяционный полином для второго набора точек данных $(x_{0},v_{0}),\ldots ,(x_{n},v_{n})$ на тех же позициях: $$v(x):=\sum _{j=0}^{k}v_{j}c_{j}(x).$$

Кроме того, коэффициенты $c_{j}(x)$ зависят только от относительных пространств $x_{i}-x_{j}$ между позициями. Таким образом, учитывая третий набор данных, точки которого задаются новой переменной $t=ax+b$ ( преобразование аффинное $x$, инвертированный $x={\tfrac {t-b}{a}}$): $$(t_{0},w_{0}),\ldots ,(t_{j},w_{j})\ldots ,(t_{n},w_{n})\qquad {\text{with}}\qquad t_{j}=ax_{j}+b,$$

мы можем использовать преобразованную версию предыдущих полиномов коэффициентов:

${\tilde {c}}_{j}(t):=c_{j}({\tfrac {t-b}{a}})=c_{j}(x),$

и запишите интерполяционный полином как:

${\textstyle w(t):=\sum _{j=0}^{k}w_{j}{\tilde {c}}_{j}(t).}$

Точки данных $(x_{j},y_{j})$ часто имеют равноотстоящие друг от друга позиции , которые можно нормализовать с помощью аффинного преобразования к $x_{j}=j$. Например, рассмотрим точки данных

$(0,y_{0}),(1,y_{1}),(2,y_{2})$.

Интерполяционный полином в форме Лагранжа представляет собой линейную комбинацию

$${\begin{aligned}y(x):=\sum _{j=0}^{2}y_{j}c_{j}(x)&=y_{0}{\frac {(x-1)(x-2)}{(0-1)(0-2)}}+y_{1}{\frac {(x-0)(x-2)}{(1-0)(1-2)}}+y_{2}{\frac {(x-0)(x-1)}{(2-0)(2-1)}}\\&={\tfrac {1}{2}}y_{0}(x-1)(x-2)-y_{1}(x-0)(x-2)+{\tfrac {1}{2}}y_{2}(x-0)(x-1).\end{aligned}}$$

Например, $y(3)=y_{3}=y_{0}-3y_{1}+3y_{2}$ и $y(1.5)=y_{1.5}={\tfrac {1}{8}}(-y_{0}+6y_{1}+3y_{2})$.

Случай равноотстоящих друг от друга точек также можно рассматривать методом конечных разностей . Первое отличие последовательности значений $v=\{v_{j}\}_{j=0}^{\infty }$ это последовательность $\Delta v=u=\{u_{j}\}_{j=0}^{\infty }$ определяется $u_{j}=v_{j+1}-v_{j}$. Повторение этой операции дает n ^й разностная операция $\Delta ^{n}v=u$, определенный явно: $$u_{j}=\sum _{k=0}^{n}(-1)^{n-k}{n \choose k}v_{j+k},$$где коэффициенты образуют знаковую версию треугольника Паскаля, треугольника коэффициентов биномиального преобразования :

								1									Строка n = 0
							1		−1								Строка n = 1 или d = 0
						1		−2		1							Строка n = 2 или d = 1
					1		−3		3		−1						Ряд n = 3 или d = 2
				1		−4		6		−4		1					Ряд n = 4 или d = 3
			1		−5		10		−10		5		−1				Ряд n = 5 или d = 4
		1		−6		15		−20		15		−6		1			Ряд n = 6 или d = 5
	1		−7		21		−35		35		−21		7		−1		Ряд n = 7 или d = 6

Полином $y(x)$ степени d определяет последовательность значений в положительных целых точках, $y_{j}=y(j)$и $(d+1)^{\text{th}}$ разность этой последовательности тождественно равна нулю:

$\Delta ^{d+1}y=0$.

Таким образом, заданные значения $y_{0},\ldots ,y_{n}$ в равноотстоящих точках, где $n=d+1$, у нас есть: $$(-1)^{n}y_{0}+(-1)^{n-1}{\binom {n}{1}}y_{1}+\cdots -{\binom {n}{n-1}}y_{n-1}+y_{n}=0.$$Например, 4 равноотстоящие друг от друга точки данных. $y_{0},y_{1},y_{2},y_{3}$ квадратичного $y(x)$ подчиняться $0=-y_{0}+3y_{1}-3y_{2}+y_{3}$и решение для $y_{3}$ дает то же интерполяционное уравнение, полученное выше с помощью метода Лагранжа.

Ошибка интерполяции: формула остатка Лагранжа

При интерполяции заданной функции f полиномом $p_{n}$ степени n в узлах x ₀ ,..., x _n получим ошибку $$f(x)-p_{n}(x)=f[x_{0},\ldots ,x_{n},x]\prod _{i=0}^{n}(x-x_{i})$$

где ${\textstyle f[x_{0},\ldots ,x_{n},x]}$ это ( n +1) ^ул. разделенная разность точек данных

$(x_{0},f(x_{0})),\ldots ,(x_{n},f(x_{n})),(x,f(x))$.

Более того, существует форма остатка Лагранжа ошибки для функции f , которая n + 1 раз непрерывно дифференцируема на замкнутом интервале. $I$и полином $p_{n}(x)$ степени не выше n , которая интерполирует f в n + 1 различных точках $x_{0},\ldots ,x_{n}\in I$. Для каждого $x\in I$ существует $\xi \in I$ такой, что

$$f(x)-p_{n}(x)={\frac {f^{(n+1)}(\xi )}{(n+1)!}}\prod _{i=0}^{n}(x-x_{i}).$$

Эта граница ошибки предполагает выбор точек интерполяции x _i, чтобы минимизировать произведение ${\textstyle \left|\prod (x-x_{i})\right|}$, что достигается узлами Чебышева .

Доказательство остатка Лагранжа

Установите термин ошибки как ${\textstyle R_{n}(x)=f(x)-p_{n}(x)}$и определим вспомогательную функцию: $$Y(t)=R_{n}(t)-{\frac {R_{n}(x)}{W(x)}}W(t)\qquad {\text{where}}\qquad W(t)=\prod _{i=0}^{n}(t-x_{i}).$$Таким образом: $$Y^{(n+1)}(t)=R_{n}^{(n+1)}(t)-{\frac {R_{n}(x)}{W(x)}}\ (n+1)!$$

Но поскольку $p_{n}(x)$ является многочленом степени не выше n , имеем ${\textstyle R_{n}^{(n+1)}(t)=f^{(n+1)}(t)}$, и: $$Y^{(n+1)}(t)=f^{(n+1)}(t)-{\frac {R_{n}(x)}{W(x)}}\ (n+1)!$$

Теперь, поскольку x _i являются корнями $R_{n}(t)$ и $W(t)$, у нас есть $Y(x)=Y(x_{j})=0$, что означает, что Y имеет как минимум n + 2 корня. По Ролля теореме $Y^{\prime }(t)$ имеет не менее n + 1 корней, и итеративно $Y^{(n+1)}(t)$ имеет хотя бы один корень в интервале I. ξ Таким образом: $$Y^{(n+1)}(\xi )=f^{(n+1)}(\xi )-{\frac {R_{n}(x)}{W(x)}}\ (n+1)!=0$$

и: $$R_{n}(x)=f(x)-p_{n}(x)={\frac {f^{(n+1)}(\xi )}{(n+1)!}}\prod _{i=0}^{n}(x-x_{i}).$$

Это соответствует рассуждениям, лежащим в основе остаточного члена Лагранжа в теореме Тейлора ; Фактически, остаток Тейлора представляет собой частный случай ошибки интерполяции, когда все узлы интерполяции x _i идентичны. Обратите внимание, что ошибка будет равна нулю, если $x=x_{i}$ для любого я . Таким образом, максимальная ошибка возникнет в какой-то момент интервала между двумя последовательными узлами.

Равноотстоящие интервалы

В случае равноотстоящих друг от друга узлов интерполяции, где $x_{i}=a+ih$, для $i=0,1,\ldots ,n,$ и где $h=(b-a)/n,$ член продукта в формуле ошибки интерполяции может быть связан как $$\left|\prod _{i=0}^{n}(x-x_{i})\right|=\prod _{i=0}^{n}\left|x-x_{i}\right|\leq {\frac {n!}{4}}h^{n+1}.$$

Таким образом, граница ошибки может быть задана как $$\left|R_{n}(x)\right|\leq {\frac {h^{n+1}}{4(n+1)}}\max _{\xi \in [a,b]}\left|f^{(n+1)}(\xi )\right|$$

Однако это предполагает, что $f^{(n+1)}(\xi )$ преобладает $h^{n+1}$, то есть $f^{(n+1)}(\xi )h^{n+1}\ll 1$. В некоторых случаях это не так, и ошибка фактически увеличивается при n → ∞ (см. феномен Рунге ). Этот вопрос рассматривается в разделе « Свойства сходимости» .

Константы Лебега

Основная статья: константа Лебега

Фиксируем узлы интерполяции x ₀ , ..., x _n и интервал [ a , b ], содержащий все узлы интерполяции. Процесс интерполяции отображает функцию f в полином p . Это определяет отображение X из пространства C ([ a , b ]) всех непрерывных функций на [ a , b ] в себя. Отображение X линейно и является проекцией на подпространство $P(n)$ многочленов степени n или меньше.

Лебега L определяется как операторная норма X Константа . Имеет место (частный случай леммы Лебега ): $$\left\|f-X(f)\right\|\leq (L+1)\left\|f-p^{*}\right\|.$$

Другими словами, интерполяционный полином не более чем в раз ( L + 1) хуже наилучшего возможного приближения. Это предполагает, что мы ищем набор узлов интерполяции, который делает L маленьким. имеем В частности, для узлов Чебышева : $$L\leq {\frac {2}{\pi }}\log(n+1)+1.$$

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

Свойства сходимости

Естественно задаться вопросом, для каких классов функций и для каких узлов интерполяции последовательность интерполяционных многочленов сходится к интерполируемой функции при n → ∞ ? Сходимость можно понимать по-разному, например, поточечно, равномерно или в некоторой интегральной норме.

Ситуация довольно плохая для эквидистантных узлов, поскольку равномерная сходимость не гарантируется даже для бесконечно дифференцируемых функций. Одним из классических примеров Карла Рунге является функция f ( x ) = 1 / (1 + x ² ) на интервале [−5, 5] . Ошибка интерполяции || ж - п _п || _∞ неограниченно растет при n → ∞ . Другой пример — функция f ( x ) = | х | на интервале [−1, 1] , для которого интерполяционные полиномы даже не сходятся поточечно, за исключением трех точек x = ±1, 0.

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

Теорема . Для любой функции f ( x ), непрерывной на интервале [ a , b ], существует таблица узлов, для которой последовательность интерполирующих полиномов $p_{n}(x)$ сходится к f ( x ) равномерно на [ a , b ].

Доказательство

Ясно, что последовательность полиномов наилучшего приближения $p_{n}^{*}(x)$ сходится к f ( x ) равномерно (в силу аппроксимационной теоремы Вейерштрасса ). Теперь нам осталось только показать, что каждый $p_{n}^{*}(x)$ могут быть получены путем интерполяции на определенных узлах. Но это справедливо благодаря особому свойству полиномов наилучшего приближения, известному из теоремы об эквиколебаниях . В частности, мы знаем, что такие многочлены должны пересекать f ( x ) не менее n + 1 раз. Выбирая точки пересечения в качестве узлов интерполяции, получаем интерполяционный полином, совпадающий с полиномом наилучшего приближения.

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

Теорема . Для любой таблицы узлов существует непрерывная функция f ( x ) на интервале [ a , b ], для которой последовательность интерполяционных полиномов расходится на [ a , b ].

Доказательство по существу использует оценку нижней оценки константы Лебега, которую мы определили выше как операторную норму X _n (где X _n — оператор проектирования на Π _n ). Теперь ищем таблицу узлов, для которых

$$\lim _{n\to \infty }X_{n}f=f,{\text{ for every }}f\in C([a,b]).$$

Согласно теореме Банаха–Штайнхауза это возможно только тогда, когда нормы X _n равномерно ограничены, что не может быть правдой, поскольку мы знаем, что

$$\|X_{n}\|\geq {\tfrac {2}{\pi }}\log(n+1)+C.$$

Например, если в качестве узлов интерполяции выбраны равноудаленные точки, функция из явления Рунге демонстрирует расходимость такой интерполяции. Обратите внимание, что эта функция не только непрерывна, но даже бесконечно дифференцируема на [−1, 1] . Однако для лучших узлов Чебышева такой пример найти гораздо сложнее из-за следующего результата:

Теорема . Для каждой абсолютно непрерывной функции на [−1, 1] последовательность интерполяционных многочленов, построенных на узлах Чебышёва, сходится к f ( x ) равномерно.

Связанные понятия

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

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

Задачи интерполяции Эрмита — это задачи, в которых заданы не только значения многочлена p в узлах, но и все производные до заданного порядка. Это оказывается эквивалентным системе одновременных сравнений полиномов и может быть решено с помощью китайской теоремы об остатках для многочленов. Интерполяция Биркгофа — это дальнейшее обобщение, в котором назначаются только производные некоторых порядков, а не обязательно всех порядков от 0 до k .

Коллокационные методы решения дифференциальных и интегральных уравнений основаны на полиномиальной интерполяции.

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

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

См. также

• серия Ньютона
• Полиномиальная регрессия

Ссылки

• Бернштейн, Сергей Н. (1912). «О порядке наилучшего приближения непрерывных функций многочленами заданной степени». Память акад. Рой. Бельгия. (на французском языке). 4 :1–104.
• Фабер, Георг (1914). «Об интерполяции непрерывных функций». Год немецкой математики. (на немецком языке). 23 : 192-210.
• Уотсон, Г. Алистер (1980). Теория приближений и численные методы . Джон Уайли. ISBN  0-471-27706-1 .

Дальнейшее чтение

• Аткинсон, Кенделл А. (1988). «Глава 3.». Введение в численный анализ (2-е изд.). Джон Уайли и сыновья. ISBN  0-471-50023-2 .
• Брутман, Л. (1997). «Функции Лебега для полиномиальной интерполяции — обзор». Энн. Число. Математика . 4 : 111–127.
• Пауэлл, MJD (1981). «Глава 4». Теория и методы приближения . Издательство Кембриджского университета. ISBN  0-521-29514-9 .
• Шацман, Мишель (2002). «Глава 4». Численный анализ: математическое введение . Оксфорд: Кларендон Пресс. ISBN  0-19-850279-6 .
• Сюли, Эндре ; Майерс, Дэвид (2003). «Глава 6». Введение в численный анализ . Издательство Кембриджского университета. ISBN  0-521-00794-1 .

Внешние ссылки

• «Процесс интерполяции» , Математическая энциклопедия , EMS Press , 2001 [1994]
• ALGLIB имеет реализации на C++/C#.
• GSL имеет код полиномиальной интерполяции на языке C.
• Демонстрация полиномиальной интерполяции .