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

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

(Различия между версиями)
Перейти к: навигация, поиск
м (викификация, оформление)
Строка 123: Строка 123:
== Вычислительный эксперимент ==
== Вычислительный эксперимент ==
-
Для реализации поставленной задачи была написана программа на языке С++, которая по заданной функции приближает её каноническим полиномом. Разумеется, необходимо указать узлы, через которые полином пройдёт, и значения функции в этих узлах.
+
Для реализации поставленной задачи была написана [[Media:Deryabin-prak1.zip| программа]] на языке С++, которая по заданной функции приближает её каноническим полиномом. Разумеется, необходимо указать узлы, через которые полином пройдёт, и значения функции в этих узлах.
Далее строится СЛАУ, которая решается методом Гаусса. На выходе получаем коэффициенты для полинома и ошибку аппроксимации.
Далее строится СЛАУ, которая решается методом Гаусса. На выходе получаем коэффициенты для полинома и ошибку аппроксимации.
Строка 145: Строка 145:
== Рекомендации программисту ==
== Рекомендации программисту ==
-
Программа написана на языке C++ с использованием библиотеки линейной алгебры UBlas, которая является частью собрания библиотек Boost.
+
Программа написана на языке C++ с использованием библиотеки линейной алгебры UBlas, которая является частью собрания библиотек [http://boost.org/ Boost.] Скачать исходный текст программы можно [[Media:Deryabin-prak1.zip| здесь [2.55Кб].]]
=== Предварительные настойки ===
=== Предварительные настойки ===
Чтобы воспользоваться программой, необходимо сделать следующее:
Чтобы воспользоваться программой, необходимо сделать следующее:
1. Определиться с функцией, которую вы собираетесь интерполировать
1. Определиться с функцией, которую вы собираетесь интерполировать
-
2. Создать текстовый файл (например, vec.txt), в первой строчке которого через пробел размещены узлы интерполяции, а во второй – значения выбранной функции в этих узлах.
+
2. Создать текстовый файл (например, '''vec.txt'''), в первой строчке которого через пробел размещены узлы интерполяции, а во второй – значения выбранной функции в этих узлах.
Например, функция y = sin(x):
Например, функция y = sin(x):
Строка 163: Строка 163:
return sin(x);
return sin(x);
</pre>
</pre>
 +
4. В функции '''int main()''' исходного кода присвоить переменной '''char* flname''' путь к входному файлу с данными. В нашем случае '''char* flname = "vec.txt";'''
=== Использование программы ===
=== Использование программы ===
В программе реализованы следующие основные функции:
В программе реализованы следующие основные функции:
* '''double f(double x)''', описание которой было дано выше
* '''double f(double x)''', описание которой было дано выше
-
* '''double errapprox(vector coef, double a, double b, double h)''' – возвращает ошибку аппроксимации полиномом исходной функции.
+
* '''int load(char *filename, vector<double> &x, vector<double> &y)''' - загрузка узлов интерполяции в переменную '''x''' и значения функции в этих узлах в переменную '''y''' текстового файла '''filename'''. В случае удачной загрузки данных из файла функция возвращает '''0'''.
 +
* '''void matrix2diag(matrix<double> &A, vector<double> &y)''' - приводит матрицу '''A''' к треугольному виду. '''y''' - столбец правой части (также изменяется вместе с матрицей '''A''').
 +
* '''void SolveSystem(matrix<double> &A, vector<double> &y, vector<double> &coef)''' - решить СЛАУ ('''A''' - треугольная матрица, '''y''' - столбец правой части, '''coef''' - в этот вектор заносится решение СЛАУ)
 +
* '''double errapprox(vector<double> coef, double a, double b, double h)''' – возвращает ошибку аппроксимации полиномом исходной функции.
 +
 
 +
На вход функции подаются следующие параметры:
-
На вход подаются следующие параметры:
 
** '''vector coef''' – вектор коэффициентов интерполяционного полинома, который получается в ходе решения СЛАУ
** '''vector coef''' – вектор коэффициентов интерполяционного полинома, который получается в ходе решения СЛАУ
** '''double a, double b''' – границы промежутка интерполяции [a, b]
** '''double a, double b''' – границы промежутка интерполяции [a, b]
** '''double h''' – шаг, с которым «пробегаем» промежуток [a, b]
** '''double h''' – шаг, с которым «пробегаем» промежуток [a, b]
-
* '''double outpolyn(char** filename, vector coef)''' – сохраняет коэффициенты полинома '''coef''' в файл '''filename'''
+
* '''int outpolyn(char** filename, vector<double> coef)''' – сохраняет коэффициенты полинома '''coef''' в файл '''filename'''. В случае удачного сохранения функция возвращает '0'.
После запуска программы на экране появляются коэффициенты интерполяционного полинома и ошибка аппроксимации.
После запуска программы на экране появляются коэффициенты интерполяционного полинома и ошибка аппроксимации.
Строка 182: Строка 187:
Также замечено, что от выбора узлов интерполяции напрямую зависит качество интерполяции. Минимальная ошибка интерполяции достигается при выборе «чебышевских» узлов.
Также замечено, что от выбора узлов интерполяции напрямую зависит качество интерполяции. Минимальная ошибка интерполяции достигается при выборе «чебышевских» узлов.
 +
 +
== Прикреплённые файлы ==
 +
 +
[[Media:Deryabin-prak1.zip| Программная реализация на языке С++ (исходный текст) [2.55Кб]]]
 +
 +
== Литература ==
 +
 +
* ''Н.С. Бахвалов, Н.П. Жидков, Г.М. Кобельков.'' Численные методы. Изд-во "Лаборатория базовых знаний". Москва. 2003.
 +
* ''И.С. Березин, Н.П. Жидков.'' Методы вычислений. Изд. ФизМатЛит. Москва. 1962.
 +
* ''Дж. Форсайт, М.Мальком, К. Моулер.'' Машинные методы математических вычислений. Изд-во "Мир". Москва. 1980.
== Смотри также ==
== Смотри также ==
Строка 187: Строка 202:
* [http://ru.wikipedia.org/wiki/%D0%98%D0%BD%D1%82%D0%B5%D1%80%D0%BF%D0%BE%D0%BB%D1%8F%D1%86%D0%B8%D0%BE%D0%BD%D0%BD%D1%8B%D0%B9_%D0%BC%D0%BD%D0%BE%D0%B3%D0%BE%D1%87%D0%BB%D0%B5%D0%BD_%D0%9B%D0%B0%D0%B3%D1%80%D0%B0%D0%BD%D0%B6%D0%B0 Интерполяционный многочлен Лагранжа]
* [http://ru.wikipedia.org/wiki/%D0%98%D0%BD%D1%82%D0%B5%D1%80%D0%BF%D0%BE%D0%BB%D1%8F%D1%86%D0%B8%D0%BE%D0%BD%D0%BD%D1%8B%D0%B9_%D0%BC%D0%BD%D0%BE%D0%B3%D0%BE%D1%87%D0%BB%D0%B5%D0%BD_%D0%9B%D0%B0%D0%B3%D1%80%D0%B0%D0%BD%D0%B6%D0%B0 Интерполяционный многочлен Лагранжа]
* [[Метод наименьших квадратов]]
* [[Метод наименьших квадратов]]
 +
 +
== Ссылки ==
 +
* [http://boost.org/ Собрание библиотек Boost]
{{Stub}}
{{Stub}}
[[Категория:Учебные задачи]]
[[Категория:Учебные задачи]]

Версия 10:34, 20 октября 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. Скачать исходный текст программы можно здесь [2.55Кб].

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

Чтобы воспользоваться программой, необходимо сделать следующее: 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);

4. В функции int main() исходного кода присвоить переменной char* flname путь к входному файлу с данными. В нашем случае char* flname = "vec.txt";

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

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

  • double f(double x), описание которой было дано выше
  • int load(char *filename, vector<double> &x, vector<double> &y) - загрузка узлов интерполяции в переменную x и значения функции в этих узлах в переменную y текстового файла filename. В случае удачной загрузки данных из файла функция возвращает 0.
  • void matrix2diag(matrix<double> &A, vector<double> &y) - приводит матрицу A к треугольному виду. y - столбец правой части (также изменяется вместе с матрицей A).
  • void SolveSystem(matrix<double> &A, vector<double> &y, vector<double> &coef) - решить СЛАУ (A - треугольная матрица, y - столбец правой части, coef - в этот вектор заносится решение СЛАУ)
  • double errapprox(vector<double> coef, double a, double b, double h) – возвращает ошибку аппроксимации полиномом исходной функции.

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

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

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

Вывод

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

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

Прикреплённые файлы

Программная реализация на языке С++ (исходный текст) [2.55Кб]

Литература

  • Н.С. Бахвалов, Н.П. Жидков, Г.М. Кобельков. Численные методы. Изд-во "Лаборатория базовых знаний". Москва. 2003.
  • И.С. Березин, Н.П. Жидков. Методы вычислений. Изд. ФизМатЛит. Москва. 1962.
  • Дж. Форсайт, М.Мальком, К. Моулер. Машинные методы математических вычислений. Изд-во "Мир". Москва. 1980.

Смотри также

Ссылки

Источник — «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»
Личные инструменты