Интерполяция функций

Задача интерполяции
Интерполяция – это способ нахождения промежуточных значений величины по имеющемуся дискретному набору известных значений.
Основное условие интерполяции: равенство функции и интерполяционного полинома в узлах интерполяции

Пример из области автоматизации проектирования
Математическая модель течения
Численная реализация
Геометрия объекта
Расчетная сетка
Визуализация результатов
Иллюстрации из: Плыкин М. FSI-технологии ANSYS // САПР и графика. – 2006. – Т. 7.

Методы построения интерполяционного полинома

Интерполяция алгебраическими полиномами
Определитель матрицы – детерминант Вандермонда. В случае различия всех узлов сетки он отличен от нуля, и, значит, существует единственное решение системы – набор коэффициентов.
Сеточная функция
Утверждение 2.1. Если заданы N + 1 узлов x0, …, xN среди которых нет совпадающих, и значения функции в этих узлах f (x0), ..., f (xN), то существует один и только один многочлен степени не выше N, принимающий в узлах xi заданные значения f (xi).

Метод Лагранжа

Примеры построенных методом Лагранжа полиномов

Метод Ньютона
Разделенная разность первого порядка
Разделенная разность второго порядка
Интерполяционный полином в форме Ньютона
Разделенная разность k-го порядка

Метод Ньютона
Разделенная разность k-го порядка
Интерполяционный полином в форме Ньютона

Примеры построенных методом Ньютона полиномов
Добавка к P1(x)

Лагранж vs Ньютон
Удобно применять, когда интерполируется одна и та же функция, но число узлов интерполяции постепенно увеличивается.
Удобно применять, когда узлы интерполяции фиксированы и интерполируется не одна, а несколько функций.
Рекомендуется использовать при программной реализации
Используют при доказательствах теорем
LN(x) и PN(x) – различные формы записи одного и того же многочлена

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

Погрешность интерполяции
Опр: Разница между функцией и интерполяционным полиномом N-ой степени в точке x называется остаточным членом интерполяции: RN(x) = f (x) – LN(x)
Утверждение 2.2. Пусть на отрезке [a,b] функция u(x) (N+1) раз непрерывно дифференцируема. Тогда:
Доказательство. Если x = xj, то утверждение верно. Иначе введем в рассмотрение функцию:
Функция g(t), имеет N + 2 нуля на [a, b]:

Погрешность интерполяции на равномерной сетке
Утверждение 2.3. Для случая равномерной сетки на отрезке [a,b]
для любого x на отрезке [a, b]
Доказательство. Положим
Тогда

Погрешность в задаче экстраполяции
Экстраполяция – аппроксимация функции вне отрезка, на котором заданы узлы интерполяции.
Экстраполяция функции менее надежна, чем интерполяция, и ее точность резко падает по мере удаления от носителя информации.
Оценка ухудшается как за счет появления множителей, пропорциональных N, так и за счет увеличения оценки производной.

N = 4
N = 5
N = 6

Увеличение числа узлов интерполяции
Сетка на [a, b]:
Рассмотрим последовательность сеток с возрастающим числом узлов:
Пусть f(x) определена и непрерывна на [a,b]. Построим последовательность интерполяционных многочленов для функции f(x) по ее значениям в узлах сетки ΩN: LN [f (x)].
Поточечная сходимость в точке x*∈ [a, b]:
Равномерная сходимость на [a, b]:

Пример С.Н. Берштейна u(x) = |x|, равномерная сетка
Нет сходимости ни в одной точке, кроме –1, 0, 1.
Колебательное свойство интерполяционных полиномов

U(x) = | x |, равномерная сетка vs оптимальная сетка
Возможное решение проблемы – сетка из нулей полинома Чебышева

Пример Рунге, u(x) = 1/( 1 + x2 ), равномерная сетка

Влияние неустранимой погрешности
Знаем как оценить
Необходимо оценить
Константа Лебега
Пример. Константа Лебега для случая линейной интерполяции

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

Недостатки глобальной интерполяции
При вычислении многочлена высокой степени могут накапливаться ошибки округления (например как при суммировании ряда Тейлора);
Интерполяционный многочлен может плохо приближать исходную функцию
(примеры Берштейна и Рунге на равномерной сетке);
Задача интерполяции может быть плохо обусловлена
(интерполяционный многочлен чувствителен к возмущениям значений в узлах сетки);
Частично указанные проблемы можно решить введением Чебышевской сетки, но не всегда такое возможно.
Ранее функцию U(x) мы искали в виде полинома от x
Рассмотрим теперь вариант, когда на каждом отрезке [xi-1, xi] функция является некоторым многочленом si(x), причем для каждого отрезка эта функция своя.
В такой постановке задача имеет множество решений. Единственность решения можно обеспечить потребовав от функции U(x) некоторой гладкости в местах стыков функций si(x), то есть в узлах интерполяции

Примеры сплайнов
На каждом отрезке функция приближается линейной. Дополнительных условия не требуется, условия гладкости на U(x) в данном случае не налагаются.
Гладкая кусочно-кубическая интерполяция
Кусочно-линейная интерполяция
На каждом отрезке функция приближается кубическим многочленом. Дополнительно требуется непрерывность первой и второй производных U(x) на всем отрезке [x0, xN].

Построение сплайна
Характеристики сплайна
Найдем выражения для функций si(x), составляющих гладкий кубический сплайн.
Гладкостью сплайна называется количество непрерывных производных, которые U(x) имеет на всем отрезке [x0, xN].
Дефектом сплайна называется разность между степенью и гладкостью сплайна.
Например, кусочно-линейный сплайн имеет степень 1, гладкость 0 и дефект 1.
Гладкий кусочно-кубический сплайн имеет степень 3, гладкость 2 и дефект 1.
Построение сплайна
Степенью сплайна называется максимальная из степеней многочленов si(x).
Поскольку сплайн имеет степень 3, все функции si(x) являются многочленами степени 3. Запишем их в виде:
Такая форма записи соответствует ряду Тейлора для si(x) в окрестности точки xi. Поскольку si(x) – кубический многочлен, его ряд Тейлора обрывается после кубического слагаемого. Из аналогии с рядом Тейлора заключаем, что
Хотя в этом можно убедиться и обычной подстановкой.

Построение сплайна
Условия непрерывности
Выпишем условия непрерывности первой и второй производной U(x) в точках xi-1:
Условия гладкости
Выразим условия непрерывности и гладкости сплайна в терминах коэффициентов ai, bi, ci, di. Для удобства введем обозначение для длины i-го отрезка hi = xi – xi-1. Запишем условие непрерывности U(x) в точке xi-1:

Построение сплайна
Основное условия интерполяции
Выпишем условия интерполирования, то есть U(xi) = u(xi):
Кроме этого, есть еще условия в точке x0,
Мы не требуем дополнительно si+1(xi) = U(xi), поскольку эти условия автоматически удовлетворяются при выполнении условия непрерывности.

Система для нахождения сплайна
В этой системе 3(N – 1) + N + 1 = 4N -2 уравнения и 4N неизвестных. Обычно, 2 недостающих условия задают в концах отрезка x0 и xN. В этом случае они называются краевыми условиями.
Объединим полученные ранее уравнения в единую систему
Краевые условия
«Естественный сплайн»
Понижение степени сплайна на краях до второй
Периодический сплайн
Рассмотрим наиболее используемый вариант

Линейная система
После добавления двух краевых условия количество уравнений совпало с количеством неизвестных. Можно было бы на этом остановиться, ведь формально, задача сведена к хорошо изученной.
Тем не менее, можно значительно упростить эту систему линейных уравнений, сведя ее к системе линейных уравнений специального трехдиагонального вида.
Начнем с исключения из системы неизвестных ai:
Удобно воспользоваться обозначениями разделенных разностей Ньютона

Упрощение системы
Подставим вместо ai значения u(xi):
Из уравнения (3) и (7) выразим dihi:
Исключим di из уравнений
и приведем подобные при сi

Упрощение системы
После приведения подобных
Выразим bi
И подставим в уравнение (2’’)
Заметим, что выражение для b1 формально совпадает с выражением для bi при i = 1, если доопределить с0 = 0.
Подставляя это в выражение (2’’) и упрощая, получаем

Трехдиагональная система
В результате серии упрощений у нас получилась система, относительно значений с1, … сN-1, причем, структура уравнений довольно специфическая. В i-е уравнение системы входят только три неизвестные.

Свойства сплайна
Оказывается, что если u(x) непрерывна, то последовательность кубических сплайнов UN(x) будет сходиться к u(x) равномерно, то есть
Построенный сплайн относится к глобальным. Если изменить значение u(xi) в какой-либо точке, это приведет к изменению всего сплайна U(x). Правда, амплитуда изменения быстро уменьшается при удалении от точки xi.

Литература
Петров И.Б., Лобанов А.И. Лекции по вычислительной математике: учеб. пособие. – М.: Интернет-Университет Информационных Технологий. Бином. Лаборатория знаний, 2006. – С. 133 – 141.
Самарский А.А., Гулин А.В. Численные методы. – М.: Наука, 1989. – С. 127 – 134.
Федоренко Р.П. Введение в вычислительную физику: учеб. пособие. – М.: Изд-во МФТИ, 1994. – С. 28 – 34.
Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. – М.: Бином. Лаборатория знаний, 2008. – С. 58 – 62.
Press W.H. et al. Numerical Recipes in C. – Cambridge University Press. – P. 120. – пример программной реализации построения интерполяционных полиномов.

Спасибо за внимание!

Оптимальный выбор узлов интерполяции. Многочлены Чебышева.

Минимизации погрешности интерполяции
Минимизируем за счет выбора
узлов интерполяции
Получили задачу на минимакс
(или задачу о построении полинома, наименее уклоняющемся от нуля на заданном отрезке):
Решение задачи – нормированный многочлен Чебышева степени N, а оптимальный выбор узлов интерполяции – нули многочлена Чебышева.

Многочлены Чебышева

Другая форма записи многочленов Чебышева
Функция cos( N arccos x ) удовлетворяет тому же разностному уравнению, что и TN(x)

Нули полиномов Чебышева
Отрезок [–1,1]
Отрезок [a,b]