ЛАБОРАТОРНАЯ РАБОТА 2

ТАБЛИЧНОЕ ЗАДАНИЕ И ИНТЕРПОЛИРОВАНИЕ ФУНКЦИЙ

1. Введение

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

2. Теоретическая справка

2.1. Задача интерполяции

Задача интерполяции состоит в нахождении обобщенного многочлена

                                                                    (2.1)

где  — фиксированные функции, а значения коэффициентов определяются из условия равенства со значением приближаемой функции в узлах интерполяции

                                                          (2.2)

Набор точек  на интервале [ab],
 в которых заданы значения функции  называют сеткой. Множество точек  иногда также называют узлами сетки или узлами интерполяции.

Мы будем называть сетку равномерной, если

           

       

Если  то соответствующая интерполяция называется алгебраической, если  — тригонометрические функции, то говорят о тригонометрической интерполяции.

Если построенный многочлен (2.2) используется для восстановления функции на всем отрезке [a, b], то говорят о глобальной интерполяции. Если же для восстановления функции между каждыми двумя соседними узлами строится многочлен заданной невысокой степени, то говорят о кусочно-многочленной интерполяции.

Если значения функции f(x) заданы в узлах  на интервале [ab],  то говорят, что функция f(x) задана таблицей.

2.2. Алгебраическая интерполяция

Теорема 1. Пусть задан n + 1 узел   …,  среди которых нет совпадающих, и значения функции в этих узлах   …,  Тогда существует один и только один алгебраический многочлен  степени не выше n, принимающий в узлах  заданные значения

Интерполяционный многочлен можно записать (и соответственно вычислить) различными способами представляя его в виде разложения по степеням xформе Лагранжа и в форме Ньютона), или в виде разложения по ортогональным многочленам.

2.3. Непосредственное вычисление коэффициентов интерполяционного полинома

Полином степени n можно записать в виде

                                         (2.3)

где  …,  — неопределенные коэффициенты. Их можно определить из n+1-го условия:

........................................................................                                    (2.4)

Определитель системы (2.4) есть детерминант Вандермонда, известный из курса линейной алгебры. Его значение в случае, когда выполняются условия теоремы 1, отлично от нуля, что доказывает существование и единственность полинома. Эта линейная система во многих случаях является плохо обусловленной. Последнее связано с тем, что последовательные степени 1, x,  …,  «почти линейно зависимы» на интервале 0 < x < 1. Напомним, что обусловленность линейной системы Ay = b определяется числом

                                                                          (2.5)

которое определяет относительную погрешность решения системы в зависимости от относительной погрешности правой части b:

                                                                               (2.6)

Матрицу A будем называть сингулярной, если в рамках системы вычислений с плавающей точкой на данной машине выполняется равенство

2.4. Интерполяционный полином в форме Лагранжа. Интерполяционный полином в форме Ньютона

Введем вспомогательные многочлены

           (2.7)

Многочлен  заданный равенством

          (2.8)

есть интерполяционный многочлен в форме Лагранжа.

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

Определим разностные отношения (иногда употребляется термин «разделенные разности»). Пусть функция  в точках     и т. д. принимает значения    

Разностные отношения нулевого порядка  функции  в точке  определяют, как значение функции в этой точке  k = a, b, c, d … Разностные отношения первого порядка  функции  для произвольной пары точек  и  определим через разностные отношения нулевого порядка:

                                                          (2.9)

Разностное отношение  порядка n определим через разностное отношение порядка  положив:

            (2.10)

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

      (2.11)

Если  то соответствующую интерполяцию называют интерполяцией вперед; в случае  интерполяцию называют интерполяцией назад.

2.5. Формула для погрешности алгебраической интерполяции

Оценим погрешность   возникающую при приближенной замене  алгебраическим многочленом  В основе оценки лежит следующая общая теорема о формуле погрешности.

Теорема 2. Пусть  — функция, определенная на некотором отрезке  и имеющая производные до некоторого порядка s + 1 включительно.

Пусть   …,  — произвольный набор попарно различных точек из отрезка    …,  — значения функции  в этих точках;  — интерполяционный многочлен степени не выше s, построенный по этим значениям. Тогда погрешность интерполяции  представляется формулой:

                            (2.12)

где  — некоторая точка интервала

2.6. О сходимости интерполяционного процесса

На отрезке a ≤ х ≤ b будем рассматривать бесконечную последовательность узлов интерполяции

.......................................................................................................... (2.13)

..................................

и соответствующую последовательность интерполяционных многочленов  построенную для некоторой функции  принимающей конечное значение.

Теорема 3. (Фабера). Какова бы ни была последовательность узлов интерполяции, существует непрерывная функция f, для которой последовательность интерполяционных многочленов расходится.

Теорема 4. Для каждой функции f, непрерывной на конечном отрезке, существует такая последовательность узлов интерполяции, что соответствующий ей интерполяционный процесс равномерно сходится к f.

Теорема 5. Не существует последовательности узлов, для которой интерполяционный процесс был бы равномерно сходящимся для всякой непрерывной на отрезке функции.

Теорема 6. Если функция f имеет ограниченную производную на отрезке, то интерполяционный процесс, в котором за узлы принимаются корни многочленов Чебышёва, сходится равномерно к f.

Основные термины. Интерполирующую функцию иногда называют интерполянтом.

Гладкий кусочно-многочленный интерполянт называется сплайном.

2.7. Обусловленность задачи построения интерполяционного многочлена для функции, заданной таблицей

Интерполяционный многочлен

                                                                   (2.14)

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

                                                                   (2.15)

для      иногда называют фундаментальными полиномами.

Придадим значениям функции  возмущения  Интерполяционный многочлен  заменится многочленом

Так как в силу линейности (2.15) по f, то возмущение  которое претерпевает интерполяционный многочлен, можно оценить как:

                               (2.16)

Это возмущение при заданных узлах интерполяции и фиксированных базисных функциях  зависит только от df.

Введем в рассмотрение функцию  которая называется функцией Лебега.

За меру чувствительности интерполяционного многочлена к возмущениям задания функции в узлах df принимается наименьшее число  при котором для каждого df выполнено неравенство

                                   (2.17)

Очевидно, что

Числа   называют константами Лебега. Эти числа растут с ростом n. Их поведение при возрастании n существенно зависит от расположения узлов интерполяции на отрезке.

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

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

Если узлами интерполяции являются корни полинома Чебышёва или точки, где этот полином достигает экстремумов, то  т. е. с ростом n константы Лебега растут очень медленно. В этом случае вычислительная неустойчивость не является препятствием для использования интерполяционных многочленов высокой степени.

2.8. Классическая кусочно-многочленная интерполяция

Пусть функция  задана таблицей. Для восстановления функции между узлами можно воспользоваться функцией, которая между каждыми двумя соседними узлами является многочленом заданной невысокой степени, например, первой, второй, третьей и т. д.

Соответствующая интерполяция называется кусочно-линейной, кусочно-квадратичной и т. п.

2.9. Оценка неустранимой погрешности при приближении функции по ее значениям в узлах интерполяции. Выбор степени кусочно-многочленной интерполяции

Пусть функция  определена на отрезке  и пусть заданы ее значения в узлах равномерной сетки  
 По таблице   …,  в принципе, нельзя восстановить функцию  точно, потому что значения различных функций могут совпадать в точках  = 1, …, n, т. е. различные функции могут иметь одинаковую таблицу.

Если, о функции известно лишь то, что она непрерывна, то ее нельзя восстановить в точке  k = 0, 1, …, n, ни с какой гарантированной точностью.

Пусть о функции  известно, что она имеет производные порядка + 1, причем

                                                 (2.18)

Укажем две функции

                                       (2.19)

для которых таблицы  k = 0, 1, …, n, совпадают (обе таблицы содержат лишь нули) Эти функции уклоняются друг от друга на величину порядка

           (2.20)

Таким образом, зная лишь оценку s + 1 производной, в принципе нельзя восстановить табличную функцию с точностью, большей, чем  Данная погрешность неустранима.

2.10. Насыщаемость (гладкостью) кусочно-многочленной интерполяции

Пусть функция  определена на отрезке [a, b], и задана ее таблица  в равноотстоящих узлах   с шагом

Погрешность кусочно-многочленной интерполяции степени s (с помощью интерполяционных многочленов  на отрезке ) в случае, если на [a, b] существует и ограничена , имеет порядок

Если о функции  известно лишь, что она имеет ограниченную производную до некоторого порядка q,  то неустранимая погрешность при ее восстановлении по таблице есть  Можно показать, что при интерполяции с помощью  порядок  достигается.

Если  имеет ограниченную производную порядка + 1, s, то погрешность интерполяции с помощью  остается  т. е. порядок погрешности не реагирует на дополнительную, сверх + 1 производной, гладкость функции  Это свойство кусочно-многочленной интерполяции называют свойством насыщаемости (гладкостью).

2.11. Кусочно-многочленная гладкая интерполяция (сплайны). Локальные сплайны

Классическая кусочно-линейная, кусочно-квадратичная и, вообще, кусочно-многочленная интерполяция заданной степени приводят к интерполирующей функции (интерполянту), которая в узлах интерполяции, вообще говоря, не имеет производной даже первого порядка. Существует два типа гладких кусочно-много­членных интерполянтов — локальные и нелокальные сплайны. Они обладают заданным числом производных всюду, включая узлы интерполяции.

Рассмотрим локальные сплайны. Пусть заданы узлы интерполяции   и значения функции  в них. Зададим натуральное число s, фиксируем натуральное число j: ≤ s. Каждой точке  сопоставим интерполяционный многочлен  построенный по значениям   …,  в узлах   …,  Кусочно-многочленную интерполирующую функцию  имеющую непрерывные производные порядка s, определим равенствами

       k = 0, ±1, ±2, …                                                                        (2.21)

где  — многочлен степени не выше 2+ 1, определяемый равенствами

                                                                                                           (2.22)

Верны следующие теоремы.

Теорема 7. Существует один и только один многочлен сте­пени не выше 2s + 1,  удовлетворяющий (2.22).

Теорема 8. Пусть  — многочлен степени не выше s. Тогда интерполянт  совпадает с этим многочленом.

Теорема 9. Кусочно-многочленная интерполирующая функция  определяемая равенством (2.22), в узлах интерполяции  совпадает с заданным в них значением  = 0, 1, …

Кроме того,  имеет всюду в области своего определения непрерывную производную порядка s.

Многочлен  можно записать в виде

                                    (2.23)

обозначив через  поправку к классическому интерполяционному многочлену

Рассмотрим здесь наиболее интересный для приложений случай, когда s = 2, j = 1, а узлы интерполяции функции составляют равномерную сетку. Тогда

          (2.24)

Эта формула справедлива только при 0 < k < n. На отрезках
 и  поправка:

2.12. Нелокальная гладкая кусочно-многочленная
интерполяция

Пусть задана таблица. Поставим задачу найти на каждом отрезке  кубический многочлен  так, чтобы возникающая при этом на отрезке  кусочно-многочленная функция совпадала с заданной функцией в узлах и имела непрерывные производные до порядка s = 2. Общее число неизвестных — 4n. Число дополнительных условий равно 4n – 2. Недостающие условия можно задавать различными способами. Наиболее употребляемыми являются следующие два:

 — «свободный сплайн»;           (2.25)

                                      (2.26)

Здесь  — единственные кубические кривые, которые проходят соответственно через четыре первые и четыре последние из заданных точек.

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

Однако сплайны Шонберга теряют свойства локальности, присущие как классической кусочно-многочленной интерполяции, так и локальной гладкой интерполяции: коэффициенты многочлена, задающие интерполянт на каком-либо отрезке  зависят от значений функции во всех узлах сетки.

2.13. Тригонометрическая интерполяция

Задача (линейной) тригонометрической интерполяции состоит в нахождении тригонометрического многочлена вида

         (2.27)

Здесь k и n — натуральные числа,  — положительное число,  — отрезок интерполяции,  и  — числовые  коэффициенты.

Теорема 10. (Первый вариант задания узлов интерполяции). Пусть N = 2(n+1), n — натуральное число. При произвольном задании значений функции  периодической с периодом L, в узлах сетки

               = 0, 1, …, N  1

существует один и только один интерполяционный тригонометрический многочлен

         (2.28)

удовлетворяющий равенствам  m = 0, …, N – 1.

Коэффициенты этого многочлена задаются формулами

                     

              (2.29)

    

Теорема 11. (Bторой вариант задания узлов интерполяции). Пусть N = 2n. При произвольном задании значений функции  периодической с периодом L, в узлах сетки

               m = 0, ±1, …, ±(N  1)

существует один и только один интерполяционный тригонометрический многочлен

         (2.30)

удовлетворяющий равенствам  m = 0, ±1, ±2, …, ± (– 1).

Коэффициенты этого многочлена задаются формулами

                           

                      (2.31)

         

Теорема 12. (Произвольное расположение узлов интерполяции на отрезке периодичности). Пусть заданы значения  i = 1, …, N, периодической с периодом L функции  в N = 2n несовпадающих точках  i = 1, …, N, принадлежащих отрезку   

Тогда существует один и только один интерполяционный тригонометрический многочлен

      (2.32)

                                                                                                                    (2.33)

(В произведениях отсутствует сомножитель, соответствующий i = k, так что ).

Построенные тригонометрические многочлены обладают определенными преимуществами перед алгебраическим многочленом, построенным по значениям функции в узлах

Во-первых, при  погрешность тригонометрической интерполяции

равномерно стремится к нулю, если  имеет хотя бы вторую производную, причем скорость убывания погрешности автоматически учитывает гладкость  т. е. возрастает с ростом числа (r+1) производных:

       

                                                                                                                    (2.34)

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

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

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

Многочлены Чебышёва можно ввести по формуле
 при  . Функции  суть многочлены степени k = 0, 1, … При этом     … вычисляются по рекуррентной формуле

                                                  (2.35)

Нули  определяются из уравнения:

                                                       (2.36)

т. е.  

Точки экстремума  определяются как:

                        1 = 0, 1, …, k.

Таким образом,  на отрезке  имеет k вещественных нулей и  точку экстремума. На оси x эти точки получаются проекцией точек пересечения полуокружности с множеством лучей, имеющих между собой равные углы.

2.15. Алгебраический интерполяционный полином на сетке из нулей полинома Чебышёва

При выборе в качестве узлов интерполяции нулей полинома Чебышева интерполяционный полином можно записать в виде:

                                                              (2.37)

где   n + 1 — число узлов интерполяции.

2.16. Алгебраический интерполяционный полином на сетке из экстремумов полинома Чебышёва

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

                                                              (2.38)

где

Здесь n + 1 — число узлов интерполяции.

2.17. Чувствительность интерполяционного тригонометрического многочлена к погрешностям задания функции в узлах интерполяции

Чувствительность интерполяционного тригонометрического многочлена к погрешности задания значений  оценивается следующим образом. Пусть вместо  задана сеточная функция  Тогда возникающая погрешность

                                          (2.39)

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

                                                       (2.40)

Теорема 13. Константы Лебега тригонометрического интерполяционного многочлена удовлетворяют оценке

Интерполяционный полином на сетке из нулей или экстремумов полиномов Чебышёва наследует от тригонометрической интерполяции слабый рост константы Лебега при увеличении n. Для этого случая справедлива оценка

2.18. Библиографическая справка

Теория интерполяции — обширный раздел вычислительной математики. В данной работе изложение ведется по книге [1]. Другие аспекты тории изложены в [9, 10], в частности, подробно рассмотрено свойство насыщаемости алгоритмов. О построении глобального сплайна см. [11]. О локальных сплайнах, кроме [1], см. также [12].

В связи с развитием алгоритмов машинной графики бурно развивается теория сплайн-интерполяции. О применении сплайнов и кривых Безье см. также [5]. Различным приложениям сплайнов посвящены книги [13–15]. В [5, 14]дается понятие о В-сплайнах, нашедших широкое применение как в алгоритмах машинной графики, так и в развитии численных методов. О применении В-сплайнов в инженерных расчетах см. в [16]. В [15] описаны алгоритмы построения сглаживающих сплайнов, находящие применение, в частности, в обработке экспериментальных данных.

3. Задание

1. Войдите в меню «Параметры/Таблица» и выберите из предложенного списка функций в пункте «Функции» какую-нибудь гладкую (т. е. имеющую непрерывные производные) функцию. В меню «Сетка» выберите пункт равномерная, число узлов, равное двум. Алгебраический интерполяционный полином какой степени может быть построен по этим данным?

Выберите в подменю «Метод/Глобальная» строку «Интерполяционный полином в форме Лагранжа», степень которого 3, 4, 5, 6, 7, ¼, 50. Постройте интерполяционный полином. Теоретически оцените погрешность интерполяции и сравните ее с фактической (т. е. разностью между исходной функцией и ее интерполянтом по данной таблице) погрешностью, выбрав пункт меню «Окна/Ошибка». Сформируйте таблицу, выбрав функцию, имеющую:

а) только две непрерывных производных

б) одну непрерывную производную

в) не имеющую непрерывной производной

Как будет меняться фактическая погрешность восстановления функции с ростом степени полинома для равномерной сетки? Почему в пункте a) для нечетных степеней k интерполяционного полинома ошибка меньше, чем для  1 и k + 1? То же для сетки из нулей полинома Чебышева. Чем объяснить существенное различие в поведении погрешности?

2. Выберите целую (разлагающуюся в сходящийся степенной ряд для любого конечного x) функцию (например, ) и постройте на сетке из равноотстоящих узлов глобальный алгебраический интерполянт. Сравните различные способы вычисления интерполяционного полинома:

а)      находя его коэффициенты по базису  решая соответствующую линейную систему;

б, в) записывая интерполяционный многочлен в форме Лагранжа; в форме Ньютона.

Какой способ требует меньшего числа арифметических действий?

Выберите число узлов  Объясните проблемы, возникающие при численном решении линейной системы. Почему результаты использования формы Ньютона и Лагранжа в начале совпадают, а затем начинают различаться?

Проанализируйте накопление погрешности при вычислении полинома в форме Ньютона при перенумерации в порядке возрастания и в порядке убывания (интерполяция вперед и назад)?

Как изменятся результаты, если вычисления производить с двойной (мантисса 52 бита) и стандартной (мантисса 24 бита) точностью?

3. Проанализируйте влияние погрешности задания функции в узлах на величину фактической погрешности восстановления функции. В силу линейности интерполяции эти эффекты удобно изучать на функции у = 0. Как задать распределение погрешности для получения максимального отклонения между узлами интерполяции? Рассмотрите случай равноотстоящих узлов и сетку из нулей полинома Чебышева. Сравните с теоретической оценкой.

4. Сравните фактическую погрешность интерполяции функций  и  (функция Рунге) на отрезке  на равномерной сетке при числе узлов n = 11, 21, 31, 41, 51. В чем причина отсутствия сходимости интерполяционного процесса для функции Рунге? Задайте теперь отрезок интерполяции  и повторите вычисления. Попытайтесь объяснить полученный результат.

Указание. Обратите внимание на то, что функция  целая, а функция Рунге имеет полюсы при

5. Посмотрите, как ведет себя глобальный интерполянт вне отрезка, на котором расположены узлы интерполяции. Для этого в меню «Окно» установите соответствующие границы изменения аргумента x для отображения  и  Что лучше использовать для экстраполяции: глобальный интерполянт или кусочно-многочленный?

6. Сравните ошибку при алгебраической интерполяции и интерполяции сплайнами на равномерной сетке из 3, 4, …, 25 узлов на отрезке [–1, 1] следующих функций:

а)              б)            в) функции Рунге.

Объясните результаты.

7. Интерполирующим кубическим сплайном (Шонберга) приближается функция

на отрезке

Почему в окрестности разрыва возникают осцилляции? Как меняется их амплитуда в зависимости от изменения числа узлов сетки? Объясните эффект. Ответьте на те же вопросы для локального сплайна.

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

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

4. Контрольные вопросы

1. Докажите, что при выборе в качестве узлов интерполяции нулей полинома Чебышева, алгебраический интерполяционный полином можно записать в виде

где

       .

Здесь n + 1 — число узлов интерполяции.

2. Докажите, что при выборе в качестве узлов интерполяции экстремумов полинома Чебышева алгебраический интерполяционный полином можно записать в виде

где коэффициенты

Здесь n + 1 — число узлов интерполяции.

3. Выведите формулы для интерполяции табличной функции кубическим сплайном (Шонберга):

а) на равномерной сетке;

б) на неравномерной сетке.

Указание. Пусть  — значение второй производной в i-м узле,  — значение второй производной в i + 1-м узле. На отрезке   — линейная функция. На этом отрезке  — непрерывная функция, которую можно получить, дважды интегрируя  по x. При этом возникают две константы интегрирования. Они находятся из условий:   Для значения вторых производных  строится система алгебраических уравнений из условия