Интерполяция каноническим полиномом
Материал из MachineLearning.
(→Постановка задачи) |
м (викификация, категория) |
||
Строка 1: | Строка 1: | ||
{{TOCright}} | {{TOCright}} | ||
- | Интерполя́ция — способ нахождения промежуточных значений величины по имеющемуся дискретному набору известных значений. | + | '''Интерполя́ция''' — способ нахождения промежуточных значений величины по имеющемуся дискретному набору известных значений. |
В настоящем исследовании будет изучена проблема интерполяции [http://ru.wikipedia.org/wiki/Функция_(математика) функции] одной переменной [http://ru.wikipedia.org/wiki/Многочлен полиномом,] будет рассмотрен вопрос точности приближения, и как, варьируя узлы, через которые пройдёт полином, достигнуть максимальной точности интерполяции. | В настоящем исследовании будет изучена проблема интерполяции [http://ru.wikipedia.org/wiki/Функция_(математика) функции] одной переменной [http://ru.wikipedia.org/wiki/Многочлен полиномом,] будет рассмотрен вопрос точности приближения, и как, варьируя узлы, через которые пройдёт полином, достигнуть максимальной точности интерполяции. | ||
Строка 55: | Строка 55: | ||
[[Изображение:Report1-fig1.gif|Зависимость числа обусловленности матрицы от количества узлов интерполяции]] | [[Изображение:Report1-fig1.gif|Зависимость числа обусловленности матрицы от количества узлов интерполяции]] | ||
- | Из-за плохой обусловленности матрицы рекомендуется применять другие методы интерполяции (например, метод Лагранжа). При этом важно понимать, что при теоретическом | + | Из-за плохой обусловленности матрицы рекомендуется применять другие методы интерполяции (например, метод Лагранжа). При этом важно понимать, что при теоретическом применении различных методов они приводят к одинаковому результату, т.е. мы получим один и тот же полином. |
Однако при практической реализации мы получим полиномы различной точности аппроксимации из-за погрешности вычислений аппаратуры. | Однако при практической реализации мы получим полиномы различной точности аппроксимации из-за погрешности вычислений аппаратуры. | ||
Строка 211: | Строка 211: | ||
* [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 Интерполяционный многочлен Лагранжа] | ||
* [[Метод наименьших квадратов]] | * [[Метод наименьших квадратов]] | ||
+ | |||
{{Stub}} | {{Stub}} | ||
+ | |||
+ | [[Категория:Учебные задачи]] |
Версия 18:49, 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). Справедлива следующая Теорема (Вейерштрасса): Для любого >0 существует полином Pn(x) степени , такой, что
В качестве аппроксимирующей функции выберем полином степени n в каноническом виде:
Коэффициенты полинома определим из условий Лагранжа , , что с учётом предыдущего выражения даёт систему линейных алгебраических уравнений с n+1 неизвестными:
Обозначим систему таких уравнений символом (*) и перепишем её следующим образом:
или в матричной форме: где - вектор-столбец, содержащий неизвестные коэффициенты , - вектор-столбец, составленный из табличных значений функции , а матрица имеет вид:
Система линейных алгебраических уравнений (*) относительно неизвестных иметь единственное решение, если определитель матрицы отличен от нуля.
Определитель матрицы называют определителем Вандермонда, его можно вычислить по следующей формуле:
Число узлов интерполяционного полинома должно быть на единицу больше его степени. Это понятно из интуитивных соображений: через 2 точки можно провести единственную прямую, через 3 - единственную параболу и т.д. Но полином может получиться и меньшей степени. Т.е. если 3 точки лежат на одной прямой, то через них пройдёт единственный полином первой степени (но это ничему не противоречит: просто коэффициент при старшей степени равен нулю).
При достаточной простоте реализации метода он имеет существенный недостаток: число обусловленности матрицы быстро растёт с увеличением числа узлов интерполяции, что можно показать на следующем графике
Из-за плохой обусловленности матрицы рекомендуется применять другие методы интерполяции (например, метод Лагранжа). При этом важно понимать, что при теоретическом применении различных методов они приводят к одинаковому результату, т.е. мы получим один и тот же полином.
Однако при практической реализации мы получим полиномы различной точности аппроксимации из-за погрешности вычислений аппаратуры.
Способ вычисления полинома в точке
Чтобы изобразить графически аппроксимирующий полином, необходимо вычислить его значение в ряде точек. Это можно сделать следующими способами.
Первый способ. Можно посчитать значение a1x и сложить с a0. Далее найти a2x2, сложить с полученным результатом, и так далее. Таким образом, на j-ом шаге вычисляется значение ajxj и складывается с уже вычисленной суммой .
Вычисление значения ajxj требует j операций умножения. Т.е. для подсчёта многочлена в заданной точке требуется (1 + 2 + ... + n) = n(n+1)/2 операций умножения и n операций сложения. Всего операций в данном случае: Op1 = n(n+1)/2 + n.
Второй способ. Полином можно также легко вычислить с помощью так называемой схемы Горнера:
Для вычисления значения во внутренних скобках 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) определяется следующим соотношением:
Здесь - производная (n+1)-го порядка функции f(x) в некоторой точке а функция определяется как
Если максимальное значение производной fn+1(x) равно то для погрешности интерполяции следует оценка:
При реализации данного метода на ЭВМ ошибкой интерполяции En(x) будем считать максимальное уклонение полинома от исходной функции на выбранном промежутке:
Выбор узлов интерполяции
Ясно, что от выбора узлов интерполируемой функции напрямую зависит, насколько точно многочлен будет являться её приближением.
Введём следующее определение: полиномом Чебышева называется многочлен вида
Известно (см. ссылки литературы), что если узлы интерполяции x0, x1,...,xn являются корнями полинома Чебышева степени n+1, то величина принимает наименьшее возможное значение по сравнению с любым другим выбором набора узлов интерполяции.
Очевидно, что в случае 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]: которые и следует выбирать в качестве узлов интерполирования.
Нетрудно видеть, что корни на [-1,1] расположены симметрично и неравномерно - чем ближе к краям отрезка, тем корни расположены плотнее. Максимальное значение модуля полинома Чебышева равно 1 и достигается в точках
Если положить то для того, чтобы коэффициент при старшей степени полинома ωk(x) был равен 1,
Известно, что для любого полинома pk(x) степени k с коэффициентом, равным единице при старшей производной верно неравенство т.е. полиномы Чебышева являются полиномами, наименее уклоняющимися от нуля.
Вычислительный эксперимент
Для реализации поставленной задачи была написана программа на языке С++, которая по заданной функции приближает её каноническим полиномом. Разумеется, необходимо указать узлы, через которые полином пройдёт, и значения функции в этих узлах.
Далее строится СЛАУ, которая решается методом Гаусса. На выходе получаем коэффициенты для полинома и ошибку аппроксимации.
Как было показано выше, и в чём мы убедимся в дальнейшем, от выбора узлов зависит точность, с которой полином будет приближать функцию.
Пример 1: Интерполяция синуса
Попробуем интерполировать функцию y = sin(x) на отрезке [1, 8.5]. Выберем узлы интерполяции: {1.1, 2, 4.7, 7.5, 8.5}
Полученный в результате интерполяции полином отображён на рисунке (синим цветом показан график y = sin(x), красным – интерполяционного полинома)
Ошибка интерполяции в этом случае: 0.1534
Давайте посмотрим, что произойдёт, если выбрать равномерно стоящие узлы {2, 3.5 5, 6.5, 8} для той же функции на том же отрезке.
На отрезке [3, 6] приближение, бесспорно, стало лучше. Однако разброс на краях очень большой. Ошибка интерполяции: 2.3466
Рекомендации программисту
Программа написана на языке 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
После запуска программы на экране появляются коэффициенты интерполяционного полинома и ошибка аппроксимации.
Вывод
Был исследован и программно реализован метод интерполяции функции каноническим полиномом. В ходе исследований установлено, что ошибка интерполяции получается как из-за ошибок компьютерных вычислений, так и из-за ошибок метода.
Также замечено, что от выбора узлов интерполяции напрямую зависит качество интерполяции. Минимальная ошибка интерполяции достигается при выборе «чебышевских» узлов.