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

Материал из MachineLearning.

Версия от 13:39, 17 октября 2008; Leo (Обсуждение | вклад)

(разн.) ← Предыдущая | Текущая версия (разн.) | Следующая → (разн.)

Содержание

1 Постановка задачи
- 1.1 1. Полином в каноническом виде
- 1.2 2. Способ вычисления полинома в точке
2 Анализ метода
3 Смотри также

Постановка задачи

Пусть задана функция f(x) на отрезке [a,b]. Задача интерполяции состоит в построении функции g(x), совпадающей с f(x) в некотором наборе точек x₀, x₁,...,x_n из отрезка [a,b]. Эти точки называются узлами интерполяции. Также должно выполняться условие: g(x_k) = y_k, k=0,...,n, где y_k = f(x_k).

1. Полином в каноническом виде

Известно, что любая непрерывная на отрезке [a,b] функция f(x) может быть хорошо приближена некоторым полиномом P_n(x). Справедлива следующая Теорема (Вейерштрасса): Для любого $\eps$ >0 существует полином P_n(x) степени $n=n(\eps)$ , такой, что $max_{x \in [a,b]}|f(x)-P_n(x)|<\eps$

В качестве аппроксимирующей функции выберем полином степени n в каноническом виде:

$f(x)=P_n(x)=c_0+c_1x+c_2x^2+ \ldots + c_nx^n$

Коэффициенты полинома $c_i$ определим из условий Лагранжа $P_n(x_i)=y_i$ , $i=1, \ldots, n$ , что с учётом предыдущего выражения даёт систему уравнений с n+1 неизвестными:

$\begin{matrix} c_0 + c_1x_0 + c_2x_0^2 + \ldots + c_nx_0^n = y_0 \\ c_0 + c_1x_1 + c_2x_1^2 + \ldots + c_nx_1^n = y_1 \\ \ldots \ldots \ldots \ldots \ldots \\ c_0 + c_1x_n + c_2x_n^2 + \ldots + c_nx_n^n = y_n \end{matrix}$

Обозначим систему таких уравнений символом (*) и перепишем её следующим образом:

$\sum_{p=0}^n c_i^p = y_i, \quad i=1, \ldots, n$

или в матричной форме: $\mathbf{Ac}=\mathbf{y},$ где $\mathbf{c}$ - вектор-столбец, содержащий неизвестные коэффициенты $c_i$ , $\mathbf{y}$ - вектор-столбец, составленный из табличных значений функции $y_i$ , а матрица $\mathbf{A}$ имеет вид:

$\mathbf{A} = \begin{pmatrix} 1 & x_0 & x_0^2 & \ldots & x_0^n \\ 1 & x_1 & x_1^2 & \ldots & x_1^n \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ 1 & x_n & x_n^2 & \ldots & x_n^n \\ \end{pmatrix}$

Система линейных алгебраических уравнений (*) относительно неизвестных $c_i$ иметь единственное решение, если определитель матрицы $\mathbf{A}$ отличен от нуля.

Определитель матрицы $\mathbf{A}$ называют определителем Вандермонда, его можно вычислить по следующей формуле:

$\mathbf{detA} = \prod_{i,j=0 \\ i \neq j }^n (x_i - x_j) \neq 0$

Число узлов интерполяционного полинома должно быть на единицу больше его степени. Это понятно из интуитивных соображений: через 2 точки можно провести единственную прямую, через 3 - единственную параболу и т.д. Но полином может получиться и меньшей степени. Т.е. если 3 точки лежат на одной прямой, то через них пройдёт единственный полином первой степени (но это ничему не противоречит: просто коэффициент при старшей степени равен нулю).

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

[Рис.1 Зависимость числа обусловленности матрицы от количества узлов интерполяции]

Из-за плохой обусловленности матрицы рекомендуется применять другие методы интерполяции (например, метод Лагранжа). При этом важно понимать, что при теоретическом принемении различных методов они приводят к одинаковому результату, т.е. мы получим один и тот же полином.

Однако при практической реализации мы получим полиномы различной точности аппроксимации из-за погрешности вычислений аппаратуры.

2. Способ вычисления полинома в точке

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

Первый способ. Можно посчитать значение a₁x и сложить с a₀. Далее найти a₂x², сложить с полученным результатом, и так далее. Таким образом, на j-ом шаге вычисляется значение a_jx^j и складывается с уже вычисленной суммой $a_0 + a_1x + \ldots + a_{j-1}x^{j-1}$ .

Вычисление значения a_jx^j требует j операций умножения. Т.е. для подсчёта многочлена в заданной точке требуется (1 + 2 + ... + n) = n(n+1)/2 операций умножения и n операций сложения. Всего операций в данном случае: Op₁ = n(n+1)/2 + n.

Второй способ. Полином можно также легко вычислить с помощью так называемой схемы Горнера: $P_n(x) = (\ldots ((a_nx + a_{n-1})x + a_{n-2})x + \ldots)x + a_0$

Для вычисления значения во внутренних скобках a_nx + a_n-1 требуется одна операция умножения и одна операция сложения. Для вычисления значения в следующих скобках (a_nx + a_n-1)x + a_n-2 требуется опять одна операция умножения и одна операция сложения, т.к. a_nx + a_n-1 уже вычислено, и т.д.

Тогда в этом способе вычисление P_n(x) потребует n операций умножения и n операций сложения, т.е. сложность вычислений Op₂ = n+n = 2n. Ясно, что Op₂ << Op₁.

Анализ метода

1. Сложность вычислений

2. Погрешность интерполяции

Предположим, что на отрезке интерполирования [a,b] функция f(x) n раз непрерывно-дифференцируема. Погрешность интерполяции складывается из погрешности самого метода и ошибок округления.

Ошибка приближения функции f(x) интерполяционным полиномом n-ой степени P_n(x) в точке x определяется разностью: R_n(x) = f(x) - P_n(x).

Погрешность R_n(x) определяется следующим соотношением:
$R_n(x) = \frac{f^{n+1}(\xi)}{(n+1)!}\omega_n(x)$

Здесь $f^{n+1}(\xi)$ - производная (n+1)-го порядка функции f(x) в некоторой точке $\xi \in [a,b],$ а функция $\omega_n(x)$ определяется как
$\omega_n(x)=(x-x_0)(x-x_1)\ldots (x-x_n).$

Если максимальное значение производной fⁿ⁺¹(x) равно $M_{n+1} = \sup_{x \in [x_0, x_n]} \left| f^{n+1}(x) \right| ,$ то для погрешности интерполяции следует оценка: $\left| R_n(x) \right| = \frac{M_{n+1}}{(n+1)!}\omega_n(x).$

3. Выбор узлов интерполяции

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

Введём следующее определение: полиномом Чебышева называется многочлен вида

T_k(x) = cos(k arccos x), |x|≤1.

Известно (см. ссылки литературы), что если узлы интерполяции x₀, x₁,...,x_n являются корнями полинома Чебышева степени n+1, то величина $\max_{x \in [a,b]} \left| \omega_{n+1}(x) \right|$ принимает наименьшее возможное значение по сравнению с любым другим выбором набора узлов интерполяции. Очевидно, что в случае k = 1 функция T₁(x), действительно, является полиномом первой степени, поскольку T₁(x) = cos(arccos x) = x. В случае k = 2 T₂(x) тоже полином второй степени. Это нетрудно проверить. Воспользуемся известным тригонометрическим тождеством: cos2θ = 2cos²θ - 1, положив θ = arccos x. Тогда получим следующее соотношение: T₂(x) = 2x² - 1. С помощью тригонометрического тождества cos(k + 1)θ = 2cosθcoskθ - cos(k - 1) легко показать, что для полиномов Чебышева справедливо реккурентное соотношение: T_k+1(x) = 2xT_k(x) - T_k-1(x)

При помощи данного соотношения можно получить формулы для полиномов Чебышева любой степени. Корни полинома Чебышева легко находятя из уравнения: T_k(x) = cos(k arccos x) = 0. Получаем, что уравнение имеет k различных корней, расположенных на отрезке [-1,1]: $x_m = \cos \frac{2m+1}{2k}\pi, \, m = 0,1, \cdots, k-1,$ которые и следует выбирать в качестве узлов интерполирования. Нетрудно видеть, что корни на [-1,1] расположены симметрично и неравномерно - чем ближе к краям отрезка, тем корни расположены плотнее. Максимальное значение модуля полинома Чебышева равно 1 и достигается в точках $\cos \frac{m}{k}\pi.$ Если положить $\omega_k(x) = \frac1{2^{k-1}}T_k(x),$ то для того, чтобы коэффициент при старшей степени полинома ω_k(x) был равен 1, $\max_{x \in [-1,1]}\omega_k(x) = \frac1{2^{k-1}}.$ Известно, что для любого полинома p_k(x) степени k с коэффициентом, равным единице при старшей производной верно неравенство $\max_{x \in [-1,1]}p_k(x) \geq \frac1{2^{k-1}},$ т.е. полиномы Чебышева являются полиномами, наименее уклоняющимися от нуля.