Полиномиальная интерполяция

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

Содержание

    1 Теория
           Барицентрическое представление полинома
           Замечания по выбору узлов
    2 Полиномиальная интерполяция
           Функции
           Примеры
    3 Полиномиальная аппроксимация
    4 Ссылки по теме

Теория

Барицентрическое представление полинома

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

Такое представление обладает рядом удобных свойств, благодаря чему оно используется в пакете ALGLIB для работы с рациональными функциями. Та же формула может быть использована и для представления полинома, являющегося частным случаем рациональной функции. Барицентрическое представление имеет три важных преимущества перед другими способами интерполяции. Во-первых, оно более устойчиво к погрешностям, чем представление в степенном базисе. Во-вторых, после построения модели значение полинома может быть вычислено за время O(N), что быстрее, чем с использованием алгоритма типа Невилля (с трудоемкостью O(N²)). В-третьих, в ряде важных случаев (равномерные и Чебышевские сетки) трудоемкость построения модели имеет сложность O(N).

Замечания по выбору узлов

Выбор узлов очень сильно влияет на точность интерполяции. При этом сетка с равноудаленными узлами, не смотря на свою простоту и удобство в использовании, не подходит для интерполяции по двум тесно связанным причинам.

Чтобы понять первую причину, рассмотрим график справа. На графике изображена функция

(черным цветом) и интерполяционный полином, построенный по 11 точкам на равномерной сетке (красным цветом). Можно видеть, что в центре полином очень хорошо интерполирует функцию, однако по мере приближения к краям точность интерполяции падает. Хотя полином проходит через все указанные ему точки, но между точками он очень сильно отклоняется от функции. Можно было бы ожидать, что с ростом числа точек ошибка будет падать. На самом деле наоборот - чем больше n, тем сильнее интерполяционный полином будет отклоняться от функции по мере приближения к границам отрезка. При неограниченном росте числа точек ошибка интерполяции на отрезке будет стремиться к бесконечности.

Доказано, что существует целый класс функций, которые не могут быть интерполированы полиномом на равномерной сетке. Это функции, имеющие полюса на комплексной плоскости в окрестностях отрезка интерполяции (в частности, приведенная выше функция имеет полюса в точках x = +i и x = -i). Следует понимать, что рост ошибки интерполяции с ростом числа точек - не недостаток алгоритма и не следствие естественных погрешностей при операциях с вещественными числами. Это фундаментальное свойство интерполяционного полинома - в точности проходя через все предписанные точки, он будет резко возрастать в интервале между ними.

Вторая причина, по которой не следует использовать сетку с равноудаленными узлами, тесно связана с первой. При интерполяции на сетке с равноудаленными узлами погрешности операций с вещественными числами могут накапливаться и привести к снижению качества интерполяции. Причиной является то, что даже если интерполируемая функция относится к классу "хороших функций", т.е. не имеет полюсов в окрестностях отрезка интерполяции, погрешности операций с вещественными числами обычно вносят небольшие искажения в её график. И эти искажения очень часто принимают облик "плохой функции", что приводит к катастрофическому росту погрешности с ростом числа точек.

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

Полиномиальная интерполяция

Функции

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

polynomialbuild - для построения интерполяционного полинома на произвольной сетке за время O(N²) (интерполяция при этом имеет трудоемкость всего лишь O(N)).
polynomialbuildeqdist - для построения интерполяционного полинома на равномерной сетке за время O(N)
polynomialbuildcheb1 - для построения интерполяционного полинома на Чебышевской сетке I рода за время O(N)
polynomialbuildcheb2 - для построения интерполяционного полинома на Чебышевской сетке II рода за время O(N)

Во всех случаях результатом работы является структура barycentricinterpolant. Полученная структура может использоваться со следующими подпрограммами модуля ratint:

barycentriccalc, вычисляющей значение интерполянта в заданной точке
barycentricdiff1 и barycentricdiff2, вычисляющими производные интерполянта в заданной точке
barycentriclintransx, осуществляющей линейное преобразование аргумента полинома x=a·t+b
barycentriclintransy, осуществляющей линейное преобразование полинома P₂(x)=a·P(x)+b

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

polynomialbar2pow - осуществляет преобразование полинома из барицентрической формы в степенной базис (используйте эту функцию, если вам нужны коэффициенты полинома)
polynomialpow2bar - осуществляет преобразование полинома из степенного базиса в барицентрическое представление
polynomialbar2cheb - осуществляет преобразование полинома из барицентрической формы в базис полиномов Чебышева
polynomialcheb2bar - осуществляет преобразование полинома из базиса полиномов Чебышева в барицентрическое представление

Наконец, следующие функции можно использовать для быстрого - O(N) - вычисления значения интерполяционного полинома на одной из специальных сеток без создания барицентрической модели:

polynomialcalceqdist - для вычисления значения интерполяционного полинома в заданной точке с использованием значений функции на равномерной сетке
polynomialcalccheb1 - для вычисления значения интерполяционного полинома в заданной точке с использованием значений функции на Чебышевской сетке I рода
polynomialcalccheb2 - для вычисления значения интерполяционного полинома в заданной точке с использованием значений функции на Чебышевской сетке II рода

Примеры

ALGLIB Reference Manual содержит следующие примеры применения полиномиальной интерполяции:

polint_d_calcdiff - демонстрирует интерполяцию и дифференцирование с использованием барицентрического представления полинома.
polint_d_conv - демонстрирует преобразование между барицентрическим представлением и степенным базисом.
polint_d_spec - демонстрирует полиномиальную интерполяцию на специальных сетках.

Полиномиальная аппроксимация

Пакет ALGLIB также может использоваться для полиномиальной аппроксимации без ограничений (функция polynomialfit) или с ограничениями на значение полинома или его производной (функция polynomialfitwc). Более подробно этот вопрос рассмотрен в главе ALGLIB User Guide, посвященной аппроксимации.