Интерполяция функций

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

Пусть задана совокупность узлов интерполяции или сетка на некотором отрезке [a, b]. В простейшем случае сетка — равномерная, т.е. расстояние между соседними узлами одинаково. В дальнейшем также рассмотрим неравномерные сетки.

Совокупность узлов $$\left\{{t_n}\right\}_{n = 0}^N,$$ $$t_{n} = a + n\tau,$$ $$\tau = (b - a)/N,$$ $$t \in \left[{a, b}\right].$$

Сеточная проекция функции f(t) на [a, b], т.е. таблица $$f_n = = \left\{{f(t_n)}\right\}_{n = 0}^N $$ ; эту таблицу задает оператор ограничения на сетку или рестрикции (от английского restriction ) $$\mathbf{R}.$$

Задача состоит в том, чтобы по таблице {fn} восстановить непрерывную функцию. Обозначим ее через F(t). Разумеется, она отличается от исходной функции f(t), причем такое восстановление неоднозначно и осуществляется оператором интерполяции $$\mathbf{I}.$$ Сама функция F(t) называется интерполирующей или интерполянтом. Необходимо оценить потерю информации при действии этого оператора, т. е. величину |f(t) - F(t)|, зависящую от типа оператора интерполяции и свойств f(t), в частности, ее гладкости. Таким образом, имеем схему:

$$f(t) \mathop \to\limits_{\mathbf{R}} \left\{{f_n}\right\}_{n = 0}^{N} \mathop \to\limits_{\mathbf{I}} F(t).$$

6.2. Кусочно - линейная интерполяция

Простейший способ интерполяции — кусочно - линейная, требующая минимальных требований на гладкость функции f(t). При таком способе интерполяции соседние точки ( tn, fn ) и ( tn + 1, fn + 1 ) соединяют отрезками прямых

$$$ F(t) = \frac{f_{n + 1} (t - t_n) + f_n (t_{n + 1} - t)}{t_{n + 1} - t_n}, t \in [t_n , t_{n + 1} ]. $$$

Теорема. Пусть f(t) — Липшиц непрерывная функция, т.е. $$| f(t_1 ) - f(t_2 ) | \le c \left|{t_1 - t_2 }\right|,$$ тогда $$$ \left|{f(t) - F(t)}\right| \le c\frac{\tau }{2} $.$$

Примечание. Если сетка неравномерная и $$\tau = \max\limits_n (t_{n + 1} - t_n),$$ то теорема верна и для этого случая.

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

Пусть $$t \in [t_n , t_{n + 1} ],$$ обозначим $$\tau = t_{n + 1} - t_{n}.$$ Тогда $$t = t_n + \alpha \cdot \tau$$ ; $$0 \le \alpha \le 1.$$ В силу линейности f(t) имеем равенство $$F(t) = \alpha f_{n + 1} + (1 - \alpha )f_n.$$

Оценим разность

$$\begin{gather*} \left|{F(t) - f(t)}\right| = \left|{\alpha f_{n + 1} + (1 - \alpha )f_n - \alpha f(t) - (1 - \alpha )f(t)}\right| \le \\ \le \alpha \left|{f_{n + 1} - f(t)\left. {}\right| + (1 - \alpha )}\right|f_n - f(t)\left. {}\right|. \end{gather*} $$

Поскольку $$f_{n + 1} = f(t_{n} + \tau ),$$ имеем

$$\begin{multline*} \left|f_{n + 1} - f(t)\right| = \left|f(t_n + \tau) - f(t_n + \alpha\tau)\right| \le \left|c(1 - \alpha )\tau\right| = \\ = c(1 - \alpha)\tau, \quad \mbox{т.к. }\quad 0 \le \alpha \le 1. \end{multline*}$$

Аналогично $$\left|{f_n - f(t)}\right| \le c\alpha \tau.$$ В таком случае $$\left|{f(t) - F(t)}\right| \le 2\alpha (1 - \alpha ) c\tau \le c\tau/2.$$

Замечание. Простой аппарат кусочно - линейной интерполяции позволяет ввести объекты, на которых базируется один из наиболее известных современных численных методов — метод конечных элементов. Сетке { tn } ставится в соответствие набор базисных функций $$\varphi_n(t),$$ каждая из которых сопоставляется своему узлу tn, причем $$\varphi_n (t_k ) = \delta_k^{n}, \varphi_n (t_{n - 1}) = \varphi_n (t_{n + 1}) = 0, \varphi_n (t_n) = 1,$$ а в остальных точках она вычисляется с помощью кусочно - линейной интерполяции.

Функция f(t) в этом случае представляется в виде

$$F(t) = \sum\limits_{n = 0}^{N}{f_n}\varphi_n (t).$$

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

6.3. Интерполяция обобщенными полиномами

Для того чтобы функция ( обобщенный полином )$$F(t) = \sum\limits_{n = 0}^{N}{u_n \varphi_n (t)}$$ была интерполирующей, необходимо выполнение условий: F(tk) = fk, $$k = 0 \div N,$$ где fk — значения функции в точках интерполяции. Для коэффициентов обобщенного полинома получаем систему уравнений:

$$\left\{ \begin{array}{l} u_0 \cdot \varphi_0 (t_0 ) + u_1 \cdot \varphi_1 (t_0 ) + \ldots + u_N \cdot \varphi_N (t_0 ) = f_0, \\ u_0 \cdot \varphi_0 (t_1 ) + u_1 \cdot \varphi_1 (t_1 ) + \ldots + u_N \cdot \varphi_N (t_1 ) = f_1, \\ \ldots \\ u_0 \cdot \varphi_0 (t_N) + u_1 \cdot \varphi_1 (t_N) + \ldots + u_N \cdot \varphi_N (t_N) = f_N , \\ \end{array} \right.$$

или в векторной форме

$${\mathbf{Au}} = {\mathbf{f}},$$

где

$${\mathbf{A}} = \left( \begin{array}{cccc} {\varphi_0 (t_0 )} {\varphi_1 (t_0 )} \ldots {\varphi_N (t_0 )} \\ {\varphi_0 (t_1 )} {\varphi_1 (t_1 )} \ldots {\varphi_N (t_1 )} \\ \ldots \ldots \ldots \ldots \\ {\varphi_0 (t_N)} {\varphi_1 (t_N)} \ldots {\varphi_N (t_N)} \\ \end{array} \right), {\mathbf{u}} = \left( \begin{array}{cccc} {u_0 } {u_1 } \ldots {u_N} \\ \end{array} \right)^{T}, \\ {\mathbf{f}} = \left( \begin{array}{cccc} {f_0 } {f_1 } \ldots {f_N} \\ \end{array} \right)^{T}.$$

Теорема (доказывается в курсе линейной алгебры.) Для того чтобы решение задачи интерполяции существовало и было единственным, необходимо и достаточно, чтобы система базисных функций $$\varphi_N (t_k)$$ была линейно независима.

Теорема (доказывается в курсе линейной алгебры.) Для того чтобы система функций $$\varphi_N (t_k)$$ была линейной независимой в точках t0, ..., tn, необходимо и достаточно, чтобы определитель матрицы Грама

$${\mathbf{C}} = {\mathbf{A}}^\ast {\mathbf{A}} = \left( \begin{array}{cccc} {(\varphi_0, \varphi_0)} {(\varphi_0, \varphi_1)} \ldots {(\varphi_0, \varphi_N)}\\ {(\varphi_1, \varphi_0)} {(\varphi_1, \varphi_1)} \ldots {(\varphi_1, \varphi_N)}\\ \ldots \ldots \ldots \ldots \\ {(\varphi_N , \varphi_0)} {(\varphi_N , \varphi_1)} \ldots {(\varphi_N , \varphi_N)} \end{array} \right),$$

был отличен от нуля. Здесь каждый элемент матрицы Грама имеет вид

$$\gamma_{jk} = (\varphi_k , \varphi_j) = \sum\limits_{i = 0}^{N}{\varphi_k (t_i) \cdot \varphi_j (t_i)}.$$

В случае, если система функций $$\left\{{\varphi_j}\right\}_0^{N}$$ ортогональна на множестве точек $$\left\{{t_j}\right\}_0^{N},$$ решение задачи интерполяции значительно упрощается (напомним, что система функций $$\left\{{\varphi_j}\right\}_0^{N}$$ является ортогональной на множестве точек $$\left\{{t_j}\right\}_0^{N},$$ если $$(\varphi_k , \varphi_j) = 0$$ при $$k \ne j$$ и $$\left({\varphi_k , \varphi_j}\right) \ne 0$$ при k = j для всех k = 0, 1, ..., N ; j = 0, 1, ..., n ).

Дело в том, что матрица Грама для ортогональной системы функций диагональна, и ее определитель отличен от нуля (всякая ортогональная система функций заведомо линейно независима). Линейная система уравнений представляется как $${\mathbf{A}}*{\mathbf{Au}} = {\mathbf{A}}*{\mathbf{f}},$$ или $$\mathbf{Cu} = \mathbf{b},$$ где $$\mathbf{C} = \mathbf{A}*\mathbf{A},$$ $$\mathbf{b} = \mathbf{A}*\mathbf{f}$$ - вектор, а ее решение в случае $$\mathbf{A}*\mathbf{A} = \mathbf{E}$$ есть $${\mathbf{u}} = {\mathbf{A}}*{\mathbf{f}}.$$

Примером ортогональной системы являются показательные функции $$e^{2\pi ikt_j}$$ на множестве точек tj = {j / N}, j = 0, 1, ..., N (на отрезке [0, 1] ).

6.4. Полиномиальная (алгебраическая) интерполяция

В этом случае (uk(t) = tk) СЛАУ для определения коэффициентов имеет вид

$$\left\{ \begin{array}{l} u_0 + u_1 t_0 + \ldots + u_N t_0^{N} = f_0, \\ u_0 + u_1 t_1 + \ldots + u_N t_1^{N} = f_1, \\ \ldots \\ u_0 + u_1 t_N + \ldots + u_N t_N^{N} = f_N , \\ \end{array} \right.$$

а ее определитель

$$\det \left( \begin{array}{cccc} 1 {t_0 } {t_0^2 } \ldots \\ 1 {t_1 } {t_1^2 } \ldots \\ \cdots \cdots \cdots \cdots \\ 1 {t_N} {t_N^2 } \ldots \\ \end{array} \right) = \mathop {\Pi\limits_{i \ne j}}\limits^{N}(t_i - t_j), 0 \le j < i \le N,$$

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

Ответ на вопрос о существовании и единственности решения СЛАУ оказывается — утвердительным - решение задачи алгебраической интерполяции всегда существует и единственно, но при больших N система оказывается плохо обусловленной. Однако решение этой задачи можно выписать в явном виде

$$L_N (t) = \sum\limits_{n = 0}^{N}{f_n} \cdot \varphi_n^{N} (t),$$

где $$$ \varphi_n^{N} (t) = {\mathop \Pi\limits^N_{\substack{i = 0 \\ i \ne n}}} \frac{t - t_j}{t_n - t_i} $$$ - базисные функции, являющиеся полиномами степени N, каждый из которых сопоставлен со своим узлом сетки так, что $$\varphi_n^{N} (t_k) = \delta_k^{n} .$$ Заметим, что правильнее было бы писать LN(t, {tn}, {fn}), т.е. интерполянт зависит от t, сетки и сеточной функции. Такой вид записи алгебраического интерполяционного полинома не единственен. Выписанный полином называется интерполяционным полиномом в форме Лагранжа. Он удобен для теоретического рассмотрения, но на практике часто оказывается более удобной другая форма представления — полином в форме Ньютона, о котором речь пойдет ниже.

6.5. Теорема об остаточном члене интерполяции

Введем понятие остаточного члена интерполяции для оценки погрешности

$$R_N (t) = f({t}) - L_N (t).$$

Теорема. Пусть функция f(t) имеет на отрезке [a, b] — N + 1 ограниченную производную. Тогда $$$ R_N (t) = \frac{1}{(N + 1)!} {\mathop \Pi\limits_{j = 0}^{N} (t - t_j)} \cdot f^{(N + 1)}(\xi ),$$ где $$\xi \in \left[{a, b}\right] $.$$

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

Рассмотрим функцию

$$$ {\psi}(x) = f(x) - L_N (x) - R_N (t)\frac{(x - t_0 )(x - t_1 ) \ldots (x - t_N)}{(t - t_0 )(t - t_1 ) \ldots (t - t_N)}, $$$

имеющую, по крайней мере, N + 1 производную. По условию, эту производную имеет f(x), а два остальных члена — полиномы.

Кроме того, $${\psi}(x)$$ на [a, b] имеет, по крайней мере, N + 2 нуля.

Их можно указать. Точки x = tn (n = 0, …, N) — нули, поскольку f(tn) = L(tn), а последнее слагаемое обращается в них в нуль. N + 2 нулем является точка x = t в силу определения остаточного члена. Далее, поскольку между каждыми двумя нулями непрерывно дифференцируемой функции имеется хотя бы один нуль ее производной, на [a, b] имеется хотя бы N + 1 нуль $${\psi}^{\prime}.$$ Применяя это рассуждение к $$\psi^{\prime\prime}, \psi^{\prime\prime\prime}, \ldots $$ можно показать, что существует точка $$\xi \in [a, b]$$ такая, что $${\psi}^{(N + 1)}(\xi ) = 0.$$

Вычислим N + 1 производную правой части выражения для f(x) с учетом того, что L(N + 1) = 0. Кроме того, в точке $$\xi$$

$$\begin{gather*} {\psi}^{(N + 1)} (\xi ) = f^{(N + 1)} (\xi ) - L^{(N + 1)} (\xi ) - \frac{d^{N + 1}}{dx^{N + 1}}\left[{R_N (t) \cdot \frac{(x - t_0 ) \ldots (x - t_N)}{(t - t_0 ) \ldots (t - t_N)}}\right]_\xi, \\ L^{(N + 1)} (\xi ) = 0; {\psi}^{(N + 1)} (\xi ) = 0; \\ \frac{d^{N + 1}}{dx^{N + 1}} \left. \left[\frac{(x - t_0 ) \ldots (x - t_N)}{(t - t_N) \ldots (t - t_N)}\right]\right|_{x = \xi} = \frac {(N + 1)!}{\mathop \Pi\limits_{j = 0}^N (t - t_j)} \end{gather*}$$ Тогда $$$ f^{(N + 1)} (\xi ) - R_N (t) \cdot \frac{{(N + 1)!}}{{ \mathop \Pi\limits_{j = 0}^{N} (t - t_j)}} = 0 $,$$ откуда получим выражение для RN(t): $$$ R_N (t) = \frac{{f^{(N + 1)} (\xi )}}{{(N + 1)!}}{ \mathop \Pi\limits_{j = 0}^N (t - t_j)} $$$

Рассмотрим некоторые важные следствия этой теоремы.

Следствие (точность интерполяции на равномерной сетке). Положим, что $$t_n = n\tau, \tau = (b - a)/N, t \in \left[{a, b}\right],$$ — сетка равномерная. В этом случае имеет место оценка

$$$ \left|{R_N (t)}\right| \le \frac{{\tau ^{N + 1}}}{{N + 1}}C, C = \max\limits_{t \in \left[{a, b}\right]} \left|{f^{(N + 1)}(t)} \right|. $$$

Доказательство. Пусть $$t = t_k + \alpha \tau , \alpha \in \left[{0, 1} \right], k = 0, 1, \ldots , N - 1.$$

Тогда $$t - t_n = k\tau + \alpha \tau - n\tau = (k + \alpha - n)\tau $$ ; откуда $${\mathop \Pi\limits_{n = 0}^{N} (t - t_n)} = \tau ^{N + 1} {\mathop \Pi\limits_{n = 0}^{N} (k + \alpha - n)}.$$ Можно показать, что $${\mathop \Pi\limits_{n = 0}^{N}} \left|{k + \alpha - n}\right| \le N$$! . Остаточный член оценивается следующим образом: $$$ R_N (t) = \frac{f^{(N + 1)} (\xi )}{(N + 1)!} {\mathop \Pi\limits_{n = 0}^{N} (t - t_n)} $,$$ поэтому с учетом приведенных оценок получим

$$$ \left|R_N(t)\right| \le \frac {\tau^{N + 1}}{N + 1} \max\limits_{\xi \in \left[{a, b}\right]} \left|{f^{(N + 1)} (\xi)}\right|. $$$

Рассмотрим, как ведет себя оценка в задаче экстраполяции при удалении точки t от интервала [t0, tN]. При $$t \in \left[{t_N , t_N + \tau }\right]$$ имеем $$\left|{R_N (t)}\right| \le \tau ^{N + 1} \cdot \max\limits_{\xi \in [t_0, t_N + \tau ]} \left|{f^{(N + 1)} (\xi )}\right|,$$ поскольку $${\mathop \Pi\limits_{n = 0}^{N + 1}} \left|{k + \alpha - n}\right| \le (N + 1)$$!. При$$t \in \left[{t_N + \tau , t_N + 2\tau }\right] \left|{R_N (t)}\right| \le (N + 2) \tau ^{N + 1} \cdot \max\limits_{\xi \in [t_0, t_N + 2\tau ]} \left|{f^{(N + 1)} (\xi )}\right|,$$ так как $${\mathop \Pi\limits_{n = 0}^{N + 2}} \left|{(k + \alpha - n)}\right| \sim (N + 2)$$!

При$$t \in \left[{t_N + 2\tau , t_N + 3\tau }\right]$$ $$$ \left|{R_N (t)}\right| \le \frac {(N + 2)(N + 3)}{2!} \tau ^{N + 1} \max\limits_{\xi \in [t_0, t_N + 3\tau ]} \left|{f^{(N + 1)}(\xi )}\right| $,$$ и так далее.

Видно, что ошибка экстраполяции растет быстро, но не сразу: экстраполяция допустима на интервалах $$\sim O(\tau ).$$

6.6. Интерполяционный полином в форме Ньютона

6.6.1. Разделенные и конечные разности

Определение. Пусть задана система узлов $$\left\{{t_n}\right\}_{n = 0}^{N}, t_n \in \left[{a, b}\right], t_0 = a, t_N = b.$$

Разделенные разности нулевого порядка в точке ti совпадают со значениями функции f(ti) ;

Разности первого порядка определяются для двух точек ti, ti + 1 равенством

$$$ f(t_i , t_{i + 1}) = \frac{{f(t_{i + 1}) - f(t_i)}}{{t_{i + 1} - t_i}}, $$$

разности второго порядка — для трех точек ti, ti + 1, ti + 2

$$$ f(t_i , t_{i + 1}, t_{i + 2}) = \frac{f(t_{i + 1}, t_{i + 2}) - f(t_i , t_{i + 1})}{t_{i + 2} - t_i}, $$$

разности порядка k — для k + 1 точки по рекуррентной формуле

$$$ f(t_i, t_{i + 1}, \ldots , t_{i + k}) = \frac{f(t_{i + 1}, \ldots , t_{i + k}) - f(t_i, \ldots , t_{i + k - 1})}{t_{i + k} - t_i}. $$$

Методом математической индукции можно показать, что

а) $$$ f(t_i, t_{i + 1}, \ldots , t_{i + k}) = \sum\limits_{j = 0}^{k}{\frac{f(t_{i + j})}{\prod\limits_{\substack{r = 0 \\ r \ne {i + j}}}^{k}{(t_{i + j} - t_{i + r})}}} $;$$

б) существует точка $$\xi \in \left[{a, b}\right]$$ такая, что $$k!f(t_{i}, t_{i + 1}, …, t_{i + k}) = f^{(k)}(\xi ).$$

Отсюда следует, что разделенная разность есть симметричная функция своих аргументов ti, ..., ti + k и она не изменяется при их перестановке.

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

t0	f(t0)			$$\ldots $$
		f(t0, t1)		$$\ldots $$
t1	f(t1)		f(t0, t1, t2)	$$\ldots $$
		f(t1, t2)		$$\ldots $$
t2	f(t2)		$$\ldots $$	$$\ldots $$	f(t0, …, tn)
$$\ldots $$	$$\ldots $$	$$\ldots $$	f(tn - 2, tn - 1, tn)	$$\ldots $$
		f(tn - 1, tn)		$$\ldots$$
tn	f(tn)			$$\ldots $$

Пусть сетка — равномерная. Тогда конечной разностью первого порядка функции f(t) в точке tk с шагом $$\tau$$ называют величину $$\Delta f_{k} = f_{k + 1} - f_{k},$$ где fk = f(tk), второго порядка — величину

$$\Delta ^{2}f_{k} = \Delta f_{k + 1} - \Delta f_{k} = f_{k + 2} - 2f_{k + 1} + f_{k},$$

третьего

$$\Delta ^{3}f_{k} = \Delta f_{k + 3} - 3f_{k + 2} + 3f_{k + 1} - f_{k} = f_{k + 3} - 3f_{k + 2} + 3f_{k + 1} - f_{k},$$

четвертого

$$\Delta ^{4}f_{k} = f_{k + 4} - 4f_{k + 3} + 6f_{k + 2} - 4f_{k + 1} + f_{k},$$

причем

$${\Delta}^{n} f_k = {\Delta}^{n - 1} f_{k + 1} - {\Delta}^{n - 1} f_k, k \ge 1, {\Delta}^0 f_k = f_k.$$

Методом математической индукции доказывается формула

$${\Delta}^{n} f_k = \sum\limits_{s = 0}^{n}{(- 1)^{n - 1} \cdot C_n^{s} f_{k + s}, } $$

где $$C_i^{s} = \frac{i!}{s!(i - s)!}$$ - биномиальные коэффициенты.

Нетрудно показать, например, используя формулу Лагранжа, что существует точка $$\xi \in \left[{a, b}\right]$$ такая, что $$\tau ^{n} f^{(n)} (\xi ) = {\Delta}^{n} f_k , t_k \in \left[{a, b}\right],$$ поэтому в вычислительных методах используется приближенная формула

$$$ f^{(n)} (t) \approx \frac{{\Delta}^{n} f_k}{\tau ^{n}}, $$$

аналогичная тем формулам численного дифференцирования, что были получены в первой лекции методом неопределенных коэффициентов.

Заметим, что введенные конечные разности называют "разностями вперед". Аналогично можно ввести "разности назад":

$$\Delta f_{k} = f_{k} - f_{k - 1}, \Delta ^{2}f_{k} = \Delta f_{k} - \Delta f_{k - 1} = f_{k} - 2f_{k - 1} + f_{k - 2}, \dots \\ \Delta ^{n}f_{k} = \Delta ^{n - 1}f_{k} - \Delta ^{n - 1}f_{k - 1}$$

и центральные разности

$$\Delta f_{k} = f_{k + 1/2} - \Delta f_{k - 1/2}, \Delta ^{2}f_{k} = \Delta f_{k + 1/2} - \Delta f_{k - 1/2} = f_{k + 1} - 2f_{k} + f_{k - 1}, \dots , \\ \Delta ^{n}f_{k} = \Delta ^{n - 1}f_{k + 1/2} - \Delta ^{n - 1}f_{k - 1/2}.$$

Иногда для обозначения первых конечных разностей вперед и назад используют обозначения $$\Delta ^{+}f_{k}, \Delta ^{-}f_{k}.$$

6.6.2. Интерполяционный полином в форме Ньютона

Интерполяционный полином может быть записан с использованием введенных выше разделенных разностей. Такая форма его записи называется интерполяционным полиномом в форме Ньютона. Полином имеет вид

Nn(t) = f(t1) + f(t1, t2)(t - t1) + ... + f(t1, ..., tn + 1)(t - t1) ... (t - tn).

То, что это интерполяционный полином, может быть доказано, например, методом математической индукции. Отметим, что полином в форме Ньютона напоминает ряд Тейлора, а остаточный член интерполяционного полинома — остаточный член этого ряда. Достоинством записи интерполянта в форме Ньютона является то, что для повышения порядка полинома нет необходимости в его полной перестройке; достаточно лишь добавить к уже полученному выражению еще одно или несколько слагаемых. С помощью разделенных разностей можно оценивать погрешность интерполяции. Читатели могут предложить способ контроля точности вычислений, основанный на использовании разделенных разностей.

6.7. Многочлены Чебышёва и минимизация остаточного члена интерполяции

Многочленом Чебышева первого рода называется функция Tn(t) = cos (n arccos t), где $$t \in \left[{- 1, 1}\right], n = 0, 1, \ldots $$

Убедимся в том, что функция Tn(t) действительно является многочленом. При n = 0 и n = 1 имеем T0(t) = 1, T1(t) = t.

Положив $$\theta = \arccos t,$$ получим $$T_{1}(t) = \cos \theta , T_{n}(t) = \cos n \theta,$$ $$T_{n - 1}(t) = \cos (n - 1)\theta , T_{n + 1}(t) = \cos (n + 1)\theta.$$ По формуле суммы косинусов $$\cos (n + 1)\theta + \cos (n - 1)\theta = 2\cos \theta \cos n\theta,$$ и справедливо рекуррентное соотношение Tn + 1(t) + Tn - 1(t) = 2T1(t)Tn(t), или Tn + 1(t) = 2t Tn(t) - Tn - 1(t). Отсюда следует вид записи полиномов Чебышева: T2 = 2t2 - 1, T3(t) = 4t3 - 3t, T4(t) = 8t4 - 8t2 + 1 и так далее. Функции Tn(t) являются многочленами степени n со старшим членом 2n - 1tn.

Введем также нормированные многочлены Чебышева $$$ \bar {T_n} (t) = \frac{T_n (t)}{2^{n - 1}} $.$$

Нули многочлена Чебышева находятся из очевидного уравнения Tn(t) = cos (n arccos t) = 0, откуда $$$ t_m = \cos \left({\frac{{2m - 1}}{n}\pi }\right), m = 1, 2, \ldots n, t \in [{- 1, 1}] $.$$ Для произвольного отрезка [a, b] нули полинома Чебышева получаются очевидным линейным преобразованием, выражения для них будут $$$ t_m = \frac{a + b}{2} + \frac{b - a}{2}\cos \left({\frac{2m - 1}{2n}\pi }\right), m = 1, 2, \ldots n $.$$ Легко отыскиваются также точки экстремумов полинома Чебышева, для них | Tn(t) | = 1 и на отрезке $$t \in [- 1, 1]$$ точки экстремумов есть $$$ t_m = \cos (\frac{m}{n}\pi ), m = 1, 2, \ldots n $.$$

Нас интересует решение следующей задачи на минимакс: найти

$$\min\limits_{\left\{t_n\right\}_{n = 0}^{N}}\left\{\max\limits_{t \in [- 1, 1] } \left|{\mathop \Pi\limits_{n = 0 }^{N}(t - t_n)}\right| \right\}$$

чтобы путем выбора узлов сетки минимизировать остаточный член интерполяции. Эта задача была решена П.Л.Чебышевым.

Теорема. (Чебышева (без доказательства)) Среди всех многочленов степени $$n \ge 1,$$ со старшим коэффициентом an равным единице, наименьшее уклонение от нуля, равное 21 - n, имеет нормированный полином Чебышева $$\bar {T_n} (t) = 2^{1 - n} T_n (t), t \in \left[{- 1, 1}\right].$$

Это свойство полиномов Чебышева, наименьшее уклонение от нуля, можно сформулировать по - другому: для любого полинома Pn(t) = tn + an - 1tn - 1 + ... + a0, отличного от $$$ \bar {T_n} (t) $$$ справедливо $$2^{1 - n} = \max\limits_{\left[{- 1, 1}\right]} \left|{\bar {T_n} (t)}\right| < \max\limits_{\left[{- 1, 1}\right]} \left|{P_n (t)}\right|, t \in \left[{- 1, 1}\right].$$

Если в качестве интерполяционных узлов выбрать нули полинома Чебышева, то произведение $${\mathop \Pi\limits_{j = 0}^{N + 1} (t - t_j)},$$ а также Rn(t) будут наименее уклоняющимися от нуля.

6.8. Обусловленность задачи интерполяции. Постоянная Лебега

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

$$L_N (t) = \sum\limits_{n = 0}^{N}{f_n}\varphi_n^{N} (t) + \sum\limits_{n = 0}^{N}{\delta f_n} \varphi_n^{N} (t),$$

причем слагаемое

$$\Delta_N (t, \delta f) = \sum\limits_{n = 0}^{N}{\delta f_n} \varphi_n^{N} (t)$$

характеризует чувствительность к ошибкам начальных данных и ошибкам вычислений. Нас интересует оценка

$$\max\limits_{t \in \left[{a, b}\right]} \left|{\Delta_N (t, \delta f)}\right| \le l_N \delta ; \delta = \max\limits_{t \in \left[{a, b}\right]} \left|{\delta f_n}\right|, l_N = \max\limits_{t \in [a, b]}\sum\limits_{n = 0}^{N} \left|\varphi_n^{N} (t)\right|,$$

коэффициент lN называется постоянной Лебега.

Введем в рассмотрение еще один объект. Пусть $$\sum {\left|{\varphi_i^{N} \left(x\right)}\right|}$$ - сумма модулей всех базисных функций. Обозначим ее $${L}\left(x\right) = \sum {\left|{l_i^{N} \left(x\right)}\right|}$$ - функция Лебега (сетки). Тогда константа Лебега $$l_N = \sup\limits_{x \in \left[{a, b}\right]}{L}\left(x\right).$$

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

Если рассматривать оператор интерполяции как оператор проекции (проектор), переводящий элемент одного пространства (пространства сеточных функций) в другое (пространство непрерывно дифференцируемых функций), то постоянная Лебега есть норма такого оператора проекции. Подробнее об этом в [6.7].

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

Приведем (без доказательства) примерные оценки роста постоянной Лебега в зависимости от числа узлов сетки. Константа Лебега растет примерно как $$l_{N} \sim 2^{N}$$ для равномерной сетки и $$l_{N} \sim ln (N)$$ для сетки с чебышевским набором узлов. Доказано, что рост константы Лебега для последней сетки асимптотически стремится к минимально возможному, и сетка с чебышевскими узлами близка к оптимальной для задач интерполяции.

6.9. Интерполяция с кратными узлами

Определение. Пусть в узлах сетки $$\left\{{t_n}\right\}_{n = 0}^{M}$$ заданы не только значения функции f(tn), но и значения ее производных f'(tn), f''(tn), ..., $$f^{(k_n - 1)}(t_n).$$ В этом случае узел tn называется кратным, а число kn, равное количеству заданных значений производных в n узле — кратностью узла.

Доказывается теорема о существовании единственного полинома PN(t), удовлетворяющего условиям

$$\begin{gather*} P_N (t_n) = f_n, P^{\prime}_N (t_n) = f^{\prime}_n, \ldots , P_N^{(k_n - 1)} (t_n) = f_n^{k_n - 1}, \\ N = k_0 + k_1 + \ldots + k_M - 1. \end{gather*} $$

Такой полином называется полиномом с кратными узлами. Отметим два частных случая.

в точке t = t0 заданы f0, f'0, ..., f0(N) (M = 0, k0 = N + 1).

Тогда многочлен PN(t), удовлетворяющий этим условиям, может быть записан как

$$$ P_N (t) = \sum\limits_{i = 0}^{N}{f^{(i)} (t_0 )}\frac{{(t - t_0)}^i}{i!}. $$$

Это — ряд Тейлора, который является интерполянтом с кратным узлом в точке t = t0 кратности N + 1.

Пусть на концах отрезка [t0, t1] заданы значения f0, f1, f'0, f'1 (M = 1, k0 = 2, k1 = 2, N = 3). Тогда P3(t0) = f0, P'3(t0) = f'0, P3(t1) = f1, P'3(t1) = f'1, а интерполянт имеет вид

$$\begin{gather*} P_3 (t) = f_0 \frac{{(t_1 - t)^2 \left[{2(t - t_0 ) + \tau }\right]}}{{\tau ^3 }} + {f^{\prime}}_0 \frac{{(t_1 - t)^2 (t - t_0 )}}{{\tau ^2 }} + \\ + f_1 \frac{{(t - t_0 )^2 \left[{2(t_1 - t) + \tau }\right]}}{{\tau ^3 }} + f^{\prime}_1 \frac{{(t - t_0)^2 (t - t_1 )}}{{\tau ^2 }}, \end{gather*}$$

здесь $$\tau = t_{1} - t_{0}.$$

Такой многочлен называется кубическим интерполяционным многочленом Эрмита.

Теорема (без доказательства). Пусть f(t) имеет N + 1 ограниченную производную на отрезке [a, b]. Тогда погрешность интерполяционного многочлена Эрмита степени N выражается формулой

$$$ R_N (t) = \frac{{f^{(N + 1)} (\xi )}}{{(N + 1)!}}(t - t_0 )^{k_0 } (t - t_1)^{k_1 } \ldots (t - t_M)^{k_M}, $$$

где tn — интерполяционные узлы, n = 0, ..., M, ki - кратность i узла, $$\xi \in \left[{a, b}\right], N = k_0 + k_1 + \ldots + k_M - 1.$$

Поставим теперь следующую задачу: построить кусочно - кубическую интерполирующую функцию, непрерывную на отрезке [a, b] со своими двумя первыми производными.

Обозначим такую функцию S(t) ; значения производных в узлах tn обозначим mn = S'(tn). Если задать в узлах tn, tn + 1 значение функции и ее первой производной, то получим эрмитов кусочно - кубический полином

$$\begin{gather*} S(t) = \frac{{(t_{n + 1} - t)^2 \left[{2(t - t_n) + \tau }\right]}}{{\tau ^3 }}f_n + \frac{{(t - t_n)^2 \left[{2(t_{n + 1} - t) + \tau }\right]}}{{\tau ^3 }}f_{n + 1} + \\ + \frac{{(t_{n + 1} - t)^2 (t - t_n)}}{{\tau ^2 }}m_n + \frac{{(t - t_n)^2 (t - t_{n + 1})}} {{\tau ^2 }}m_{n + 1}, \end{gather*}$$

или

$$S(z) = f_{n}(1 - z)^{2}(1 + 2z) + f_{n + 1}z^{2}(3 - 2z) + m_{n}\tau _{n}z(1 - \tau )^{2} - m_{n + 1}\tau _{n}z^{2}(1 - z),$$

где $$\tau _{n} = t_{n + 1} - t_{n}, z = (t - t_{n})/\tau _{n}, m_{n} = S'(t_{n}).$$

6.9.1. Замечание о тригонометрической интерполяции

Для периодической функции f(t) с периодом T естественно строить приближение с использованием функций $$$ \varphi_n(t) = a_n \cos {\frac{\pi nt}{T}} + b_n \sin{\frac{\pi nt}{T}} $$$ Тригонометрическая интерполяция состоит в замене f(t) тригонометрическим многочленом $$$ F_N (t) = \sum\limits_{n = 0}^{N}{\varphi_n (t) = a_0 + \sum\limits_{n = 1}^{N}{(a_n \cos \frac{{\pi nt}}{T} + b_n \sin \frac{\pi nt}{T})}} $,$$ коэффициенты которого находятся при решении СЛАУ FN(tk) = f(tk), k = 1, ..., 2N + 1, t2N + 1 - t0 = T, здесь $$\left\{{t_k}\right\}_{k = 0}^{2N + 1}$$ — последовательность узлов интерполяции.

6.10. Кусочно - многочленная глобальная интерполяция (сплайны)

Определение. Пусть на отрезке [a, b] задана система узловых точек $$\{t_n\}_{n = 0}^{N - 1}.$$ Сплайном Sm(t) называется определенная на [a, b] функция, имеющая l непрерывных производных и являющаяся на каждом интервале (tn - 1, tn) многочленом степени m.

Определение. Дефектом сплайна называется разность d = m - l между степенью сплайна и показателем его гладкости l.

Замечание. Для сплайнов также используется обозначение Sm, d(t). Если сплайн строится так, чтобы выполнялись условия Sm(tn) = f(tn), где f(t) — интерполируемая функция, то он называется интерполяционным сплайном. В соответствии с определением, кусочно - линейная функция является интерполяционным сплайном первой степени дефекта 1, кусочно - квадратичная функция с первой непрерывной производной — интерполяционным сплайном второй степени дефекта 1. Наиболее известным в приложениях является интерполяционный кубический сплайн дефекта 1 (естественный сплайн ), который будем обозначать S(t).

Определение. Кубическим сплайном дефекта 1, интерполирующим на отрезке [a, b] заданную функцию f(t), называется функция S(t), удовлетворяющая следующим условиям:

S(tn) = f(tn) — условие интерполяции в узлах сетки $$\left\{{t_n}\right\}_{n = 0}^{N - 1}.$$

$$S(t) \in C^2 [a, b],$$ т.е. является непрерывной вместе с двумя первыми производными.

На каждом отрезке [tn, tn + 1], S(t) является кубическим многочленом; n = 0, ..., N - 1.

На краях отрезка [a, b] заданы краевые условия. Наиболее часто употребляются следующие:

S'(a) = f'(a), S'(b) = f'(b) ;

S''(a) = f''(a), S''(b) = f''(b) ; часто полагают S''(a) = S''(b) = 0 ;

S(a) = S(b), S'(a) = S'(b) ; эти условия называются периодическими, т.е. интерполируемая функция является периодической с периодом b - a.

Покажем, что эта задача имеет единственное решение.

Теорема. Интерполяционный кубический сплайн S(t), удовлетворяющий условиям 1 — 3 и одному из краевых условий 4, существует и единственен.

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

Пусть S(z) — эрмитов кубический многочлен, который на каждом отрезке [ tn, tn + 1 ], n = 0, …, N - 1, представлен как

$$S(z) = f_n (1 - z)^2 (1 + 2z) + f_{n + 1} \cdot z^2 (3 - 2z) + m_n \tau_n z(1 - z)^2 - m_{n + 1} \tau_n z^2 (1 - z),$$

где $$\tau _{n} = t_{n + 1} - t_{n},$$ $$z = (t - t_{n})/\tau _{n}, m_{n} = S'(t_{n}).$$ Тогда

$$\begin{gather*} S^{\prime\prime}(t) = \frac{{(f_{n + 1} - f_n)(6 - 12z)}}{{\tau_n^2 }} + m_n \frac{{6z - 4}} {{\tau_n}} + m_{n + 1} \frac{{6z - 2}}{{\tau_n}}, \\ S^{\prime\prime}(t_n + 0) = 6\frac{{f_{n + 1} - f_n}}{{\tau_n^2 }} - \frac{{4m_n}}{{\tau_n}} - \frac{{2m_{n + 1}}}{{\tau_n}}, \\ S^{\prime\prime}(t_n - 0) = - 6\frac{{f_n - f_{n - 1}}}{{\tau_{n - 1}^2 }} + \frac{{2m_{n - 1}}} {{\tau_{n - 1}}} + m_{n + 1} \frac{{4m_n}}{{\tau_{n - 1}}}. \end{gather*} $$

Условие непрерывности второй производной S''(tn + 0) = S''(tn - 0) будет

$$\begin{gather*} r_n m_{n - 1} + 2m_n + s_n m_{n + 1} = c_n , \\ c_n = 3\left({s_n \frac{{f_{n + 1} - f_n}}{{\tau_n}} + r_n \frac{{f_n - f_{n - 1}}}{{\tau_{n - 1}}}}\right), s_n = \frac{{\tau_{n - 1}}}{{\tau_{n - 1} + \tau_n}}, r_n = 1 - s_n, \\ n = 1, \ldots , N - 1 \end{gather*} $$

После добавления краевых условий получаем систему из N + 1 уравнение с N + 1 неизвестным mn. Для краевых условий первого типа (заданы первые производные) система выглядит как

m0 = f'0, 
rn mn - 1 + 2mn + snmn + 1 = cn , 
mn = f'n

Для условий второго типа (заданы вторые производные)

$$\begin{gather*} 2m_0 + m_1 = 3\frac{{f_1 - f_0 }}{{\tau_0 }} - \frac{{\tau_0 }}{2}f^{\prime\prime}_0, \\ r_n m_{n - 1} + 2m_n + s_nm_{n + 1} = c_n , \\ m_{N - 1} + 2m_N = 3\frac{{f_N - f_{N - 1}}}{{\tau_{N - 1}}} + \frac{{\tau_{N - 1}}}{2}f^{\prime\prime}_N \end{gather*}$$

Аналогично получается СЛАУ для третьего типа краевых условий.

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

Приведем еще одно доказательство этой же теоремы.

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

Рассмотрим неравномерную сетку: $$t_{n} - t_{n - 1} = \tau _{n - 1},$$ $$t_{n + 1} - t_{n} = \tau _{n}.$$ В узлах сетки определены значения функции: fn - 1, fn , fn + 1. Пусть mn — значение второй производной в точке tn (пока неизвестное!). На отрезке [tn, tn + 1] для второй производной кусочно - кубического сплайна имеем

$$$ {S^{\prime\prime}_{tt} = \frac {1}{\tau_n}(m_n (t_{n + 1} - t) + m_{n + 1} (t - t_n)).} $$$

Так как сплайн — полином третьей степени, то его вторая производная — линейная функция. Интегрируем (6.2) по t, получаем (на отрезке [ tn, tn + 1] )

$$$ S^{\prime}_{t} = \frac {1}{\tau_n} \left(m_{n + 1}\frac{{(t_{n + 1} - t)}^2}{2} - m_n \frac{{(t - t_n)}^2}{2}\right) + A_n. $$$

bi - 1Si - 1 + biSi + bi + 1Si + 1= fi,

axA + byA + c = fA, 
axB + byB + c = fB, 
axD + byD + c = fD,

Ln(x) = a0 + a1t + ... + antn, 
Ln(x) = a0 + t(a1 + t(a2 + t(... (an - 1 + cnt) ...))) (схема Горнера)

Pn(t) = a0 + a1(t - t0 ) + ... + an(t - t0) ... (t - tN - 1).

f(t) - LN(t) = (t - t0) ... (t - tN)f(t, t0, ..., tN),

a0 + a1t1 - d0f1 = t1f1, 
a0 + a1t2 - d0f2 = t2f2, 
a0 + a1t3 - d0f2 = t3f3.