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

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

(Различия между версиями)
Перейти к: навигация, поиск
Строка 1: Строка 1:
 +
{{TOCright}}
== Постановка задачи ==
== Постановка задачи ==
-
Пусть задана функция ''f(x)'' на отрезке [''a,b'']. Задача интерполяции состоит в построении функции ''g(x)'', совпадающей с ''f(x)'' в некотором наборе точек ''x<sub>0</sub>, x<sub>1</sub>,...,x<sub>n</sub>'' из отрезка [''a,b'']. Эти точки называются узлами интерполяции. Также должно выполняться условие: ''g(x<sub>k</sub>) = y<sub>k</sub>, k=0,...,n,'' где ''y<sub>k</sub> = f(x<sub>k</sub>).''
+
Пусть задана функция ''f(x)'' на отрезке [''a,b'']. '''Задача интерполяции''' состоит в построении функции ''g(x)'', совпадающей с ''f(x)'' в некотором наборе точек ''x<sub>0</sub>, x<sub>1</sub>,...,x<sub>n</sub>'' из отрезка [''a,b'']. Эти точки называются узлами интерполяции. Также должно выполняться условие: ''g(x<sub>k</sub>) = y<sub>k</sub>, k=0,...,n,'' где ''y<sub>k</sub> = f(x<sub>k</sub>).''
=== Полином в каноническом виде ===
=== Полином в каноническом виде ===
Строка 27: Строка 28:
<tex>\sum_{p=0}^n c_i^p = y_i, \quad i=1, \ldots, n</tex>
<tex>\sum_{p=0}^n c_i^p = y_i, \quad i=1, \ldots, n</tex>
-
или в матричной форме: <tex>\mathbf{Ac}=\mathbf{y},</tex> где <tex>\mathbf{c}</tex> - вектор-столбец, содержащий неизвестные коэффициенты <tex>c_i</tex>, <tex>\mathbf{y}</tex> - вектор-столбец, составленный из табличных значений функции <tex>y_i</tex>, а матрица <tex>\mathbf{A}</tex> имеет вид:
+
или в [http://ru.wikipedia.org/wiki/Матрица_(математика) матричной] форме: <tex>\mathbf{Ac}=\mathbf{y},</tex> где <tex>\mathbf{c}</tex> - [http://ru.wikipedia.org/wiki/Вектор вектор-столбец], содержащий неизвестные коэффициенты <tex>c_i</tex>, <tex>\mathbf{y}</tex> - вектор-столбец, составленный из табличных значений функции <tex>y_i</tex>, а матрица <tex>\mathbf{A}</tex> имеет вид:
<tex> \mathbf{A} =
<tex> \mathbf{A} =
Строка 38: Строка 39:
</tex>
</tex>
-
Система линейных алгебраических уравнений (*) относительно неизвестных <tex>c_i</tex> иметь единственное решение, если определитель матрицы <tex>\mathbf{A}</tex> отличен от нуля.
+
Система линейных алгебраических уравнений (*) относительно неизвестных <tex>c_i</tex> иметь единственное решение, если [http://ru.wikipedia.org/wiki/Определитель определитель] матрицы <tex>\mathbf{A}</tex> отличен от нуля.
Определитель матрицы <tex>\mathbf{A}</tex> называют определителем Вандермонда, его можно вычислить по следующей формуле:
Определитель матрицы <tex>\mathbf{A}</tex> называют определителем Вандермонда, его можно вычислить по следующей формуле:
Строка 81: Строка 82:
=== Погрешность интерполяции ===
=== Погрешность интерполяции ===
-
Предположим, что на отрезке интерполирования [''a,b''] функция ''f(x)'' ''n'' раз непрерывно-дифференцируема. Погрешность интерполяции складывается из погрешности самого метода и ошибок округления.
+
Предположим, что на отрезке интерполирования [''a,b''] функция ''f(x)'' ''n'' раз непрерывно-дифференцируема. Погрешность интерполяции складывается из погрешности самого метода и [[Ошибки вычислений | ошибок округления.]]
Ошибка приближения функции ''f(x)'' интерполяционным полиномом ''n''-ой степени ''P<sub>n</sub>(x)'' в точке ''x'' определяется разностью: ''R<sub>n</sub>(x) = f(x) - P<sub>n</sub>(x).''
Ошибка приближения функции ''f(x)'' интерполяционным полиномом ''n''-ой степени ''P<sub>n</sub>(x)'' в точке ''x'' определяется разностью: ''R<sub>n</sub>(x) = f(x) - P<sub>n</sub>(x).''
Строка 156: Строка 157:
[[Изображение:Report1-fig3.gif]]
[[Изображение:Report1-fig3.gif]]
 +
 +
== Рекомендации программисту ==
 +
 +
Программа написана на языке C++ с использованием библиотеки линейной алгебры UBlas, которая является частью собрания библиотек Boost.
 +
 +
=== Предварительные настойки ===
 +
 +
Чтобы воспользоваться программой, необходимо сделать следующее:
 +
 +
1. Определиться с функцией, которую вы собираетесь интерполировать
 +
2. Создать текстовый файл (например, vec.txt), в первой строчке которого через пробел размещены узлы интерполяции, а во второй – значения выбранной функции в этих узлах.
 +
 +
Например, функция y = sin(x):
 +
 +
<pre>
 +
0.74 2 -3.5
 +
0.6743 0.9093 0.351
 +
</pre>
 +
 +
3. В .cpp файле программы в функцию double f(double x) вместо строки return прописать возвращаемое исходной функцией значение. Например, для функции y = sin(x):
 +
<pre>
 +
return sin(x);
 +
</pre>
 +
 +
=== Использование программы ===
 +
 +
В программе реализованы следующие основные функции:
 +
* '''double f(double x)''', описание которой было дано выше
 +
* '''double errapprox(vector coef, double a, double b, double h)''' – возвращает ошибку аппроксимации полиномом исходной функции.
 +
 +
На вход подаются следующие параметры:
 +
** '''vector coef''' – вектор коэффициентов интерполяционного полинома, который получается в ходе решения СЛАУ
 +
** '''double a, double b''' – границы промежутка интерполяции [a, b]
 +
** '''double h''' – шаг, с которым «пробегаем» промежуток [a, b]
 +
 +
* '''double outpolyn(char** filename, vector coef)''' – сохраняет коэффициенты полинома '''coef''' в файл '''filename'''
 +
 +
После запуска программы на экране появляются коэффициенты интерполяционного полинома и ошибка аппроксимации.
 +
 +
= Вывод =
 +
 +
Был исследован и программно реализован метод интерполяции функции каноническим полиномом. В ходе исследований установлено, что ошибка интерполяции получается как из-за ошибок компьютерных вычислений, так и из-за ошибок метода.
 +
 +
Также замечено, что от выбора узлов интерполяции напрямую зависит качество интерполяции. Минимальная ошибка интерполяции достигается при выборе «чебышевских» узлов.

Версия 16:31, 19 октября 2008

Содержание

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

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

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

Известно, что любая непрерывная на отрезке [a,b] функция f(x) может быть хорошо приближена некоторым полиномом Pn(x). Справедлива следующая Теорема (Вейерштрасса): Для любого \eps>0 существует полином Pn(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 точки лежат на одной прямой, то через них пройдёт единственный полином первой степени (но это ничему не противоречит: просто коэффициент при старшей степени равен нулю).

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Сложность решения СЛАУ (системы линейных алгебраических уравнений), например, методом Гаусса для матрицы размера nxn: 2n3/3, т.е. O(n3).

Для нахождения полинома в заданной точке требуется n умножений и n сложений. В результате сложность метода: O(n3).

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

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

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

Погрешность Rn(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).

Если максимальное значение производной fn+1(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).

При реализации данного метода на ЭВМ ошибкой интерполяции En(x) будем считать максимальное уклонение полинома от исходной функции на выбранном промежутке: E_n(x) = \max_{x \in [a, b]} \left| f(x) - P_n(x) \right| .

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

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


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


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


Известно (см. ссылки литературы), что если узлы интерполяции x0, x1,...,xn являются корнями полинома Чебышева степени n+1, то величина \max_{x \in [a,b]} \left| \omega_{n+1}(x) \right| принимает наименьшее возможное значение по сравнению с любым другим выбором набора узлов интерполяции.


Очевидно, что в случае k = 1 функция T1(x), действительно, является полиномом первой степени, поскольку T1(x) = cos(arccos x) = x.


В случае k = 2 T2(x) тоже полином второй степени. Это нетрудно проверить. Воспользуемся известным тригонометрическим тождеством: cos2θ = 2cos²θ - 1, положив θ = arccos x.

Тогда получим следующее соотношение: T2(x) = 2x² - 1.


С помощью тригонометрического тождества cos(k + 1)θ = 2cosθcoskθ - cos(k - 1) легко показать, что для полиномов Чебышева справедливо реккурентное соотношение:

Tk+1(x) = 2xTk(x) - Tk-1(x)


При помощи данного соотношения можно получить формулы для полиномов Чебышева любой степени.


Корни полинома Чебышева легко находятя из уравнения: Tk(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}}.


Известно, что для любого полинома pk(x) степени k с коэффициентом, равным единице при старшей производной верно неравенство \max_{x \in [-1,1]}p_k(x) \geq \frac1{2^{k-1}}, т.е. полиномы Чебышева являются полиномами, наименее уклоняющимися от нуля.

Вычислительный эксперимент

Для реализации поставленной задачи была написана программа на языке С++, которая по заданной функции приближает её каноническим полиномом. Разумеется, необходимо указать узлы, через которые полином пройдёт, и значения функции в этих узлах.

Далее строится СЛАУ, которая решается методом Гаусса. На выходе получаем коэффициенты для полинома и ошибку аппроксимации.

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

Пример 1: Интерполяция синуса

Попробуем интерполировать функцию y = sin(x) на отрезке [1, 8.5]. Выберем узлы интерполяции: {1.1, 2, 4.7, 7.5, 8.5}

Полученный в результате интерполяции полином отображён на рисунке (синим цветом показан график y = sin(x), красным – интерполяционного полинома)

Изображение:Report1-fig2.gif

Ошибка интерполяции в этом случае: 0.1534

Давайте посмотрим, что произойдёт, если выбрать равномерно стоящие узлы {2, 3.5 5, 6.5, 8} для той же функции на том же отрезке.

На отрезке [3, 6] приближение, бесспорно, стало лучше. Однако разброс на краях очень большой. Ошибка интерполяции: 2.3466

Изображение:Report1-fig3.gif

Рекомендации программисту

Программа написана на языке C++ с использованием библиотеки линейной алгебры UBlas, которая является частью собрания библиотек Boost.

Предварительные настойки

Чтобы воспользоваться программой, необходимо сделать следующее:

1. Определиться с функцией, которую вы собираетесь интерполировать 2. Создать текстовый файл (например, vec.txt), в первой строчке которого через пробел размещены узлы интерполяции, а во второй – значения выбранной функции в этих узлах.

Например, функция y = sin(x): 

0.74 2 -3.5
0.6743 0.9093 0.351

3. В .cpp файле программы в функцию double f(double x) вместо строки return прописать возвращаемое исходной функцией значение. Например, для функции y = sin(x): 

return sin(x);

Использование программы

В программе реализованы следующие основные функции:

  • double f(double x), описание которой было дано выше
  • double errapprox(vector coef, double a, double b, double h) – возвращает ошибку аппроксимации полиномом исходной функции.

На вход подаются следующие параметры:

    • vector coef – вектор коэффициентов интерполяционного полинома, который получается в ходе решения СЛАУ
    • double a, double b – границы промежутка интерполяции [a, b]
    • double h – шаг, с которым «пробегаем» промежуток [a, b]
  • double outpolyn(char** filename, vector coef) – сохраняет коэффициенты полинома coef в файл filename

После запуска программы на экране появляются коэффициенты интерполяционного полинома и ошибка аппроксимации.

Вывод

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

Также замечено, что от выбора узлов интерполяции напрямую зависит качество интерполяции. Минимальная ошибка интерполяции достигается при выборе «чебышевских» узлов.


Смотри также

Источник — «http://machinelearning.ru/wiki/index.php?title=%D0%98%D0%BD%D1%82%D0%B5%D1%80%D0%BF%D0%BE%D0%BB%D1%8F%D1%86%D0%B8%D1%8F_%D0%BA%D0%B0%D0%BD%D0%BE%D0%BD%D0%B8%D1%87%D0%B5%D1%81%D0%BA%D0%B8%D0%BC_%D0%BF%D0%BE%D0%BB%D0%B8%D0%BD%D0%BE%D0%BC%D0%BE%D0%BC»
Личные инструменты