Главная       Скачать       Коммерческая поддержка       FAQ       Forum       О нас       Английская версия

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

Рациональная интерполяция, т.е. интерполяция рациональными функциями, заключается в представлении интерполируемой функции f(x) в виде отношения двух полиномов:

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

Содержание

    1 История: первые алгоритмы рациональной интерполяции
    2 История: первые барицентрические алгоритмы
    3 Интерполяция без полюсов: алгоритм Флоатера-Хорманна
    4 Настройка алгоритма Флоатера-Хорманна
    5 Барицентрические модели: конструирование
    6 Барицентрические модели: интерполяция
    7 Барицентрические модели: аппроксимация
    8 Барицентрические модели: вычисление/дифференцирование
    9 Барицентрические модели: прочие операции
    10 Manual entries
    11 Ссылки по теме

История: первые алгоритмы рациональной интерполяции

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

Одним из первых способов построения рационального интерполянта было решение системы уравнений относительно коэффициентов полиномов p(x) и q(x) (числителя и знаменателя соответственно). Однако если число точек слишком велико, то полученная система линейных уравнений будет плохо обусловленной и полученные коэффициенты будут вычислены с ошибкой. Вообще, ошибка при вычислении коэффициентов будет присутствовать всегда, даже если система хорошо обусловлена, и даже самой небольшой погрешности в коэффициентах достаточно, чтобы рациональный интерполянт не проходил через требуемые точки. По этой причине вышеприведенный алгоритм не применяется на практике.

Для решения этой проблемы Булирш и Штер (Bulirsch and Stoer) предложили обобщение алгоритма Невилля на случай рациональной интерполяции. Алгоритма Булирша-Штера строит рациональную функцию со степенями числителя и знаменателя, равными N/2. У этого алгоритма есть два недостатка, которые мешают его применению на практике. Во-первых, для некоторых наборов точек в принципе невозможно построить интерполянт того вида, который строится этим алгоритмом, причем алгоритм не в состоянии обнаружить подобные случаи и сообщить об ошибке. Вторым недостатком является отсутствие механизма для обнаружения полюсов. Спецификой рациональной интерполяции является то, что даже если интерполируемая функция не имеет полюсов, интерполирующая функция может иметь неизвестное количество полюсов в неизвестных нам точках. Алгоритм Булирша-Штера не в состоянии проверить наличие полюсов в области интерполяции, что также говорит не в его пользу.

На диаграмме выше приведен немного искусственный пример, демонстрирующий оба недостатка алгоритма Булирша-Штера. Пять точек на графике соответствуют пяти точкам: x = [-2, -1, 0, 1, 2], y = [1, 2, -1, 0, 1]. На основе этих точек построен рациональный интерполянт со степенями числителя и знаменателя, равными 2 (именно такой интерполянт был бы построен алгоритмом Булирша-Штера). Можно видеть, что он имеет полюс в точке x = -0.5 - хотя, глядя на пять исходных точек, нельзя сказать, что эта функция имеет полюса. Также можно видеть, что интерполирующая функция не проходит через одну из предписанных точек - просто потому что не существует такой функции со степенями числителя и знаменателя, равными 2, которая могла бы пройти через эти точки. Ни одна из этих проблем не может быть автоматически выявлена при интерполяции алгоритмом Булирша-Штера.

История: первые барицентрические алгоритмы

Длительное время алгоритм Булирша-Штера при всех его недостатках оставался единственным алгоритмом рациональной интерполяции, который был доступен для практического применения. Даже в настоящее время многие специалисты по численному анализу не знают о последних достижениях в этой области. Однако ещё в 1986 году Шнайдер (Schneider C.) и Вернер (Werner W.) опубликовали статью, в которой был изложен новый алгоритм рациональной интерполяции, использующий барицентрическое представление рационального интерполянта:

При выборе весов барицентрической формулы w в соответствии с алгоритмом Шнайдера и Вернера выражение, приведенное выше, становится рациональной функцией с требуемой степенью знаменателя M (и степенью числителя N-M). Этот алгоритм не только осуществляет построение рационального интерполянта в барицентрической форме, но и позволяет проверить его на наличие полюсов и недостижимых точек. Впоследствии этот алгоритм был усовершенствован Беррутом и др. (см. 'Recent developments in barycentric rational interpolation', Jean–Paul Berrut, Richard Baltensperger and Hans D. Mittelmann, 2005 [2]).

Алгоритм, предложенный Беррутом и др., решает практически все проблемы, препятствовавшие использованию алгоритма Булирша-Штера. К сожалению, мир не идеален - вместо решенных проблем появляются новые. При росте N и M ухудшается обусловленность задачи, что приводит к росту погрешности w и к недостоверности выводов о наличии/отсутствии полюсов и недостижимых точек. Таким образом, надежность и устойчивость алгоритма недостаточно хороши, чтобы рекомендовать его в качестве универсального решения проблемы рациональной интерполяции.

Интерполяция без полюсов: алгоритм Флоатера-Хорманна

Практически одновременно со статьей Беррута и др. в печати появилась статья Флоатера и Хорманна[3], в которой был описан алгоритм построения рациональной интерполирующей функции, не имеющей полюсов на действительной оси. Важнами чертами алгоритма Флоатера-Хорманна являются высокая скорость работы, устойчивость и надежность. По этим параметрам он сравним со сплайн-интерполяцией. Алгоритм детально описан в статье, которую можно скачать в PDF-формате с сайта одного из авторов (см. ссылку), здесь же приведено краткое описание.

Пусть даны точки x, ..., x и даны значения функции в них: f = f(x). Пусть d - порядок интерполяционной схемы Флоатера-Хорманна (0 ≤ d ≤ n) и пусть p(x) обозначает полином, интерполирующий функцию f в точках x, ..., xi+d . Тогда рациональный интерполянт Флоатера-Хорманна r(x) определяется, как

Приведенная выше формула для r(x) дает понимание смысла параметра d, но она не подходит для применения на практике из-за высокой вычислительной сложности. Более быстрой и удобной является барицентрическая формула:

Определенная таким образом, функция r(x) обладает рядом важных свойств. Во-первых, это рациональная функция со степенями числителя и знаменателя не больше, чем N. Во-вторых, она не имеет полюсов на действительной оси. В-третьих, если h = max(x-xi-1 ), то ошибка аппроксимации имеет порядок O(h d+1). Можно видеть, что алгоритм является для рациональной интерполяции тем же, чем кубический сплайн является для полиномиальной интерполяции.

Настройка алгоритма Флоатера-Хорманна

Параметр d задает порядок используемой интерполяционной схемы и, как следствие, её точность. Интуитивно кажется, что чем больше d, тем меньше должна быть погрешность интерполяции, однако это верно лишь до определенных пределов. Ниже приведены две таблицы, сравнивающие ошибку интерполяции при разных значениях N и d на равномерной сетке на интервале [-5, 5]. Для сравнения использовались функция f = sin(x) и функция Рунге f = 1/(1+x 2) (похожий анализ был проведен в статье Флоатера и Хорманна).

f=sin(x)
N     D=3       D=5       D=8       D=N       Optimal D
 10   1.28e-02  5.21e-03  1.62e-03  6.30e-03  1.62e-03 (D=8)
 20   1.23e-03  1.26e-04  1.96e-05  1.72e-09  1.72e-09 (D=19)
 40   8.47e-05  3.17e-06  3.03e-08  9.53e-07  2.45e-13 (D=15)
 80   5.42e-06  5.50e-08  3.46e-11  2.05e+03  3.24e-14 (D=12)
160   3.39e-07  8.85e-10  4.00e-14  1.48e+02  2.78e-15 (D=10)


f=1/(1+x^2) (Runge example)
N     D=3       D=5       D=8       D=N       Optimal D
 10   6.91e-02  2.42e-01  1.41e+00  1.92e+00  3.61e-02 (D=0)
 20   2.83e-03  9.94e-03  7.04e-02  5.98e+01  1.54e-03 (D=1)
 40   4.31e-06  1.82e-05  1.30e-04  1.05e+05  4.31e-06 (D=3)
 80   5.12e-08  7.98e-10  4.54e-10  1.31e+10  2.03e-10 (D=7)
160   2.98e-09  1.14e-11  8.25e-15  4.57e+01  1.67e-15 (D=11)

Первая из рассматриваемых функций хорошо интерполируется обычными полиномами. Рациональная интерполяция также хорошо показывает себя при умеренных d. С ростом d с какого-то момента погрешность начинает расти, так что оптимальное d обычно находится где-то в диапазоне от 8 до 20. В статье Флоатера и Хорманна показано, что при d = 3 интерполяционная схема не уступает по точности кубическим сплайнам.

Функция Рунге является классическим примером функции, не поддающейся интерполяции полиномами на равномерной сетке. Рациональная интерполяция по методу Флоатера-Хорманна успешно справляется с задачей. Сравнивая зависимость ошибки от N и d с ошибкой при интерполяции f = sin(x), мы можем отметить, что оптимальное d в данном случае принимает меньшие значения.

Интересным является поведение ошибки при больших значениях d, в частности - при d = N. Можно видеть, что ошибка увеличивается, причем в некоторых случаях она может достигать очень больших значений (вплоть до 10 10). Анализ, проведенный с использованием арифметики высокой точности, показал, что в арифметике бесконечной точности с ростом d погрешность падает. Однако в арифметике конечной точности при d, близком к N, погрешности операций с вещественными числами практически полностью разрушают точность интерполяционной схемы. При 3 ≤ d ≤ 8 интерполяционная функция вычисляется с высокой точностью, близкой к машинной.

Подводя итог сказанному выше, можно сформулировать следующие правила выбора d. Если требуется высокая точность интерполяции и имеется такая возможность, то следует подобрать оптимальное d, обеспечивающее минимальную ошибку. Если подбор оптимального d невозможен, то следует выбирать d из диапазона от 3 до 8. Значение d = 3 одинаково хорошо в обоих рассмотренных случаях, так что оно может использоваться в качестве "значения по умолчанию". Для сравнения, d = 8 намного лучше, чем d = 3, если речь идет об интерполяции "хороших" функций (синус/косинус и т.п.), но плохо подходит для функций, подобных функции Рунге. Поэтому его не следует выбирать, если мы не знаем, к какому классу относится интерполируемая функция.

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

Если пользователь заранее подготовил точки x/y и весовые коэффициенты w, можно построить барицентрическую модель, используя функцию barycentricbuildxyw. Результатом является структура barycentricinterpolant, которая может использоваться для вычисления значения рациональной функции в заданной точке или для других операций.

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

Для решения задачи интерполяции следует использовать функцию barycentricbuildfloaterhormann. Эта функция строит барицентрическую модель, используя алгоритм Флоатера-Хорманна. Входными данными в этом случае являются набор точек и параметр d, результатом - опять структура barycentricinterpolant.

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

Пакет ALGLIB поддерживает рациональную аппроксимацию с использованием базиса Флоатера-Хорманна. Аппроксимация осуществляется двумя подпрограммами. Подпрограмма barycentricfitfloaterhormann решает задачу аппроксимации без ограничений и индивидуальных весов. Подпрограмма barycentricfitfloaterhormannwc решает более сложную задачу: аппроксимацию с индивидуальными весовыми коэффициентами и ограничениями на значения построенной модели. Поддерживается произвольное количество ограничений вида f(xc)=yc или df(xc)/dx=yc.

Замечание #1
Следует учитывать, что избыточные ограничения могут оказаться несовместимы. Более подробно этот вопрос рассмотрен в комментариях к подпрограмме barycentricfitfloaterhormannwc.

Обе подпрограммы используют базис Флоатера-Хорманна, построенный на равномерной сетке с M узлами. Функции перебирают разные значения D в интервале [0,9] в поисках оптимального, т.е. дающего минимальную взвешенную ошибку. Поскольку в такой формулировке задача линейна, то для решения используется линейный МНК-солвер и трудоемкость решения равна O(N*M^2), где N - число точек, M - размерность базиса.

Важно понимать, что для аппроксимации используется множество рациональных функций, гарантированно не имеющих полюсов - т.е. подмножество множества рациональных функций. С одной стороны, это делает невозможными резкие осцилляции интерполянта. С другой - снижает точность полученных моделей. Базис Флоатера-Хорманна позволяет аппроксимировать гладкую функцию с точностью 10 -3...10 -5, используя умеренное количество точек, однако для дальнейшего увеличения точности может потребоваться очень большое количество точек.

Барицентрические модели: вычисление/дифференцирование

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

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

Прочие операции включают в себя:

Ссылки по теме

  1. 'Recent developments in barycentric rational interpolation', Jean–Paul Berrut, Richard Baltensperger and Hans D. Mittelmann
  2. 'Barycentric rational interpolation with no poles and high rates of approximation', Michael S. Floater. and Kai Hormann

Manual entries

C++ ratint subpackage   
C# ratint subpackage   

This article is intended for personal use only.

Скачать ALGLIB

C#

Исходный код на C#

Downloads page

 

C++

Исходный код на C++

Downloads page

 

C++, арифметика высокой точности

Исходный код на C++, использующий библиотеки MPFR/GMP.

Исходный код GMP доступен на сайте gmplib.org. Исходный код MPFR доступен на сайте www.mpfr.org.

Downloads page

 

FreePascal

Исходный код на Free Pascal.

Downloads page

 

Delphi

Исходный код на Delphi.

Downloads page

 

VB.NET

Исходный код на VB.NET.

Downloads page

 

VBA

Исходный код на VBA.

Downloads page

 

Python

Исходный код на Python (CPython и IronPython).

Downloads page

 

 

ALGLIB® - numerical analysis library, 1999-2017.
ALGLIB is registered trademark of the ALGLIB Project.