Сопоставление интерполяционных формул Лагранжа и Ньютона. Погрешность интерполяции
Интерполяционные многочлены Лагранжа и Ньютона, построенные для одних и тех же узлов интерполяции, тождественно равны между собой, хотя и имеют различную форму записи. Это вытекает из единственности интерполяционного многочлена заданной степени. Разница в построении алгоритмов может учитываться в связи с особенностью решаемой задачи.
Коэффициенты Лагранжа зависят от выбора узлов и точки , но не зависят от вида функции . Это удобно, когда по заданной системе узлов надо интерполировать несколько различных функций.
Выбор способа интерполяции определяется различными соображениями: точностью, временем вычислений, погрешностью округлений и т.д. В ряде случаев более выгодной может оказаться локальная интерполяция, а не построение многочлена высокой степени. Во многих случаях интерполяционный многочлен
Ньютона более удобен, чем интерполяционный многочлен Лагранжа. Особенность этого многочлена заключается в том, что при переходе от многочлена -ой степени к многочлену -й степени первые членов не меняются, а только добавляется новый член, который равен нулю при всех предыдущих значениях аргумента. Формула Лагранжа этого делать не позволяет, так как в ней добавление нового узла заставляет заново пересчитывать все коэффициенты .
В точках, отличных от узлов интерполирования, значения функций и не совпадают: . Эта разность – погрешность интерполяции – называется остаточным членом.
Если интерполируемая функция имеет непрерывные производные до порядка включительно, погрешность при замене функции многочленом , т.е. величина , удовлетворяет неравенству
Если отрезок конечен, то
Таким образом, если при растут не слишком быстро, то в этом случае абсолютная величина погрешности стремиться к нулю для каждого . В этом случае функцию можно приблизить сколь угодно точно полиномом Лагранжа сразу для всех , если степень многочлена достаточна велика. Как правило, увеличение числа узлов улучшает приближение, но существуют и отклонение от этого правила. В некоторых случаях точность может быть повышена за счет расположения узлов интерполяции.
Аналогично оценивается погрешность и для интерполяционной формулы Ньютона.
2.6. Интерполирование функции кубическими сплайнами
Пусть отрезок разбит на частей точками :
Сплайном -й степени называется функция, представляющая собой многочлен не выше -й степени на каждом из последовательно примыкающих друг к другу интервалов , причем в точках стыка двух интервалов функция непрерывна вместе со своими производными до порядка не выше .
Например, непрерывная кусочно-линейная функция (ломаная) является сплайном первой степени с производной, терпящей разрыв в точках излома.
Пусть на отрезке определена функция , значения которой в точках равны .
Задача интерполяции функции на отрезке (сплайном третьей степени) состоит в нахождении функции , равной многочлену третьей степени на каждом отрезке , т. е.
причем значения сплайна в узлах интерполяции равны соответствующим значениям заданной функции и сплайн-функция непрерывна в узлах интерполяции вместе с производными первого и второго порядков
Условия (2.17) — (2.20) дают линейных алгебраических уравнений для определения неизвестных коэффициентов при соответствующих степенях в многочленах .
Можно показать, что интерполяционный кубический сплайн для функции существует и является единственным, если вместе с уравнениями (2.17)-(2.20) удовлетворяется какая-либо пара дополнительных условий (краевых условий) следующего типа
Рассмотрим случай разбиения отрезка на равных частей с шагом , для которого и . Разберем построение интерполяционного кубического сплайна отдельно для условий I и II типов.
При построении сплайна, удовлетворяющего краевым условиям I типа, введем величины , называемые иногда наклонами сплайна в точках (узлах) .
Интерполяционный кубический сплайн вида
удовлетворяет условиям (2.17), (2.18), (2.19) для любых . Из условий (2.20) и краевых условий I типа можно определить параметр .
Действительно, легко проверить, что .
Кроме того, вычисления показывают, что
Если учесть, что
а также краевые условия I типа и условия (2.20), то получим систему из линейных уравнений относительно неизвестных
Решение этой системы позволяет найти значения неизвестных и определить интерполяционный сплайн в виде соотношений (2.21).
Матрица А системы (2.21) имеет порядок и является трехдиагональной
Метод Гаусса (метод исключения неизвестных) для системы (2.22) значительно упрощается и носит название метода прогонки. Прямой прогонкой находят так называемые прогоночные коэффициенты
Обратной прогонкой последовательно определяют неизвестные
Пример. На отрезке построить кубический сплайн с шагом , удовлетворяющий на концах отрезка краевым условиям I типа и интерполирующий функцию . С помощью интерполяционной формулы вычислить приближенное значение и сравнить его с точным.
Решение. Будем искать кубическую параболу , удовлетворяющую следующим условиям на концах отрезка и
Подставим значения в формулу (2.18) и получим сплайн вида
Тогда (точное значение равно 0.5).
При построении сплайна, удовлетворяющего краевым условиям II типа, введем величину — значение второй производной сплайна в узле .
Уравнения (2.17), (2.18), (2.20) будут удовлетворены, если интерполяционный кубический сплайн представить в виде
и используя краевые условия II типа и условия (2.19), получим систему из линейных уравнений относительно неизвестных
Системы (2.22) и (2.24) являются частными случаями системы линейных алгебраических уравнений следующего вида
Для функции , имеющей на отрезке непрерывные производные до третьего порядка включительно, точность интерполяции ее кубическим сплайном по точкам равномерного разбиения отрезка с шагом при любых указанных ранее краевых условиях оценивается следующим неравенством для любых на отрезке
Неравенство (2.26) дает завышенную оценку точности приближения функции сплайном в точке.
Интерполяционные полиномы Лагранжа и Ньютона
где значения хп и уп заданы. Такой полином определяется единственным образом и может быть представлен в следующем виде:
где сомножители 1п(х) задаются как
Полином (8.2) называется интерполяционным полиномом Лагранжа степени N. Для вычисления значения интерполяционного полинома в некорой точке х, не прибегая к явному вычислению коэффициентов полинома, применяется схема Невилла (Neville). На первом этапе зададим значения
Далее значения вспомогательных полиномов в точке х определяются согласно рекуррентному соотношению
и окончательно мы получаем
Пример 8.1 (интерполяционный полином Лагранжа)
В табл. 8.1 представлены значения уп в соответствующих точках хп и требуется провести гладкую кривую через эти точки.
Определим сомножители 1п(х), которые для данного набора данных записываются как
Тогда интерполяционный полином Лагранжа имеет вид
Предположим, что нам нужно вычислить значение Р$(х) в промежуточной точке х = 3, тогда схема Невилла (8.3) представляется в виде следующей таблицы:
Pf >(.v), я = 0. 3 — к
Рис. 8.1. Гладкая кривая, описывающая данные из табл. 8.1:
О — данные;—интерполяционный полином Лагранжа
Представление интерполяционного полинома в форме (8.2) имеет простой вид, и поэтому оно часто используется для теоретического анализа. Однако с точки зрения практических вычислений такое представление применимо только для небольших значений N. С увеличением N сомножители 1п(х) возрастают и начинают сильно осциллировать. Существует другое представление интерполяционного полинома, которое лучше подходит для вычислений, и оно основано на разделенных разностях. Пусть заданы различные точки х0. xN и значения в этих точках у0, . уjV. Разделенные разности нулевого порядка есть просто значения уп:
Разделенные разности порядка k в точке хп вычисляются но следующим рекурентным формулам:
Тогда, используя введенные обозначения, интерполяционный полином, удовлетворяющий условию (8.1), записывается в виде
Полином (8.4) называется интерполяционным полиномом Ньютона. Его можно представить и в другой форме:
Используя это представление, значение интерполяционного полинома Ньютона в некоторой точке х может быть вычислено простым и эффективным способом. Этот алгоритм называется схемой Горнера (Horner) и задается следующими рекурентными формулами:
Тогда требуемое значение полинома равно последнему элементу этой последовательности, т.е. Р^(х) = я0.
Пример 8.2 (интерполяционный полином Ньютона) Рассмотрим данные, представленные в табл. 8.1. Разделенные разности, вычисленные на основе этих данных, показаны в табл. 8.2.
Методы функциональной интерполяции
Пусть на множестве [math]\Omega=[a,b][/math] задана сетка [math]\Omega_n= \
В некоторых случаях [math]y_i=f(x_i),
i=\overline<0,n>[/math] является сеточным представлением заданной формульной функции [math]y=f(x)[/math] . Сеточная функция может задаваться совокупностью пар: [math](x_0,y_0),(x_1,y_1),\ldots,(x_n,y_n)[/math] .
Требуется найти функцию [math]y=F(x)[/math] , принимающую в точках [math]x_0,x_1,\ldots, x_n[/math] те же значения, что и функция [math]y_i=f(x_i),
i=\overline<0,n>[/math] , то есть [math]F(x_i)=y_i,
Точки [math]x_0,x_1,\ldots, x_n[/math] называются узлами интерполяции, а искомая функция [math]y=F(x)[/math] — интерполирующей.
Геометрически это означает, что нужно найти кривую, проходящую через заданное множество точек [math](x_i,y_i),
Одной из целей задачи интерполяции является вычисление значения функции в произвольной точке [math]x_<\ast>[/math] . При этом различаются собственно интерполирование, когда точка [math]x_<\ast>\in [x_0,x_n][/math] и экстраполирование, когда [math]x_<\ast>\notin [x_0,x_n][/math] .
Заметим, что можно провести бесчисленное множество "плавных" кривых, проходящих через заданное множество точек. Поэтому задача интерполяции в общей постановке не имеет единственного решения.
Если в качестве интерполирующей функции выбрать алгебраический многочлен, степень которого связана с числом заданных узлов интерполяции (на единицу меньше), решение задачи является единственным. Покажем это. Воспользуемся сначала кусочным способом. Выделим из отрезка [math][x_0,x_n][/math] частичный отрезок [x_i,x_
В качестве интерполирующей функции выберем алгебраический многочлен k-й степени (степень многочлена на единицу меньше количества узлов):
Будем искать неизвестные коэффициенты [math]a_0,a_1,\ldots,a_k[/math] из условия интерполяции (4.3) , т.е. [math]\delta \widetilde
Доказательство. Запишем условия интерполяции (4.9) с учетом (4.8) и обозначения [math]y_i= f(x_i)= f_i\colon[/math]
Эта система линейных алгебраических уравнений относительно коэффициентов [math]a_0,a_1, \ldots,a_k[/math] имеет единственное решение, так как определитель матрицы системы
не равен нулю (доказательство последнего факта содержится в курсе линейной алгебры, где этот определитель называется определителем Вандермонда). Следовательно, задача интерполяции также имеет единственное решение.
k=n[/math] , приходим к глобальному способу решения поставленной задачи. А именно, если задана сеточная функция в (n+1)-м узле [math]x_0,x_1,\ldots,x_n[/math] и требуется найти алгебраический многочлен с использованием условий интерполяции, то единственным решением задачи интерполяции является интерполяционный многочлен n-й степени:
коэффициенты которого находятся из системы (4.10) при [math]i=0,
1. Матрица (4.11) имеет дискретный (точечный) характер, так как ее элементы вычисляются по дискретным значениям [math]x_i,x_
2. При решении поставленной задачи предполагается, что исходная сеточная функция задана своими точными значениями, хотя класс задач, для которых используются такие функции, ограничен.
3. Имеются и другие формы записи интерполяционных многочленов. По теореме 4.1 все эти многочлены степени [math]n[/math] , удовлетворяющие функциональным условиям интерполяции и построенные по одним и тем же точкам, являются одним многочленом, записанным в разных формах.
4. При большом числе узлов решение системы (4.10) затруднительно. Искомый интерполяционный многочлен можно построить, не решая этой системы. Многочлены могут быть построены так, чтобы в самой структуре формулы многочлена условие интерполяции учитывалось.
5. При решении задачи функциональной интерполяции и в ее приложениях требуется:
а) выбрать наиболее удобную форму и степень интерполяционного многочлена. При этом можно использовать многочлены Лагранжа или Ньютона, а также формулу (4.8);
б) оценить погрешность интерполяции;
в) определить значения функции в точках, не совпадающих с узлами;
г) вычислить значения производных или определенных интегралов с использованием полученных интерполяционных многочленов.
Методика решения задачи интерполяции
1. По заданной сеточной функции составить интерполяционный многочлен определенной степени. При выборе степени многочлена следует руководствоваться желаемой точностью интерполяции.
2. Вычислить значения интерполяционного многочлена в заданных точках [math]x_<\ast j>,
j=1,\ldots,p[/math] , путем их подстановки в формулу многочлена.
Многочлен Лагранжа
Пусть исходная сеточная функция задана в (n+1)-й точках сетки [math]\Omega_n\colon\, y_i=f(x_i),\, i=\overline<0,n>[/math] , где [math]x_i\in [a,b]=[x_0,x_n][/math] — в общем случае неравноотстоящие узлы, определяемые шагами [math]h_
Воспользуемся сначала кусочным способом. Здесь и далее будем использовать обозначение [math]f_i=f(x_i)[/math] .
Выделим "окно" или частичный отрезок [math][x_i,x_
Действительно, легко убедиться в том, что [math]L_1(x)[/math] — алгебраический многочлен первой степени, который удовлетворяет условиям функциональной интерполяции (4.3) , т.е. [math]L_1(x_i)=f_i,
Выделим "окно" в виде двойного частичного отрезка [math][x_i,x_
Легко проверить, что (4.14) — многочлен второй степени и также удовлетворяет условиям функциональной интерполяции:
Обобщив запись многочлена на "окно" для к -кратного частичного отрезка [math][x_i,x_
где [math]P_
Легко проверить, что [math]P_
Если положить [math]i=0,\, k=n[/math] , то приходим к глобальному способу решения задачи. Тогда интерполяционный многочлен Лагранжа n-й степени имеет вид
где коэффициенты Лагранжа [math]P_
Очевидно, многочлен [math]L_n(x)[/math] , заданный равенством (4.17), является многочленом степени [math]n[/math] и удовлетворяет функциональным условиям интерполяции (4.3): [math]L_n(x_i)=f_i,
i=\overline<0,n>[/math] . Для записи интерполяционного многочлена Лагранжа удобно пользоваться табл. 4.1.
4.1>>\\\hline x-x_0& x_0-x_1& x_0-x_2& \cdots& x_0-x_n& D_0& f_0 \\\hline x_1-x_0& x-x_1& x_1-x_2& \cdots& x_1-x_n& D_1& f_1 \\\hline x_2-x_0& x_2-x_1& x-x_2& \cdots& x_2-x_n& D_2& f_2 \\\hline \vdots& \vdots& \vdots& \ddots& \vdots& \vdots& \vdots \\\hline x_n-x_0& x_n-x_1& x_n-x_2& \cdots& x-x_n& D_n& f_n \\\hline \multicolumn<5><|c|><\begin
Здесь [math]D_i[/math] — произведение элементов i-й строки, [math]\Pi_
i=\overline<0,n>[/math] . Тогда многочлен Лагранжа может быть записан в форме
1. Если заданная сеточная функция такая, что [math]f_i=1
(i=\overline<0,n>)[/math] , то из (4.17) следует, что [math]L_n(x)\equiv1[/math] , и поэтому справедливо равенство [math]\textstyle<\sum\limits_
2. Коэффициенты Лагранжа [math]P_
(i=\overline<1,n>)[/math] , то коэффициенты Лагранжа для всех исходных функций подсчитываются только один раз.
3. При введении дополнительных узлов интерполяции все коэффициенты многочлена Лагранжа необходимо пересчитывать заново, что неудобно на практике. От этого недостатка свободны многочлены Ньютона.
Перейдем к рассмотрению примеров решения задачи интерполяции на основе вышеизложенной методики.
Построить многочлен Лагранжа третьей степени для сеточной функции, заданной табл. 4.2. Вычислить значение функции в точке [math]x_<\ast>=2,\!5[/math] . Для записи многочлена использовать формулу (4.18).
Решение. 1. Построим многочлен Лагранжа. Для этого составим табл. 4.3, соответствующую табл. 4.1.
По формуле (4.18) получаем:
2. Вычислим значение функции в заданной точке: [math]L_3(2,\!5)= 4,\!8125[/math] .
Погрешность интерполяции многочленами Лагранжа
При определении значения [math]f(x),
x\ne x_[/math] для функции [math]y_i=f(x_i)
(i=\overline<0,n>)[/math] с помощью многочлена Лагранжа возникает погрешность или остаточное слагаемое [math]R_n(x)\colon[/math]
Здесь предполагается, что используется глобальный способ интерполяции и что [math]f(x)\in C_
На основе указанных предположений доказано, что при интерполяции функции [math]y_i=f(x_i)[/math] , заданной в общем случае на неравномерной сетке [math]\Omega_n[/math] , интерполяционным многочленом Лагранжа [math]\textstyle
где [math]\omega_n(x)= (x-x_0)(x-x_1)\ldots (x-x_n)[/math] — многочлен (n+1)-й степени, а [math]\xi\in (a,b)[/math] . Поскольку точно найти [math]R_n(x)[/math] нельзя (из-за неопределенности точки [math]\xi[/math] ), то при проведении вычислений обычно находятся только приближенные оценки погрешностей интерполяции, которые являются априорными.
Оценка погрешности интерполяции в некоторой произвольной фиксированной точке [math]x_<\ast>\in[a,b][/math] имеет вид, где [math]M_
Оценка максимальной погрешности интерполяции в любой точке [math]x\in[a,b][/math] , т.е. на всем отрезке [math][a,b]\colon[/math]
Замечание. Для сеточных функций с фиксированными узлами сетки (узлами интерполяции) также можно проводить оценки погрешности по формулам (4.21), (4.22), однако для этого необходимо численно определять [math]M_
С какой точностью можно вычислить значение [math]f(x)= \Bigl. <\sqrt
Решение. Так как требуется вычислить погрешность только в одной точке [math]x_<\ast>=85[/math] , то необходимо использовать формулу (4.21).
Найдем оценку погрешности в точке [math]x_<\ast>=85[/math] для [math]L_1(x)[/math] . В качестве "окна" линейной интерполяции [math](n=1)[/math] выбираем отрезок [math][x_1,x_2]=[36;100][/math] . Определим величину [math]M_2\colon[/math]
Найдем оценку погрешности в точке [math]x_<\ast>[/math] для [math]L_2(x)[/math] . В качестве "окна" квадратичной интерполяции [math](n=2)[/math] выбираем отрезок [math][x_0,x_2]= [16;100][/math] . Определим величину [math]M_3\colon[/math]
Погрешность для многочлена [math]L_2(x)[/math] получилась больше погрешности для многочлена [math]L_1(x)[/math] , что обусловлено величиной [math]\omega_2(x_<\ast>)[/math] .
Сходимости функционального интерполяционного процесса для непрерывных функций
Как отмечалось выше, выбор сетки и соответствующей степени интерполяционного многочлена при интерполяции сеточных функций является одной из важных задач, решить которую можно, рассмотрев проблему сходимости интерполяционных процессов [6,40] для непрерывных функций [math]y=f(x)\in C_
Будем считать, что интерполяция проводится на последовательности сеток
с возрастающим разбиением [math]k[/math] отрезка [math][a,b]\colon\, k_1=1; k_2=2; \ldots; k_n=n[/math] и т.д.
Если при данных разбиениях при возрастающем [math]k=1,2,\ldots,n,\ldots[/math] определяются значения [math]L_k(x_<\ast>)[/math] в некоторой промежуточной точке [math]x_<\ast>[/math] , то реализуется интерполяционный процесс, характеризующийся последовательностью значений многочленов: [math]L_1(x_<\ast>), L_2(x_<\ast>), \ldots, L_n(x_<\ast>), \ldots[/math] .
Интерполяционный процесс для функции [math]f(x)[/math] сходится в точке [math]x_<\ast>\in [a,b][/math] , если существует [math]\lim\limits_
Для отрезка [math][a,b][/math] существует понятие равномерной сходимости в некоторой норме, например max, [math]\max_ <\substack
Характер сходимости или расходимости интерполяционного процесса зависит как от гладкости и поведения функции [math]f(x)[/math] , так и от выбора последовательности сеток [math]\Omega_n
(n=1,2,\ldots)[/math] . Так, показано, что если f(x) непрерывна на [math][a,b][/math] , то найдется такая последовательность [math]\Omega_n
(n=1,2,\ldots)[/math] , для которой интерполяционный процесс сходится равномерно на [math][a,b][/math] (теорема Марцинкевича). Однако для дискретных функций, рассматриваемых в данном разделе, эта теорема не применима. Отметим также, что построены расходящиеся интерполяционные процессы и для формульных функций, например [math]f(x)=|x|,
x\in[-1;1][/math] . Кроме того, применение многочленов высоких степеней приводит к так называемым "провалам" между узлами интерполяции, часто называемым осцилляциями. Указанные свойства интерполяционных процессов обусловливают нецелесообразность применения интерполяционных многочленов высоких степеней. В связи с этим в вычислительной практике для сеточных функций степень [math]n[/math] не берут выше [math]5\div 8
(n\leqslant 5\div 8)[/math] и задание частичного отрезка согласуют с выбранной степенью многочлена.
Рассмотрим часто использующиеся на практике линейную и параболическую интерполяцию.
Линейная и параболическая интерполяция с помощью многочлена Лагранжа
В прикладных расчетах часто применяется простейшая кусочная интерполяция, основанная на многочленах первой степени [math]L_1(x)[/math] или второй степени [math]L_2(x)[/math] . В этом случае функциональная интерполяция называется линейной или параболической (квадратичной) соответственно. Рассмотрим возможные способы их реализации и найдем оценки их погрешностей.
Пусть для сеточной функции [math]y_i=f(x_i)[/math] , заданной на сетке [math]\Omega_n= \
(i=\overline<0,n>)[/math] , и оценить погрешности. Для обоснованного выбора степени интерполяционного многочлена необходимо указать, какую погрешность имеют значения исходной функции [math]f(x_i)[/math] в узлах. Если эта погрешность составляет величину [math]O(h_
Методика решения задачи линейной интерполяции
1. По расположению заданной точки [math]x_<\ast>[/math] на оси [math]Ox[/math] выбрать из всех частичных отрезков [math][x_0,x_1], [x_1,x_2],\ldots, [x_i,x_
2. Для отрезка [math][x_i,x_
3. Вычислить искомое значение [math]f(x_<\ast>)[/math] согласно (4.13):
Как показано ниже, порядок этой аппроксимации равен двум, т.е. [math]O(h_
Геометрическая интерпретация линейной интерполяции при известной формульной функции [math]y=f(x)[/math] (штриховая линия) изображена на рис. 4.3. Здесь прямая [math]AB[/math] соответствует графику функции [math]y=L_1(x)[/math] на отрезке [math][x_i, x_
Методика решения задачи параболической интерполяции
1. Из всей совокупности спаренных частичных отрезков [math][x_0,x_2], [x_1,x_3], \ldots, [x_
O_2\equiv [x_i, x_
3. Используя значения коэффициентов Лагранжа, вычислить значения [math]L_2^<(1)>(x_<\ast>),\, L_2^<(2)>(x_<\ast>)[/math] по формуле (4.14).
Если в расчетах не требуется высокая точность интерполирования, то можно ограничиться выбором одного "окна", например [math]O_2[/math] , и тогда [math]f(x_<\ast>)\approx L_2^<(2)>(x_<\ast>)[/math] . Для достижения повышенной точности интерполяцию провести для двух "окон" и результаты [math]L_2^<(1)>(x_<\ast>),\, L_2^<(2)>(x_<\ast>)[/math] усреднить:
При этом порядок интерполяции повышается на единицу, т.е. [math]L_<2\text
Замечание. Если [math]x_<\ast>\in [x_0,x_2][/math] или [math]x_<\ast>\in [x_
Геометрическая интерпретация параболической интерполяции изображена на рис. 4.4. Параболе [math]y= L_2^<(1)>(x_<\ast>)[/math] соответствует кривая [math]A_1A_2B_1[/math] , параболе [math]y= L_2^<(2)>(x_<\ast>)[/math] — кривая [math]A_2B_1A_2[/math] . Точка [math]C[/math] соответствует значению [math]L_2^<(1)>(x_<\ast>)[/math] , точка [math]D[/math] — значению [math]L_2^<(2)>(x_<\ast>)[/math] точка [math]E[/math] — значению [math]L_<2\text
Приведем оценки погрешностей линейной и параболической интерполяции. Вначале предположим, что сетка [math]\Omega_n[/math] равномерная (это имеет значение только для параболической интерполяции). Формулы (4.13), (4.14) и оценки (4.22), записанные для "окон" интерполяции, упрощаются, если ввести в рассмотрение новую переменную — фазу интерполяции и [math]u=\frac
Здесь учтено, что [math]x-x_
x-x_
Таким образом, реализуются следующие оценки погрешностей линейной и параболической интерполяции, справедливых для соответствующих "окон":
M_<3,i>= \max_<[x_i, x_
f(x)\in C_3[a,b][/math] соответственно для линейной и параболической интерполяции. Таким образом, из оценок (4.23), (4.24) следует, что линейная интерполяция обеспечивает на частичном отрезке [math][x_i, x_
1. Оценка (4.23) для линейной интерполяции инвариантна по отношению к виду сетки [math]\Omega_n[/math] (равномерной или неравномерной). Параболическая интерполяция также сохраняет указанную погрешность при выполнении интерполяции на неравномерной сетке [math]\Omega_n[/math] . В некоторых источниках для [math]L_2(x)[/math] при выполнении условия [math]f(x)\in C_3[a,b][/math] приведена оценка
2. Если гладкость функции [math]f(x)[/math] не достигает вышеуказанной и класс гладкости понижен на единицу [math]\bigl(f(x)\in C_2[a,b]\bigr)[/math] , то порядок параболической интерполяции также понижается на единицу:
3. Повышение класса гладкости функции [math]f(x)[/math] выше [math]C_3[a,b][/math] не приводит к увеличению порядка параболической интерполяции относительно Н , т.е. происходит как бы его "замораживание" или "насыщение". Указанные свойства, вытекающие из оценок (4.25), (4.26), носят общий характер и имеют место в других оценках аппроксимации многочленов, сплайнов, производных и интегралов.
4. Для произвольной степени интерполяционного многочлена при [math]h=text
f(x)\in C_
\omega_n(x)= u(u-1)\cdot \ldots\cdot (u-n)[/math] . Для многочленов с [math]n \leqslant 5[/math] величины [math]\Omega^n[/math] или их оценки являются такими:
Из (4.27) следует, что на отрезке [math][x_0,x_n][/math] величина [math]\|R_n(x)\|[/math] есть [math]O(h^
5. В широком классе задач математической физики применяются расчетные схемы в основном второго (и иногда третьего и выше) порядка точности. При их реализации, как правило, используются встроенные интерполяционные алгоритмы, основанные на многочленах и сплайн-функциях. Степени интерполяционных многочленов при этом должны выбираться из условия соответствия порядков их аппроксимации порядкам точности схем. Если эти порядки одинаковы, то порядок точности схем сохраняется, хотя константа в оценке погрешности схемы изменяется. Если же порядок встроенных интерполяционных алгоритмов хотя бы на единицу выше порядка точности схемы, то вместе с порядком точности схемы сохраняется и указанная константа. Отсюда следует, что необходимо выбирать такую степень интерполяционного многочлена, которая либо обеспечивает равенство порядка аппроксимации порядку точности схемы, либо на единицу превышает последний. Таким образом, использование параболической интерполяции в качестве встроенных алгоритмов или для восполнения численных решений, полученных по схемам второго порядка, позволяет сохранить требуемую точность расчета, а также не дает избыточный порядок и, следовательно, не усложняет алгоритм. Это замечание носит общий характер и относится к любым аппроксимационным алгоритмам, выполняющим функцию восполнения или интерполирования.
Перейдем к рассмотрению примеров решения задач линейной и параболической интерполяции.
Дана сеточная функция, являющаяся сеточным представлением формульной функции [math]f(x)=x^3[/math] (табл. 4.4). Найти значение [math]f(x_<\ast>)[/math] при [math]x_<\ast>=2[/math] с помощью линейной и параболической интерполяции.
Решение. Применение линейной интерполяции. Воспользуемся методикой.
1. Выбираем "окно" интерполяции [math][x_i,x_
2. Вычислим коэффициенты Лагранжа: [math]P_<1\,2>= \frac
P_<1\,3>= \frac
3. Определим искомое значение [math]f(2)[/math] по формуле (4.13):
Применение параболической интерполяции с осреднением. Также воспользуемся соответствующей методикой.
1. Выбираем "окна" интерполяции [math]O_1\equiv [x_1;x_3]= [0;3];
O_2\equiv [x_2;x_4]= [1;4][/math] .
Вычисления выполнены правильно, так как [math]\textstyle<\sum\limits_
Искомое значение получим в результате осреднения: [math]f(2)\approx \frac
Многочлены Ньютона
Разделенные и конечные разности . В практике функционального интерполирования иногда удобнее использовать многочлены Ньютона, степень которых можно последовательно повышать путем добавления очередных слагаемых, имеющих более высокую степень. Такие несимметричные многочлены, альтернативные симметричным многочленам Лагранжа, основаны на разделенных и конечных разностях, вычисляемых по интерполируемой сеточной функции.
Разделенные разности вводятся для функции [math]y_i= f(x_i)=f_i,
i=\overline<0,n>[/math] , заданной на неравномерной сетке [math](h_
i=\overline<0,n>[/math] , определенной на равномерной сетке [math](h_
Выбрав внутри неравномерной или равномерной сетки соответствующие шаблоны интерполяции [math](x_i,x_
– разделенная разность первого порядка: [math]f(x_i,x_
– конечная разность первого порядка: [math]\Delta f_i= f_
– конечная разность второго порядка: [math]\Delta^2f_i= \Delta (\Delta f_i)= \Delta f_
Последовательность получения разделенных и конечных разностей при [math]k=3[/math] для произвольной функции наглядно представляют табл. 4.5 и 4.6.
Для гладких функций числовые значения [math]f(x_j, x_
Связь между разделенными и конечными разностями k-го порядка при [math]h=\text
Действительно, при [math]k=1[/math] для разделенной разности [math]f(x_i,x_
Пусть [math]k=2[/math] . Тогда получаем
Таким образом, связь (4.28) выполняется и при [math]k=2[/math] . Справедливость (4.28) при произвольном к можно доказать методом математической индукции.
Интерполяционный многочлен Ньютона для неравномерной сетки
Пусть исходная (интерполируемая) сеточная функция [math]y_i=f(x_i),
i=\overline<0,n>[/math] , задана на неравномерной сетке [math]\Omega_n\equiv \
Воспользуемся сначала кусочным способом. Из всей совокупности узлов выбираем шаблон [math](x_i,x_
Тогда для функциональной интерполяции может быть использован многочлен Ньютона, основанный на разделенных разностях:
Действительно, [math]N_k(x)[/math] — многочлен k-й степени, что определяется сомножителями последнего слагаемого (разделенные разности, входящие в качестве одного из сомножителей в эти произведения, есть числа). Кроме того, для многочлена [math]N_k(x)[/math] удовлетворяются функциональные условия интерполяции: [math]N_k(x_j)= f_j,
j=i,\ldots,i+k[/math] . Проверим их справедливость при [math]k=1[/math] (шаблон [math](x_i,x_
Таким образом, условия интерполяции для [math]N_1(x)[/math] выполнены, следовательно, многочлен (4.30) может быть использован для линейной интерполяции кусочным способом. Пусть [math]k=2[/math] . Тогда
Таким образом, условия интерполяции для многочлена [math]N_2(x)[/math] также выполнены и он может использоваться для параболической интерполяции кусочным способом.
Для произвольного [math]k[/math] справедливость равенств [math]N_k(x_j)=f_i,
j=\overline
k=n[/math] , приходим к глобальному способу. Тогда интерполяционный многочлен Ньютона n-й степени имеет вид
2. Интерполяционный многочлен Ньютона (4.29) или (4.32) (так же, как и многочлен Ньютона, выражаемый ниже через конечные разности) записан не через значения функции, как это имеет место для многочлена Лагранжа, а через разделенные разности. Поэтому при изменении степени [math]k[/math] в процессе интерполирования у многочлена Ньютона [math]N_k(x)[/math] требуется только добавить или отбросить соответствующее число слагаемых. Это иногда упрощает алгоритм интерполирования.
3. При интерполяции на основе (4.29) или (4.32) узлы интерполяции [math]x_i,x_
4. Остаточное слагаемое многочлена (4.32) совпадает с остаточным слагаемым многочлена Лагранжа, и оценки (4.21), (4.22), справедливые для точки [math]x_<\ast>[/math] и всего отрезка [math][x_0,x_n][/math] , сохраняются.
Интерполяционные многочлены Ньютона для равномерной сетки
Сначала рассмотрим решение задачи кусочной интерполяции (применение кусочного способа). Если функция [math]y_i= f(x_i),
i=\overline<0,n>[/math] задана на равномерной сетке [math]\Omega_n[/math] , характеризующейся [math]h_
где [math]q=\frac
(j=1,2,\ldots,k)[/math] — конечные разности.
В соответствии с формулой для фазы интерполяции [math]q[/math] точка начала ее отсчета расположена в узле [math]x_i[/math] и входящие в (4.33) конечные разности относятся к этой же точке [math]x_i[/math] . В связи с этим (4.33) удобно применять в начале выделенного шаблона [math](x_i,x_
(q<0)[/math] . Поэтому многочлен [math]N_k^<(I)>(q)[/math] называется интерполяционным многочленом для интерполяции вперед (в начале таблицы ) или для экстраполяции назад. В (4.33) этот многочлен обозначен цифрой I , указанной в скобках вверху. Для определенности назовем его первым интерполяционным многочленом Ньютона.
k=n[/math] , получаем решение задачи глобальной интерполяции на всем отрезке [math][x_0,x_n]\colon[/math]
где [math]q=\frac
где [math]\xi\in(a,b)[/math] — некоторое промежуточное значение между узлами [math]x_0,x_1,\ldots,x_n[/math] и точкой [math]x[/math] .
Если фазу интерполяции определить относительно [math]x_= \frac
Данный многочлен удобно применять в конце выделенного шаблона или всей таблицы [math]y_i= f(x_i),
Если [math]i+k=n[/math] , а [math]x_<\ast>>x_n[/math] , то (4.35) используется для экстраполяции правее точки [math]x_n
(\widehat)[/math] называется многочленом для интерполяции назад (в конце таблицы) или экстраполяции вперед.
i+k=n[/math] , получаем решение задачи глобальной интерполяции — второй интерполяционный многочлен Ньютона n-й степени (где [math]\widehat= \frac
Остаточное слагаемое многочлена (4.36) имеет вид
Схема выбора узлов интерполяции при изменении степеней интерполяционных многочленов [math]N_
(k=1,2,3)[/math] в одном алгоритме показана на рис. 4.5. Если точка [math]x_<\ast>[/math] находится в начале отрезка [math][x_0,x_n][/math] , например на частичном отрезке [math][x_0,x_1][/math] , то применяется первый интерполяционный многочлен Ньютона, а если в конце, например на частичном отрезке [math][x_
1. Из формул интерполяционных многочленов [math]N_)[/math] видно, что повышение их степеней в процессе реализации алгоритма не требует пересчета предыдущих слагаемых, входящих в многочлены с меньшими степенями.
2. Для гладких функций при повышении порядка конечных разностей справедливо свойство: [math]\Delta^kf_i\to0[/math] при [math]k\to\infty[/math] , и поэтому, как только очередное слагаемое в рассматриваемых многочленах становится меньше требуемой точности интерполяции, увеличение степени [math]N_
N_)[/math] следует прекратить. Это замечание справедливо также и для [math]N_k(x)[/math] .
Пример 4.4. Построить многочлен Ньютона третьей степени для сеточной функции, заданной табл. 4.7. Вычислить значение функции в точке [math]x_<\ast>=2,\!5[/math] . Решить задачу интерполяции при включении одного дополнительного значения сеточной функции: [math]f(1)=5[/math] .
Построим многочлен Ньютона (4.32), справедливый для произвольного расположения узлов. Для этого составим табл. 4.8, аналогичную табл. 4.5:
По формуле (4.32) для [math]n=3[/math] имеем
Поскольку в данной задаче заданы равностоящие узлы, воспользуемся также формулой (4.34) для первого интерполяционного многочлена Ньютона:
где [math]q=\frac
Запишем также второй интерполяционный многочлен Ньютона (4.36):
Сравнивая с результатом примера 4.1, можно заключить, что [math]L_3(x)= N_3^<(I)>(x)= N_3^<(II)>(x)[/math] . Это еще раз подтверждает единственность решения задачи интерполяции в классе многочленов, удовлетворяющих условиям теоремы 4.1.
Предположим, что в ходе некоторого эксперимента получен новый результат [math]f(1)=5[/math] , дополняющий заданную сеточную функцию. Тогда для решения задачи интерполяции с помощью многочлена Ньютона можно использовать уже полученный результат. Для этого дополним табл. 4.8. Новый узел и соответствующее значение функции поместим в конце табл. 4.10.
По формуле (4.34) имеем
Легко проверить, что новый узел можно добавить и в начале таблицы, т.е. будет найден тот же многочлен [math]N_4(x)[/math] . На рис. 4.6 изображены полученные интерполяционные многочлены.
Вычислим значение функции в точке [math]x_<\ast>=2,\!5\colon[/math]
Пример 4.5. Для сеточной функции из примера 4.3 найти линейный и параболический многочлены Ньютона и на их основе подсчитать значение функции в точке [math]x_<\ast>=2,\!5[/math] . Оценить погрешность интерполяции.
Так как точка [math]x_<\ast>=2,\!5[/math] находится вблизи узла [math]x_i=x_2[/math] и шаблон [math](x_2,x_3,x_4)[/math] имеет неравномерность по шагу, то для выполнения интерполяции выбираем многочлен Ньютона (4.29).
1. Определим разделенные разности [math]f(x_i,x_
Сравнивая [math]N_2(x)[/math] с соответствующим многочленом Лагранжа L_2(x) (см. пример 4.3), видим, что они совпадают. Это подтверждает единственность решения задачи о построении интерполяционного многочлена.
Итак, линейная интерполяция приводит к отличию от точного значения на [math]\left|\frac<20-15,\!625><15,\!625>\right|\cdot 100\%=28\%[/math] , а квадратичная — на [math]\left|\frac<14,\!5-15,\!625><15,\!625>\right|\cdot 100\%=7\%[/math] .
Таким образом, найдена фактическая погрешность, которая для квадратичной интерполяции получилась в 4 раза меньше.
3. Найдем априорную оценку погрешности линейной интерполяции на отрезке [math][x_2,x_3][/math] по формуле (4.23). Для этого необходимо сначала оценить производную [math]M_<2,i>[/math] . При этом, как и в реальных задачах, будем определять [math]f''(x)[/math] численным дифференцированием. С этой целью используем многочлен [math]N_2(x)\colon\, N'_2(x)=16x-9,
N''_2(x)=16[/math] . Приближенно можно принять, что [math]M_<2,i>=16[/math] . Тогда получаем: [math]|R_1(x)|\leqslant \frac<2^2><8>\cdot16=8[/math] или [math]\frac<8><15,\!625>\cdot 100\%= 51,\!2\%[/math] . Заметим, что фактическая погрешность получилась примерно в два раза меньше априорной.
Пример 4.6. Для сеточной функции, заданной в примере 4.4, построить интерполяционные многочлены первой и второй степени для нахождения значений в точках
Для подсчета значения функции в точке [math]x_<\ast 1>=2,\!5[/math] (согласно рис. 4.5,а) выбираем шаблон [math](x_0,x_1,x_2)= (2;3;4)[/math] для параболической интерполяции [math](i=0,
k=2)[/math] и шаблон [math](x_0,x_1)= (2;3)[/math] Для линейной интерполяции [math](i=0,
k=1)[/math] . Тогда по формуле (4.34) имеем
где [math]q=\frac
Для подсчета значения в точке [math]x_<\ast 2>=4,\!5[/math] выбираем шаблон [math](3;4;5)= (x_1,x_2,x_3)[/math] для параболической интерполяции [math](i=1,
k=2)[/math] и шаблон [math](4;5)=(x_2,x_3)[/math] для линейной интерполяции [math](i=2,
k=1)[/math] . Тогда по формуле (4.35) имеем
где [math]\widehat=\frac
Многочлен [math]N_1^<(II)>(x)[/math] можно использовать для экстраполяции, т.е. для подсчета функции в точке [math]x_<\ast 4>=5,\!5\colon\, N_2^<(II)>(5,\!5)=5[/math] .
Чтобы подсчитать значения в точке [math]x_<\ast3>=3,\!5[/math] , сначала выберем шаблон [math](3;4;5)=(x_1,x_2,x_3)[/math] для параболической интерполяции [math](i=1,
k=2)[/math] и шаблон [math](3;4)=(x_1,x_2)[/math] для линейной интерполяции [math](i=1,
k=1)[/math] . Согласно рис. 4.5,б, применим формулу (4.33):
где [math]q=\frac
Подсчитать значения в точке [math]x_<\ast3>=3,\!5[/math] (согласно рис. 4.5,б) можно, использовав шаблон [math](2;3;4)= (x_0,x_1,x_2)[/math] для параболической интерполяции [math](i=0,
k=2)[/math] и шаблон [math](3;4)=(x_1,x_2)[/math] для линейной интерполяции [math](i=1,
Интерполяция полиномами Лагранжа и Ньютона
Пусть задана функция .
Пусть заданы точки из некоторой области .
Пусть значения функции известны только в этих точках.
Точки называют узлами интерполяции.
— шаг интерполяционной сетки.
Задача интерполяции состоит в поиске такой функции из заданного класса функций, что
Метод решения задачи
Полином Лагранжа
Представим интерполяционную функцию в виде полинома
где — полиномы степели n вида:
Очевидно, что принимает значение 1 в точке и 0 в остальных узлах интерполяции. Следовательно в точке исходный полином принимает значение
Таким образом, построенный полином является интерполяционным полиномом для функции на сетке .
Полином Ньютона
Интерполяционный полином в форме Лагранжа не удобен для вычислений тем, что при увеличении числа узлов интерполяции приходится перестраивать весь полином заново.
Перепишем полином Лагранжа в другом виде:
где — полиномы Лагранжа степени i ≤ n.
Пусть
. Этот полином имеет степень i и обращается в нуль при .
Поэтому он представим в виде:
, где — коэффициент при . Так как не входит в , то совпадает с коэффициентом при в полиноме . Таким образом из определения получаем:
Препишем формулу в виде
Рекуррентно выражая пролучам окончательную формулу для полинома:
Такое представление полинома удобно для вычисления, потому что увеличение числа узлов на единицу требует добавления только одного слагаемого.
Погрешность интерполирования
Поставим вопрос о том, насколько хорошо интерполяционный полином приближает функцию на отрезке [a,b].
Рассмотри м остаточный член:
, x ∈ [a, b].
По определению интерполяционного полинома
поэтому речь идет об оценке при значениях .
Пусть имеет непрерывную (n+1) производную на отрезке [a, b].
Тогда погрешность определяется формулой:
,
где ,
— точка из [a, b].
Так как точка наизвестна, то эта формула позволяет только оценить погрешность:
где
Из вида множетеля следует, что оценка имеет смысл только при . Если это не так, то при интерполяции используются полиномы низких степеней (n = 1,2).
Выбор узлов интерполяции
Так как от выбора узлов завист точность интерполяции, то возникает вопрос о том, как их выбирать. С помощью выбора узлов можно минимизировать значение в оценке погрешности. Эта задача решается с помощью многочлена Чебышева [1]:
В качестве узлов следут взять корни этого многочлена, то есть точки:
Пример
В качастве примера рассмотрим интерполяцию синуса. Возьмем равномерную решетку x = [-3,-1.5,0,1.5,3];
Интерполяция полиномом Лагранжа:
Ошибка(максимальное отклонение от sin(x) на отрезке):0.1423
Интерполяция полиномом Ньютона:
Ошибка:
Возьмем решетку x с узлами в корнях полинома Чебышева= [-2.8531,-1.7632,0,1.7634,2.8532];
Интерполяция полиномом Лагранжа:
Ошибка: 0.0944
Интерполяция полиномом Ньютона:
Ошибка: