Интерполяционные многочлены

Пусть даны значения $f_k$ функции $f$ в конечном наборе точек $x_k$, $k=0..K$. Задача интерполяции заключается в восстановлении значений функции в промежуточных точках. Так как существует бесконечное число функций, принимающих данные значения, то восстановить значения однозначно невозможно. Однако если предположить, что функция имеет некий специальный вид, то оказывается возможным однозначно восстановить в промежуточных точках, т.е. построить интерполяционную функцию. Из функционального анализа известно, что любую непрерывную функцию на отрезке можно приблизить многочленом, поэтому в качестве интерполяционной функции естественно взять многочлен (однако это не обязательно лучший выбор). Интерполяционный многочлен степени $K$ можно записать в виде: \(P(x)=a_0+a_1x+a_2x^2+\ldots+a_Kx^K,\) где $a_k$ - коэффициенты многочлена. Так как существует единственный многочлен степени $K$, проходящий через $K+1$ различных точек, то по данным $x_k$ и $f_k$ можно однозначно восстановить многочлен, т.е. найти коэффициенты $a_k$. Мы можем исходить из того, что значения многочлена в точках $x_n$ равны $f_n$, что приводит нас к следующей системе: \(\begin{cases} f_0=P(x_0)=a_0+a_1x_0+a_2x_0^2+\ldots+a_Kx_0^K,\\ f_1=P(x_1)=a_0+a_1x_1+a_2x_1^2+\ldots+a_Kx_1^K,\\ \cdots\\ f_K=P(x_K)=a_0+a_1x_K+a_2x_K^2+\ldots+a_Kx_K^K.\\ \end{cases}\) В матрично виде систему можно записать как $F=WA$, где $A=(a_0,\ldots,a_K)^T$ - столбец неизвестных, $F=(f_0,\ldots,f_K)^T$ - столбец свободных членов, а $W$ - матрица системы, составленная из степенй чисел $x_k$: \(W=\begin{pmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^K \\ 1 & x_1 & x_1^2 & \cdots & x_1^K \\ \cdots & \cdots & \cdots & \cdots & \cdots \\ 1 & x_K & x_K^2 & \cdots & x_K^K \end{pmatrix}.\) Матрица такого вида называется матрицей Вандермонда. Для нахождения коэффициентов многочлена достаточно решить эту линейную систему, однако эта задача часто плохо обусловленна, что не позволяет нам точно восстановить коэффициенты. Также решение системы требует больше вычислительных ресурсов, чем необходимо.

Многочлены Лагранжа

Представленный выше вид многочлена $P$ в виде суммы степеней аргумента $x$ с коэффициентами $a_k$ называется разложением по одночленам. Такой вид чаще всего встречается в литературе как определение многочлена. Однако для вычислений могут оказаться более удобны другие представления многочлена, например, если многочлен задан своими нулями $\xi_k$, $k=1..K$, то его удобнее представить в виде \(P(x)=a_K(x-\xi_1)\cdots(x-\xi_K).\) При построении интерполяционного многочлена по точкам $(x_k,f_k)$ удобно представить многочлен в виде интерполяционного многочлена Лагранжа: \(P(x)=\sum_{k=0}^K y_k l_k(x),\) где се многочлены $l_k$ имеют степень $K$ и зависят только от узлов интерполяции $x_k$: \(l_k(x)=\prod_{j\neq k}\frac{x-x_k}{x_j-x_k}.\) Многочлены $l_k$ выбраны так, чтобы $l_k(x_k)=1$ и $l_k(x_j)=0$ для всех $j\neq k$. Благодаря этому свойству легко находятся коэффициенты $y_k$ разложения по $l_k$: \(P(x_j)=\sum_{k=0}^K y_k l_k(x_j)=y_j l_j(x_j)=y_j,\) но так как по определение интерполяционного многочлена $f_j=P(x_j)$, то $y_j=f_j$. Это значит, что построение интерполяционного многочлена Лагранжа производится за постоянное время по данным $f_k$. Однако вычисление многочленов $l_k$ в точках может потребовать больше времени, чем вычисление одночленов. Чтобы ускорить расчеты, мы можем разложить $l_k(x)$ по степеням $x$ и проссумировать их с коэффициентами $y_k=f_k$.

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

Ортогональные многочлены

Многочлены Лагранжа позволяют легко построить многочлен, принимающий заданные значения в узловых точках. Однако часто нас интересуют не значения многочлена в точках, а его поведение на всем интервале, например, нас может интересовать многочлен, дающий наименьшую среднеквадратическую ошибку приближения заданной функции на неком интервале. В этом случае удобнее воспользоваться ортогональными многочленами. Ортогональными многочленами называют набор многочленов $p_n$, $n=0\ldots$, таких что:

1. Многочлен $p_n$ имеет степень $n$.
2. Многочлены ортогональны, т.е. $\langle p_n|p_m\rangle=0$ для $n\neq m$, относительно некоторого скалярного произведения:

Функция $\rho$ называется весом скалярного произведения, числа $a$ и $b$ задают интервал $[a,b]$, на котором рассматриваются многочлены.

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

Свойство 1. Ортогональные многочлены можно задать с помощью реккурентного соотношения: \(p_{n+1}(x)=(Ax+B)p_n(x)+Cp_{n-1}(x),\) для некоторых констант $A$, $B$ и $C$ Это соотношение позволяет находить значения многочленов $p_n$ для всех $n$ в одной точке за линейное время по числу многочленов.

Свойство 2. Ортогональные многочлены образуют полный набор, т.е. любой многочлен $s(x)$ степени $N$ можно представить в виде \(s(x)=\sum_{n=0}^N \alpha_n p_n(x).\) Если в качестве $p_n$ взять одночлены $x^n$, то такое разложение превратится в обычное определение многочлена из учебников \(s(x)=\sum_{n=0}^N s_n x^n.\) Разложение по ортогональные многочленам может быть предпочтительнее, чем разложение по одночленам, так как может уменьшить погрешность вычислений. Сложность вычисления значения многочлена разложенного по ортогональным в точке линейна по степени многочлена, так как с помощью реккурентного соотношения можно за линейное время найти значения всех $p_n$. Более эффективный способ вычислений значения многочлена дает аналог схемы Горнера для ортогональных многочленов.

Свойство 3 Каждый ортогональный многочлен $p_n$ имеет ровно $n$ различный нулей $x_n$, все они простые и лежат на интервале $[a,b]$.

Многочлены Чебышева

Существует много видов ортогональных многочленов, но в нашем введение мы рассмотрим только один тип: многочлены Чебышева. Многочлены Чебышева $T_n$ заданы на симметричном интервале $[-1,1]$ и меньше других многочленов отличаются от нуля на этом интервале, поэтому часто используются для аппроксимации.

Многочлен Чебышева на интервале $[-1,1]$ можно задать через тригонометрические функции: \(T_n(x)=cos(n\cdot\arccos x),\) для $n=0,1,2,\ldots$. Эти многочлены также подчинены рекуррентному соотношению: \(T_0(x)=1,\quad T_1(x)=x,\quad T_{n+1}(x)=2xT_n(x)-T_{n-1}(x),\) которое удобнее для вычисления значений многочлена, чем определение через тригонометрию. Нули многочлена $T_n(x)$ расположены в точках \(\zeta_j^{(n)}=\cos\frac{\pi(k-1/2)}n.\) Многочлены Чебышева ортогональны относительно следующего скалярного произведения, совпадающего с обчным евклидовым скалярным произведением векторов, составленных из значений многочелнов в нулях $\zeta_j^{(n)}$: \begin{equation} \langle f_1|f_2\rangle=\sum_{j=1}^n f_1(\zeta_j^{(n)}f_2(\zeta_j^{(n)}). \end{equation} Указанная ортогональность позволяет легко находить коэффициенты разложения по многочленам Чебышева через это скалярное произведение. Действительно, предположим, что многочлен $p$ разложен по многочленам Чебышева с некоторыми коэффициентами $p_k$: \(p(x)=\sum_{k=1}^n p_k T_k(x).\) Умножая левую и правую часть равенства скалярно на $T_j$ и замечая, что в силу ортогональности многочленов Чебышева $\langle T_j|T_k\rangle=0$ для всех $j\neq k$, получаем следующее тождество: \(\langle T_l|p\rangle=\langle T_l|\sum_{k=1}^n p_k T_k\rangle =\sum_{k=1}^n p_k \langle T_l|T_k\rangle =p_l \langle T_l|T_l\rangle,\) из которого легко выразить коэффициенты разложения многочлена $p$ по многочленам Чебышева: \(p_l=\left(\sum_{j=1}^n p(\zeta_j^{(n)})T_l(\zeta_j^{(n)})\right)/\left(\sum_{j=1}^n T_l(\zeta_j^{(n)})^2\right).\) Выше мы предполагали, что $p$ является многочленом степени $n$, что позволило нам получить точное равенство в разложении. Если же $p$ произвольная функция, то построенный по предложенным формулам многочлен будет совпадать с функцией $p$ вообще говоря только в точках $\zeta_j^{(n)}$, т.е. мы получим интерполяционный многочлен. Так как существует единственный интерполяционный многочлен данной степени $n$, принимающий заданные значения, то этот же многочлен мы получим используя интерполяционный многочлен Лагранжа с узлами интерполяции в точках $\zeta_j^{(n)}$.

Приближение функции интерполяционным многочленом дает нулевую погрешность в узлах интерполяции, однако в остальных точках ошибка может быть большой. Метод Ланцоша дает другой способ построения разложения функции по полиномам Чебышева, который дает меньшую среднюю ошибку на всем интеравле $[-1,1]$. Рассмотрим сначала дифференциальное уравнение, решением которого является искомый логарифм \(x\frac{dy}{dx}-1=0,\quad y(1)=0.\) Так как логарифм не определен на всем интервале $[-1,1]$, то мы введем новую переменную $u$, такую что \(x=\frac{1+2u/3}{1-2u/3}.\) Изменение $u$ в интервале $[-1,1]$ соответствует изменению $x$ в интервале $(1/5, 5)$. Относительно переменной $u$ дифференциальное уравнение принимает вид \((1-4u^2/9)\frac{dy}{du}-4/3=0.\) Будем решать это уравнение, раскладывая производную решения по первым $N+1$ многочленам Чебышева: \(\frac{dy}{du}=\sum_{k=0}^Na_kT_k(x).\) Выражая производную из дифференциального уравнения, и пользуясь следующим свойством многочленов Чебышева: \(u^2T_k(u)=\frac14(T_{k-2}(u)+2T_k(u)+T_{k+2}(u)),\) мы преобразуем дифференциальное уравнение в алгебраическое \(\sum_{k=0}^{N+2}\tau_k T_k(u)=0,\) где $\tau_k$ константы, выражаемые через коэффициенты $a_0,\ldots,a_N$ разоложения. Разложение по Ланцошу будет получено, если все константы равны нулю, кроме двух последних: \(\tau_0=\tau_1=\cdots=\tau_N=0,\) что дает нам систему уравнений, из которой можно найти коэффициенты $a_k$. Зная $a_k$, мы знаем и $dy/du$. Интегрируя $dy/du$ мы находим разложение для $y(u)$ по многочленам Чебышева, воспользовавшись аналитическим выражением для первообразной многочлена Чебышева: \(\int T_k(u)du=\frac1{2(k+1)}T_{k+1}(u)-\frac1{2(k-1)}T_{k-1}(u) +\frac{\sin(\pi k/2)}{k^2-1}.\)