Погрешности алгебраической интерполяции

Погрешность метода

Хотим понять, насколько точно мы восстанавливаем исходную функцию \(f(x)\) по дискретному набору \(\left\{ \left(x_i, y_i = f(x_i)\right) \right\}_{i=0}^{n}\) с помощью интерполяционного многочлена \(L_n(x)\). Для этого введём функцию погрешности как

\[ R_n(x) = f(x) - L_n(x) \]

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

Пусть даны \(f(x) \in \text{C}^{\,n+1}[a, b]\) и \(x_0, x_1, \ldots, x_n\) - строго возрастающий набор точек отрезка \([a, b]\). Пусть \(L_n(x)\) - интерполяционный многочлен построенный на наборе \(\left\{ \left(x_i, y_i = f(x_i)\right) \right\}_{i=0}^{n}\)

Тогда погрешность интерполяции \(R_n(x)=f(x)-L_n(x)\) представляется формулой

\[ R_n(x)=\frac{1}{(n+1)!} f^{(n+1)}(\xi) \cdot \omega (x) \]

\[ \omega(x) = \left(x-x_0\right)\left(x-x_1\right) \ldots\left(x-x_n\right) \]

где \(\xi=\xi(x)\) - некоторая точка отрезка \([a, b]\)

Тогда несложно оценить погрешность (метода) интерполяции как

\[ \varepsilon_{\text{method}} \leq \frac{\max_{x \in [a,b]} |f^{(n+1)}(x)|}{(n+1)!} \cdot \max_{x \in [a,b]} |\omega(x)| = \cfrac{M_{n+1}}{(n+1)!} \cdot ||\omega(x)||_\text{C} \]

Заметим, что \(M_{n+1}\) зависит только от функции, а \(||\omega(x)||_\text{C}\) только от расположения узлов. Таким образом возникает идея о выборе оптимального расположения узлов \(\{x_k\}_{k=0}^{n}\) для минимизации нормы \(||\omega(x)||_\text{C}\), а следовательно и \(\varepsilon_{\text{method}}\).

Несложно показать, что этот оптимальный набор \(\{x_k\}_{k=0}^{n}\), минимизирующий оценку на погрешность метода \(\varepsilon_\text{method}\), является корнями многочленов Чебышева:

\[ x_k=\frac{a+b}{2}+\frac{b-a}{2} \cos \left(\frac{2 (n - k) + 1}{2 (n + 1)} \cdot \pi\right), \quad k \in \overline{0,...,n} \]

Более того, верна следующая теорема.

Теорема о сходимости интерполяционного процесса

Если функция \(f(x)\) имеет ограниченную производную на отрезке, то интерполяционный процесс, в котором за узлы принимаются корни многочленов Чебышёва, сходится равномерно к \(f(x)\).

\[ L_n(x)\rightrightarrows f(x) \]

Заметим, что из равномерной сходимости функций не следует равномерная сходимость производных. Вообще, оценивать производные неизвестной \(f(x)\) по известной полиномиальной интерполяции \(L_n(x)\) может быть плохой идеей из-за осцилляций.

Погрешность данных (обусловленность интерполяции)

Естественно, что для произвольной функции мы никогда не можем абсолютно точно знать её значения \(y_i\) - есть какая-то ошибка \(\Delta y\). Тогда интерполяционный многочлен в форме Лагранжа

\[ L_n(x)=\sum_{i=0}^n y_i \cdot l_i(x) \]

А тогда, навешивая погрешность слева и справа, получаем

\[ \Delta L_n(x) = \Delta y \cdot \sum_{i=0}^n |l_i(x)| \]

Функция \(L(x) = \sum_{i=0}^n |l_i(x)|\) называется функцией Лебега и зависит только от расположения узлов. Таким образом, для погрешности интерполянта, вызванного ошибкой в данных (шумом) верна оценка

\[ \varepsilon_{\text{noise}} \leq \Delta y \cdot \max_{x \in [a,b]} \sum_{i=0}^n |l_i(x)| = \Delta y \cdot L \]

где введена т.н. константа Лебега, ограничивающая ошибку интерполяции сверху и зависящюю только от расположения узлов \(\{x_k\}_{k=0}^{n}\):

\[ L = \max_{x \in [a,b]} \sum_{i=0}^n |l_i(x)| \]

Для равномерной сетки константа Лебега зависит только от числа узлов

• Для линейной интерполяции \((n=2)\) константа Лебега \(L=1\).

• Для квадратичной интерполяции ( \(n=3\) ) константа Лебега \(L=1.25\)

• При \(n=10\) константа Лебега \(L \approx 19\).

• При \(n=20\) константа Лебега \(L \approx 6900\).

• При \(n=30\) константа Лебега \(L \approx 4 \cdot 10^6\).

• При \(n \gg 1\) константа Лебега растет как \(L \sim \frac{2^n}{e(n-1) \ln (n-1)}\)

Напротив, для сетки из нулей многочленов Чебышева (которая минимизирует ошибку метода), константа Лебега оказывается довольно малой:

• Для линейной интерполяции \((n=1)\) константа Лебега \(L=\sqrt{2}\).

• Для квадратичной интерполяции \((n=2)\) константа Лебега \(L=5 / 3 \approx 1.6667\).

• При \(n=10\) константа Лебега \(L \approx 2.4288\).

• При \(n=20\) константа Лебега \(L \approx 2.8698\).

• При \(n=30\) константа Лебега \(L \approx 3.1278\).

• В общем случае \(\frac{2}{\pi} \ln n+0.96<L<\frac{2}{\pi} \ln n+1\)

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

Мы показали, что интерполировать, используя большое количество узлов некорректно. Начиная с \(n=50\) не хватает даже машинной точности, да и производные строго говоря мы не ловим (а хотелось бы приближать и их).

Интерполяция и экстраполяция высокого порядка точности*

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

\[y_{n+1}(x)=F\left(x ; y_{n}(x), y_{n}^{\prime}(x), y_{n}^{\prime \prime}(x), \ldots\right) \ \ \ \ (1) \]

где \(F\left(x ; y(x), y^{\prime}(x), y^{\prime \prime}(x), \ldots\right)\) — функционал от \(y(x), y^{\prime}(x), y^{\prime \prime}(x), \ldots\), который зависит также от параметра x,

\(n=0,1,2\) - номер итерации, а начальные функции \(y_{0}(x), y_{0}^{\prime}(x), y_{0}^{\prime \prime}(x), \ldots\) пусть заданы.

Уравнения типа (1) удобно решать на сетке \(x=x_{j}, j \in Z-\) целое, вычисляя правую часть (1) последовательно для каждого \(x_{j}\) при каждом следующем значении \(n . \mathrm{B}\) связи с этим возникает необходимость решения задачи о численном нахождении функции \(f(x)\) и ее производных \(f^{(k)}(x)\) за пределами точек сетки \(\left(x \neq x_{j} \forall j\right)\) по известным значениям функции на сетке \(f\left(x_{j}\right), j \in Z .\)

Решение:

Пусть функция \(f(x)\) аналитична и задана на сетке \(x_{j}, j \in Z\) - целое, как \(f\left(x_{j}\right)=y_{j}\) с ошибкой \(\Delta y_{j}\). Сетка может быть как эквидистантной \(x_{j}=j \Delta x, \Delta x>0-\) шаг, так и неэвидистантной \(x_{j}=x\left(t_{j}\right)\), где значения \(t_{j}=j \Delta t-\) эквидистантны, \(\Delta t>0 .\) При этом мы полагаем, что функция \(x(t)\) аналитична и представляет собой несложное аналитическое выражение, которое легко обращается и дифференцируется в аналитическом виде.

Целью настоящей разработки является численное нахождение следующих величин:

a) Значения функции \(f(x)\) в промежуточных точках \(x_{j}<x<\) \(x_{j+1}\) между точками сетки.

б) Значения производных \(f^{\prime}(x), f^{\prime \prime}(x), \ldots\) в точках сетки \(x=x_{j}\) и между точками сетки \(x_{j}<x<x_{j+1}\).

в) Аналитическое продолжение функции \(f(x)\) и ее производных \(f^{\prime}(x), f^{\prime \prime}(x), \ldots\) в плоскость комплексных значений \(x\).

г) Аналитическое продолжение функции \(f(x)\) за пределы сетки (т. е. за пределы области ее прямого численного вычисления).

Случай (г) может потребоваться, например, в случае, когда в (1) функция \(y(x)\) численно вычисляется на отрезке \(x \in[a ; b]\), a нужно знать ее значения при некоторых \(x<a\) или \(x>b\).

Задачи (а)-(г) формально решаются путем разложения \(f(x)\) в ряд Тейлора вблизи некоторого узла сетки \(x_{j}\)

\[ f^{(k)}(x)=\sum_{n=0}^{\infty} \frac{f^{(n+k)}\left(x_{j}\right)}{n !}\left(x-x_{j}\right)^{n} \ \ \ \ (2) \]

Здесь \(x_{j}\) - ближайшая к \(x\) точка сетки (т. е. \(\left.0 \leq\left|x-x_{j}\right| \leq \Delta x / 2\right)\), и мы здесь полагаем сетку эквидистантной (неэквидистантные сетки будут рассмотрены ниже).

Из (2) видно, что для фактического решения задач (а)-(г) необходимо вычислить производные \(f^{(n)}\left(x_{j}\right)\) в точках сетки \(x_{j}\), \(j \in Z\), через значения функции на сетке \(y_{j}, j \in Z .\) При этом необходимо также рассчитать ошибку вычисления этих производных, которая связана как с ошибкой выражения \(f^{(n)}\left(x_{j}\right), j \in Z\), через \(y_{j}, j \in Z\), так и с исходной ошибкой \(\Delta y_{j}\) самих \(y_{j} .\) В результате, полная ошибка вычисления правой части (2) будет ограничивать значения \(x\), до которых экстраполяция функции \(f(x)\) и ее производных за пределы сетки численно возможна.

Для вычисления входящих в (2) производных на сетке \(f^{(n)}\left(x_{j}\right)\), \(j \in Z\), выразим производные \(f^{(n)}(x)\) в точке \(x=x_{j}\) через значения функции \(f(x)\) в точках \(x=x_{j}\).

В главном (нулевом) порядке по шагу \(\Delta x\) для производных \(f^{(n)}(x)\) имеет место следующая формула

\[ f^{(n)}(x)=\frac{1}{2^{n} \Delta x^{n}} \sum_{k=0}^{n}(-1)^{k} C_{n}^{k} f(x+(n-2 k) \Delta x)+O\left(\Delta x^{2}\right) \ \ \ \ (3) \]

\[\begin{split} C_{n}^{k}= \begin{cases}\frac{n !}{k !(n-k) !}, & 0 \leq k \leq n, \\ 0, & k<0 \text { или } k>n\end{cases} \ \ \ \ \ (4) \end{split}\]

биноминальные коэффициенты, симметричность которых \(C_{n}^{k}=C_{n}^{n-k}-\) обеспечивает обращение в нуль всех нечетных порядков по \(\Delta x\), т. е. \(O\left(\Delta x^{2 m+1}\right)=0\) в (3).

Формула (3) справедлива до \(O\left(\Delta x^{2}\right)\). Нам же нужно знать производные \(f^{(n)}(x)\) в произвольном порядке по шагу \(\Delta x\). Чтобы их вычислить, будем последовательно итерировать (3) с помощью разложения функции \(f(x)\) в ряд Тейлора

\[ f(x+a)=\sum_{n=0}^{\infty} \frac{f^{(n)}(x)}{n !} a^{n} .\ \ \ \ \ (5) \]

Именно, подставляя (5) в правую часть (3) при \(a=(n-2 k) \Delta x\) и учитывая, что разностная схема (3) не содержит нечетных порядков по \(\Delta x\), распишем входящую в правую часть (3) величину \(O\left(\Delta x^{2}\right)\) до следующего порядка

\[\begin{split} \begin{gathered} O\left(\Delta x^{2}\right)=-\frac{1}{2^{n} \Delta x^{n}} \sum_{k=0}^{n}(-1)^{k} C_{n}^{k} \frac{f^{(n+2)}(x)}{(n+2) !}(n-2 k)^{n+2} \Delta x^{n+2}+ \\ \quad+O\left(\Delta x^{4}\right)= \\ =-D_{n}^{(1)} f^{(n+2)}(x)(2 \Delta x)^{2}+O\left(\Delta x^{4}\right) \end{gathered} \ \ \ \ \ \ (6) \end{split}\]

где

\[ D_{n}^{(1)}=\sum_{k=0}^{n}(-1)^{k} C_{n}^{k} \frac{(n-2 k)^{n+2}}{2^{n+2}(n+2) !} . \]

Подставляя в (6) производную \(f^{(n+2)}(x)\) из (3), для входящей в правую часть (3) величины \(O\left(\Delta x^{2}\right)\) получаем

\(O\left(\Delta x^{2}\right)=\frac{1}{2^{n} \Delta x^{n}} \sum_{k=0}^{n+2}(-1)^{k-1} D_{n}^{(1)} C_{n+2}^{k} f(x+(n+2-2 k) \Delta x)+O\left(\Delta x^{4}\right) \ \ \ \ (7) \),

где учтено, что \(-(-1)^{k}=(-1)^{k-1}\).

Перепишем (3), сделав замену \(k \rightarrow k-1 .\)

Имеем,

\(f^{(n)}(x)=\frac{1}{2^{n} \Delta x^{n}} \sum_{k=1}^{n+1}(-1)^{k-1} C_{n}^{k-1} f(x+(n+2-2 k) \Delta x)+O\left(\Delta x^{2}\right)=\) \(=\frac{1}{2^{n} \Delta x^{n}} \sum_{k=0}^{n+2}(-1)^{k-1} C_{n}^{k-1} f(x+(n+2-2 k) \Delta x)+O\left(\Delta x^{2}\right) \ \ \ \ (8) \)

где \(O\left(\Delta x^{2}\right)-\) то же самое, что и в (3) и в левой части (7) и где учтено, что \(C_{n}^{k}=0\) при \(k<0\) или \(k>n\) (см. (4)).

Подставляя \(O\left(\Delta x^{2}\right)\) в виде (7) в формулу (8), получаем производную \(f^{(n)}(x)\) с точностью до второго порядка по \(\Delta x\) включительно

\[\begin{split} \begin{gathered} f^{(n)}(x)=\frac{1}{2^{n} \Delta x^{n}} \sum_{k=0}^{n+2}(-1)^{k-1}\left(C_{n}^{k-1}+D_{n}^{(1)} C_{n+2}^{k}\right) f(x+(n+2-2 k) \Delta x)+ \\ +O\left(\Delta x^{4}\right) . \end{gathered} \end{split}\]

Итерируя эту процедуру далее, получаем производные \(f^{(n)}(x)\) с точностью до \(2 m\)-го порядка включительно

\(f^{(n)}(x)=\frac{1}{2^{n} \Delta x^{n}} \sum_{k=0}^{n+2 m} A_{k n}^{m} f(x+(n+2 m-2 k) \Delta x)+O\left(\Delta x^{2 m+2}\right) \ \ \ \ \ (9) \)

где

\[ A_{k n}^{m}=\sum_{l=0}^{m}(-1)^{k-m} D_{n}^{(l)} C_{n+2 l}^{k-m+l} \ \ \ \ (10) \]

коэффициенты разложения, а

\[ D_{n}^{(j+1)}=\sum_{k=0}^{n+2 j} \sum_{l=0}^{j}(-1)^{k} D_{n}^{(l)} C_{n+2 l}^{k-j+l} \frac{(n+2 j-2 k)^{n+2 j+2}}{2^{n+2 j+2}(n+2 j+2) !} \ \ \ \ (11) \]

реккурентная формула для входящих в (10) величин \(D_{n}^{(l)}\), в которой положено \(D_{n}^{(0)} \equiv 1 .\)

Формула (9) является основной вспомогательной формулой для расчета сетки для производных \(f^{(n)}\left(x_{j}\right)\) в \((2)\).

Ради примера, с помощью (9) выпишем \(f^{\prime}(x)\) с точностью до 0 -го, 2-го и 4-го порядка по \(\Delta x\), включительно,

\[\begin{split} \begin{gathered} f^{\prime}(x)=\frac{f(x+\Delta x)-f(x-\Delta x)}{2 \Delta x}+O\left(\Delta x^{2}\right)= \\ =\frac{27(f(x+\Delta x)-f(x-\Delta x))-(f(x+3 \Delta x)-f(x-3 \Delta x))}{48 \Delta x}+ \\ +O\left(\Delta x^{4}\right)= \\ =\frac{1}{80 \cdot 48 \Delta x}[2250(f(x+\Delta x)-f(x-\Delta x))-125(f(x+3 \Delta x)- \\ -f(x-3 \Delta x))+9(f(x+5 \Delta x)-f(x-5 \Delta x))]+O\left(\Delta x^{6}\right) . \end{gathered} \end{split}\]

Случай эквидистантной сетки.

Рассмотрим в (9) случай эквидистантной сетки: \(x_{j}=j \Delta x\), \(j \in Z .\) Подставляя \(x=x_{j}\) в (9) и учитывая, что \(f\left(x_{j}\right)=y_{j}\) (см. разд. 1) и \(x_{j}+(n+2 m-2 k) \Delta x=x_{j+n+2 m-2 k}\), получаем входящие в (2) искомые выражения для производных на сетке

\[ f^{(n)}\left(x_{j}\right)=\frac{1}{2^{n} \Delta x^{n}} \sum_{k=0}^{n+2 m} A_{k n}^{m} y_{j+n+2 m-2 k}+O\left(\Delta x^{2 m+2}\right) \ \ \ \ \ (12) \]

где коэффициенты \(A_{k n}^{m}\) даются формулами (10) и (11). Через выражения (12) для производных на сетке \(f^{(n)}\left(x_{j}\right)\) peшаются по формуле (2) поставленные нами задачи (а)-(г).