Как найти погрешность интерполяции функции - Сайт, где вы сможете решить свои вопросы

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

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

Теорема.

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

(6)

где.

Доказательство.

Будем искать погрешность в виде
,
где–
функция, ограниченная на.
(При этом гарантируется, чтообращается в ноль в точках интерполяции).
Чтобы получить представление о,
рассмотрим вспомогательную функцию

где
– некоторое фиксированное значение.
Очевидно, что нафункция имеетнуля. Это узлы интерполяции и точка.
По теореме Ролля,.
Продифференцировав (*) (n+1) раз и подставив,
получим

Отсюда
.
Ясно, чтов теореме Ролля зависит от расположения
нулей функции, т.е.– некоторая неявная зависимость. Обозначаячерезполучим (6).

Следствие. Из формулы (6) следует
оценка интерполяции

где
.

Конкретная величина погрешности в точке
зависит, очевидно, от значения полинома
в этой точке. Качественный характеризображен на Рис. 2.

Рис. 2 Характер
.

За пределами отрезка интерполяции (т.е.
при экстраполяции)
быстро растет, экстремальные значения
меньше в середине отрезка интерполяции.
Для равноотстоящих узловдляимеет место

где

Поэтому на отрезке интерполяции

(8)

Оценка (8) – сильно завышенная оценка
ошибки. Для получения точной оценки
надо искать экстремумы
.

Оценка (7) для погрешности интерполяции
не является завышенной. Можно показать,
например, что она достигается при
интерполировании полиномом
-ой
степени полиномомстепени.

Пример.

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

Решение.

ИП Ньютона
.
Для коэффициентов:

В данном случае:

;

;;

Замечания о минимизации ошибки при
полиномиальной интерполяции.

,
где.

Ясно, что величина ошибки зависит от
расположения узлов
.
Если узлы можно выбирать, то ситуация
значительно улучшается. Можно
минимизировать ошибку интерполяции.

Полиномы Чебышева. Интерполяция по Чебышевским узлам.

Определение. На отрезкеопределим многочлены Чебышева:

(1)

Найдем несколько первых многочленов:

Т.к.
,
то

(2)

Полагая в формуле (2) ,
получим

(3)

Получена рекуррентная формула для
полиномов Чебышева. Отсюда следует, что
‑полином–
ой степени. Последовательно получаем:

;

и т.д.

Свойства многочленов Чебышева.

Системаортогональна на отрезкес весом.
Норма
–
четные функции;–
нечетные функции.
Коэффициент при старшей степени
многочленаравен(доказательство по индукции).
Многочлен
имеет на интервалеровноразличных действительных корней,
определяемых формулой:,.
,
причем максимум достигается в точках.
При этом.
Из определения следует, что.
При этом
Многочлен
среди всех многочленов–
ой степени с коэффициентом(при старшем члене)
обладает тем свойством, что(без доказательства).

Благодаря свойству 6. многочлен
называется многочленом, наименее
отклоняющимся от нуля. Полиномы Чебышева,
нормированные таким образом, чтобы
коэффициент при старшей степенибыл
равен 1:

;;;;;
и т.д.

Можно записать полином степени
,
наименее отклоняющийся от нуля на
произвольном отрезке.

Применение полиномов Чебышева к задаче
интерполяции.

Задача.Оптимизировать интерполяцию
полиномом Лагранжа с помощью выбора
узлов интерполяции.

Решение. Выбор узлов интерполяции.

Пусть.
Погрешность интерполяции,
причем,
где,
а полином–
многочлен степени,
с коэффициентамипри старшем члене;– его нули.

но по свойству 6. многочленов Чебышева,

Если узлы интерполяции выбрать как у
полиномов Чебышева, т.е. в точках
,
то нулиисовпадут, а так как в обоих многочленахпри старшем члене,
следовательно,и
достигается:

Вывод.

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

Геометрическая интерпретация корней
полинома Чебышева

Если верхнюю полуокружность единичного
радиуса разделить на
частей, то середины дуг – координаты
нулей, экстремумы – точки деления (см.
Рис. 3)

Рис. 3 Геометрическая интерпретация
корней полинома чебышева.

Равномерное приближение функций на
отрезке.

Пусть,.
Расстояние между двумя функциями.
Пусть–
достаточно гладкая, т.е..
Тогда найдется такое(степень интерполяционного полинома),
что выполняется условие:

Пример.

.
Приблизить многочленом–
ой степени так, чтобы выполнялось
условие:.
Найти.

Решение.

.
Возьмем для интерполяции интерполяционный
полином Лагранжа, построенный по нулям
полинома Чебышева.
Имеем:

где.

Производные функции
:

,
,
…
.

Тогда

Отсюда

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

Соседние файлы в папке Учебники

Источник

Макеты страниц

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

Теорема 11.4. Пусть функция дифференцируема раз на отрезке содержащем узлы интерполяции Тогда для погрешности интерполяции в точке справедливо равенство

в котором некоторая точка, принадлежащая интервалу

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

Следствие. В условиях теоремы 11.4 справедлива оценка погрешности интерполяции в точке имеющая вид

а также оценка максимума модуля погрешности интерполяции на отрезке имеющая вид

Здесь

Пример 11.4. Оценим погрешность приближений к полученных в примере 11.3 с помощью интерполяции многочленами первой и второй степени. В этих случаях неравенство (11.26) примет вид

Заметим, что для имеем Поэтому здесь и 1.5. Тогда в силу неравенств (11.28) и (11.29) получаем следующие оценки погрешности:

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

Источник

Материал из MachineLearning.

Перейти к: навигация, поиск

Содержание

1 Полином в каноническом виде
2 Способ вычисления полинома в точке
3 Анализ метода
- 3.1 Сложность вычислений
- 3.2 Погрешность интерполяции
- 3.3 Выбор узлов интерполяции
4 Вычислительный эксперимент
- 4.1 Пример: Интерполяция синуса
5 Рекомендации программисту
- 5.1 Предварительные настойки
- 5.2 Использование программы
6 Вывод
7 Прикреплённые файлы
8 Литература
9 Смотри также
10 Ссылки

Интерполяция — приближение одной функции другой функцией.

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

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

Полином в каноническом виде

Известно, что любая непрерывная на отрезке [a,b] функция f(x) может быть хорошо приближена некоторым полиномом P_n(x).
Справедлива следующая Теорема (Вейерштрасса): Для любого >0 существует полином P_n(x) степени , такой, что $max_{x in [a,b]}|f(x)-P_n(x)|<eps$

В качестве аппроксимирующей функции выберем полином степени n в каноническом виде:

Коэффициенты полинома c_i определим из условий Лагранжа , , что с учётом предыдущего выражения даёт систему линейных алгебраических уравнений с n+1 неизвестными:

$begin{matrix} c_0 + c_1x_0 + c_2x_0^2 + ldots + c_nx_0^n = y_0 \ c_0 + c_1x_1 + c_2x_1^2 + ldots + c_nx_1^n = y_1 \ ldots ldots ldots ldots ldots \ c_0 + c_1x_n + c_2x_n^2 + ldots + c_nx_n^n = y_n end{matrix}$

Обозначим систему таких уравнений символом (*) и перепишем её следующим образом:

$sum_{p=0}^n c_i^p = y_i, quad i=1, ldots, n$

или в матричной форме: где – вектор-столбец, содержащий неизвестные коэффициенты c_i , – вектор-столбец, составленный из табличных значений функции y_i , а матрица имеет вид:

$mathbf{A} = begin{pmatrix} 1 & x_0 & x_0^2 & ldots & x_0^n \ 1 & x_1 & x_1^2 & ldots & x_1^n \ ldots & ldots & ldots & ldots & ldots \ 1 & x_n & x_n^2 & ldots & x_n^n \ end{pmatrix}$

Система линейных алгебраических уравнений (*) относительно неизвестных c_i иметь единственное решение, если определитель матрицы отличен от нуля.

Определитель матрицы называют определителем Вандермонда, его можно вычислить по следующей формуле:

$mathbf{detA} = prod_{i,j=0 \ i neq j }^n (x_i - x_j) neq 0$

Число узлов интерполяционного полинома должно быть на единицу больше его степени. Это понятно из интуитивных соображений: через 2 точки можно провести единственную прямую, через 3 – единственную параболу и т.д. Но полином может получиться и меньшей степени. Т.е. если 3 точки лежат на одной прямой, то через них пройдёт единственный полином первой степени (но это ничему не противоречит: просто коэффициент при старшей степени равен нулю).

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

Из-за плохой обусловленности матрицы рекомендуется применять другие методы интерполяции (например, интерполяция полиномами Лагранжа). При этом важно понимать, что при теоретическом применении различных методов они приводят к одинаковому результату, т.е. мы получим один и тот же полином.

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

Способ вычисления полинома в точке

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

Первый способ. Можно посчитать значение a₁x и сложить с a₀. Далее найти a₂x², сложить с полученным результатом, и так далее. Таким образом, на j-ом шаге вычисляется значение a_jx^j и складывается с уже вычисленной суммой $a_0 + a_1x + ldots + a_{j-1}x^{j-1}$ .

Вычисление значения a_jx^j требует j операций умножения. Т.е. для подсчёта многочлена в заданной точке требуется (1 + 2 + … + n) = n(n+1)/2 операций умножения и n операций сложения. Всего операций в данном случае: Op₁ = n(n+1)/2 + n.

Второй способ. Полином можно также легко вычислить с помощью так называемой схемы Горнера: $P_n(x) = (ldots ((a_nx + a_{n-1})x + a_{n-2})x + ldots)x + a_0$

Для вычисления значения во внутренних скобках a_nx + a_n-1 требуется одна операция умножения и одна операция сложения. Для вычисления значения в следующих скобках (a_nx + a_n-1)x + a_n-2 требуется опять одна операция умножения и одна операция сложения, т.к. a_nx + a_n-1 уже вычислено, и т.д.

Тогда в этом способе вычисление P_n(x) потребует n операций умножения и n операций сложения, т.е. сложность вычислений Op₂ = n+n = 2n.
Ясно, что Op₂ << Op₁.

Анализ метода

Сложность вычислений

Оценка сложности интерполирования функции складывается из количества операций для решения системы линейных алгебраических уравнений (СЛАУ) и нахождения значения полинома в точке.

Сложность решения СЛАУ, например, методом Гаусса для матрицы размера nxn: 2n³/3, т.е. O(n³).

Для нахождения полинома в заданной точке требуется n умножений и n сложений. В результате сложность метода: O(n³).

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

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

Ошибка приближения функции f(x) интерполяционным полиномом n-ой степени P_n(x) в точке x определяется разностью: R_n(x) = f(x) – P_n(x).

Погрешность R_n(x) определяется следующим соотношением:

$R_n(x) = frac{f^{n+1}(xi)}{(n+1)!}omega_n(x)$

Здесь $f^{n+1}(xi)$ – производная (n+1)-го порядка функции f(x) в некоторой точке а функция определяется как

Если максимальное значение производной fⁿ⁺¹(x) равно $M_{n+1} = sup_{x in [x_0, x_n]} left| f^{n+1}(x) right| ,$ то для погрешности интерполяции следует оценка: $left| R_n(x) right| = frac{M_{n+1}}{(n+1)!}omega_n(x).$

При реализации данного метода на ЭВМ ошибкой интерполяции E_n(x) будем считать максимальное уклонение полинома от исходной функции на выбранном промежутке: $E_n(x) = max_{x in [a, b]} left| f(x) - P_n(x) right| .$

Выбор узлов интерполяции

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

Введём следующее определение: полиномом Чебышева называется многочлен вида

T_k(x) = cos(k arccos x), |x|≤1.

Известно (см. ссылки литературы), что если узлы интерполяции x₀, x₁,…,x_n являются корнями полинома Чебышева степени n+1, то величина $max_{x in [a,b]} left| omega_{n+1}(x) right|$ принимает наименьшее возможное значение по сравнению с любым другим выбором набора узлов интерполяции.

Очевидно, что в случае k = 1 функция T₁(x), действительно, является полиномом первой степени, поскольку T₁(x) = cos(arccos x) = x.

В случае k = 2 T₂(x) тоже полином второй степени. Это нетрудно проверить. Воспользуемся известным тригонометрическим тождеством: cos2θ = 2cos²θ – 1, положив θ = arccos x.

Тогда получим следующее соотношение: T₂(x) = 2x² – 1.

С помощью тригонометрического тождества cos(k + 1)θ = 2cosθcoskθ – cos(k – 1) легко показать, что для полиномов Чебышева справедливо реккурентное соотношение:

T_k+1(x) = 2xT_k(x) – T_k-1(x)

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

Корни полинома Чебышева легко находятя из уравнения: T_k(x) = cos(k arccos x) = 0. Получаем, что уравнение имеет k различных корней, расположенных на отрезке [-1,1]: $x_m = cos frac{2m+1}{2k}pi, , m = 0,1, cdots, k-1,$ которые и следует выбирать в качестве узлов интерполирования.

Нетрудно видеть, что корни на [-1,1] расположены симметрично и неравномерно – чем ближе к краям отрезка, тем корни расположены плотнее. Максимальное значение модуля полинома Чебышева равно 1 и достигается в точках $cos frac{m}{k}pi.$

Если положить $omega_k(x) = frac1{2^{k-1}}T_k(x),$ то для того, чтобы коэффициент при старшей степени полинома ω_k(x) был равен 1, $max_{x in [-1,1]}omega_k(x) = frac1{2^{k-1}}.$

Известно, что для любого полинома p_k(x) степени k с коэффициентом, равным единице при старшей производной верно неравенство $max_{x in [-1,1]}p_k(x) geq frac1{2^{k-1}},$ т.е. полиномы Чебышева являются полиномами, наименее уклоняющимися от нуля.

Вычислительный эксперимент

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

Далее строится СЛАУ, которая решается методом Гаусса. На выходе получаем коэффициенты для полинома и ошибку аппроксимации.

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

Пример: Интерполяция синуса

Попробуем интерполировать функцию y = sin(x) на отрезке [1, 8.5]. Выберем узлы интерполяции: {1.1, 2, 4.7, 7.5, 8.5}

Полученный в результате интерполяции полином отображён на рисунке (синим цветом показан график y = sin(x), красным – интерполяционного полинома)

Ошибка интерполяции в этом случае: 0.1534

Давайте посмотрим, что произойдёт, если выбрать равномерно стоящие узлы {2, 3.5 5, 6.5, 8} для той же функции на том же отрезке.

На отрезке [3, 6] приближение, бесспорно, стало лучше. Однако разброс на краях очень большой. Ошибка интерполяции: 2.3466.

Наконец, выберем узлы интерполяции в соответствии с Чебышевским алгоритмом. Получим их по следующей формуле (просто сделаем замену переменной):
$x_m = frac{a+b}{2} + frac{b-a}{2}cos frac{pi(2m-1)}{2n}, , m=1,2, cdots, n.$

В нашем случае [a,b] – отрезок [1, 8.5], y = cosx, n+1 – количество узлов.

Остаётся открытым вопрос, какое количество узлов выбрать.

При значении n меньше 3 ошибка аппроксимации получается более 10.6626.
При n = 4: приближение лучше (ошибка равна 1.0111),
при n = 5: ошибка аппроксимации 0.2797

График функций при n = 4 выглядит следующим образом:

Продолжим исследования.

n = 6: ошибка аппроксимации 1.0233.

При n = 7 ошибка аппроксимации принимает наименьшее из полученных ранее значений (для данного промежутка): 0.0181. График синуса (обозначен синим цветом) и аппроксимационного полинома (обозначен красным цветом) представлены на следующем графике:

Что интересно, если при этом же количестве узлов выбирать их на отрезке [1, 8], то ошибка аппроксимации становится ещё меньше : 0.0124. График в этом случае выглядит так:

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

Ошибка аппроксимации будет только расти с увеличением числа узлов.

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

Вывод

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

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

Прикреплённые файлы

Программная реализация на языке С++ (исходный текст) [2.55Кб]

Литература

Н.С. Бахвалов, Н.П. Жидков, Г.М. Кобельков. Численные методы. Изд-во “Лаборатория базовых знаний”. Москва. 2003.
И.С. Березин, Н.П. Жидков. Методы вычислений. Изд. ФизМатЛит. Москва. 1962.
Дж. Форсайт, М.Мальком, К. Моулер. Машинные методы математических вычислений. Изд-во “Мир”. Москва. 1980.

Смотри также

Проблема выбора узлов для интерполяции
Ошибки вычислений
Метод наименьших квадратов
Интерполяционный многочлен Лагранжа
Практикум ММП ВМК, 4й курс, осень 2008

Ссылки

Собрание библиотек Boost

Источник

Методы функциональной интерполяции

Постановка задачи

Пусть на множестве [math]Omega=[a,b][/math] задана сетка [math]Omega_n= {x_i,, i=overline{0,n}}[/math], определяемая [math]n+1[/math] точкой [math]x_0,x_1,ldots,x_n[/math], а на сетке задана сеточная функция [math]y_i=f(x_i),~ i=overline{0,n}[/math]

[math]y_0=f(x_0),~ y_1=f(x_1),~ ldots,~ y_n=f(x_n).qquad mathsf{(4.6)}[/math]

В некоторых случаях [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,~ i=overline{0,n}[/math].

Точки [math]x_0,x_1,ldots, x_n[/math] называются узлами интерполяции, а искомая функция [math]y=F(x)[/math] — интерполирующей.

Геометрически это означает, что нужно найти кривую, проходящую через заданное множество точек [math](x_i,y_i),~ i=overline{0,n}[/math] (рис. 4.2).

Одной из целей задачи интерполяции является вычисление значения функции в произвольной точке [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_{i+k}] и рассмотрим сеточную функцию [math]y_i=f(x_i)[/math], заданную в (k+1)-м узле [math]x_i,x_{i+1},ldots,x_{i+k}[/math] (узлы не совпадают):

[math]y_i=f(x_i),~ y_{i+1}=f(x_{i+1}),~ ldots,~ y_{i+k}= f(x_{i+k}).[/math]

(4.7)

В качестве интерполирующей функции выберем алгебраический многочлен k-й степени (степень многочлена на единицу меньше количества узлов):

[math]F(x)= widetilde{f}_k(x,overline{a})= sumlimits_{j=0}^{k} a_jx^j= a_0+a_1x+ ldots+a_kx^k.[/math]

(4.8)

Будем искать неизвестные коэффициенты [math]a_0,a_1,ldots,a_k[/math] из условия интерполяции (4.3) , т.е. [math]delta widetilde{f}_k(x_i,overline{a})=0colon[/math]

[math]widetilde{f}_k(x_i,overline{a})=y_i,~ widetilde{f}_k (x_{i+1}, overline{a})= y_{i+1},~ ldots,~ widetilde{f}_k(x_{i+k}, overline{a})= y_{i+k}.[/math]

(4.9)

Теорема 4.1 (о единственности решения задачи интерполяции). Задача о нахождении интерполяционного многочлена [math]textstyle{widetilde{f}_k(x,overline{a})= sumlimits_{j=0}^{k} a_jx^j}[/math], удовлетворяющего условиям (4.9), на частичном отрезке [math][x_i,x_{i+k}][/math] по заданной сеточной функции (4.7) имеет единственное решение.

Доказательство. Запишем условия интерполяции (4.9) с учетом (4.8) и обозначения [math]y_i= f(x_i)= f_icolon[/math]

[math]begin{aligned}& a_0+a_1x_i+ a_2x_i^2+ ldots+ a_kx_i^k=f_i,\ & a_0+a_1x_{i+1}+ a_2x_{i+1}^2+ ldots+ a_kx_{i+1}^k=f_{i+1},\ & quadvdots\ & a_0+a_1x_{i+k}+ a_2x_{i+k}^2+ ldots+ a_kx_{i+k}^k=f_{i+k}. end{aligned}[/math]

(4.10)

Эта система линейных алгебраических уравнений относительно коэффициентов [math]a_0,a_1, ldots,a_k[/math] имеет единственное решение, так как определитель матрицы системы

[math]C_k(x_i,x_{i+1},ldots,x_{i+k})= begin{pmatrix}1& x_i& x_i^2& cdots& x_i^k\ 1& x_{i+1}& x_{i+1}^2& cdots& x_{i+1}^k\ vdots& vdots& vdots& ddots& vdots\ 1& x_{i+k}& x_{i+k}^2& cdots& x_{i+k}^kend{pmatrix}[/math]

(4.11)

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

Полагая [math]i=0,~ k=n[/math], приходим к глобальному способу решения поставленной задачи. А именно, если задана сеточная функция в (n+1)-м узле [math]x_0,x_1,ldots,x_n[/math] и требуется найти алгебраический многочлен с использованием условий интерполяции, то единственным решением задачи интерполяции является интерполяционный многочлен n-й степени:

[math]F(x)= widetilde{f}_n(x,overline{a})= sumlimits_{j=0}^{n} a_jx^j= a_0+a_1x+ ldots+ a_nx^n,[/math]

(4.12)

коэффициенты которого находятся из системы (4.10) при [math]i=0,~ k=n[/math].

Замечания

1. Матрица (4.11) имеет дискретный (точечный) характер, так как ее элементы вычисляются по дискретным значениям [math]x_i,x_{i+1},ldots,x_{i+k}[/math].

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_ncolon, y_i=f(x_i),, i=overline{0,n}[/math], где [math]x_iin [a,b]=[x_0,x_n][/math] — в общем случае неравноотстоящие узлы, определяемые шагами [math]h_{i+1}= x_{i+1}-x_i~ (h_{i+1}=text{var})[/math].

Воспользуемся сначала кусочным способом. Здесь и далее будем использовать обозначение [math]f_i=f(x_i)[/math].

Выделим “окно” или частичный отрезок [math][x_i,x_{i+1}][/math], содержащий только две точки (шаблон [math](x_i,x_{i+1})[/math]). Тогда многочлен Лагранжа, интерполирующий исходную функцию на данном шаблоне, имеет вид (где [math]P_{1i}(x),,P_{1i+1}(x)[/math] — коэффициенты)

[math]L_1(x)= frac{(x-x_{i+1})}{(x_{i}-x_{i+1})}f_i+ frac{(x-x_{i})}{x_{i+1}-x_i}f_{i+1}= P_{1i}(x)cdot f_i+ P_{1i+1}(x)cdot f_{i+1}.[/math]

(4.13)

Действительно, легко убедиться в том, что [math]L_1(x)[/math] — алгебраический многочлен первой степени, который удовлетворяет условиям функциональной интерполяции (4.3) , т.е. [math]L_1(x_i)=f_i,~ L_1(x_{i+1})= f_{i+1}[/math].

Выделим “окно” в виде двойного частичного отрезка [math][x_i,x_{i+2}][/math] c шаблоном [math](x_i,x_{i+1},x_{i+2})[/math]. Тогда многочлен Лагранжа записывается в виде (где [math]P_{2i}(x),, P_{2i+1}(x),, P_{2i+2}(x)[/math] коэффициенты)

[math]begin{aligned}L_2(x)&= frac{(x-x_{i+1})(x-x_{i+2})}{(x_{i}-x_{i+1})(x_{i}-x_{i+2})}f_i+ frac{(x-x_{i})(x-x_{i+2})}{(x_{i+1}-x_{i})(x_{i+1}-x_{i+2})}f_{i+1}+ frac{(x-x_{i})(x-x_{i+1})}{(x_{i+2}-x_{i})(x_{i+2}-x_{i+1})}f_{i+2}=\ &= P_{2i}(x)cdot f_{i}+ P_{2i+1}(x)cdot f_{i+1}+ P_{2i+2}(x)cdot f_{i+2}.end{aligned}[/math]

(4.14)

Легко проверить, что (4.14) — многочлен второй степени и также удовлетворяет условиям функциональной интерполяции:

[math]L_2(x_i)= f_i;qquad L_2(x_{i+1})=f_{i+1},qquad L_2(x_{i+2})=f_{i+2}.[/math]

Обобщив запись многочлена на “окно” для к -кратного частичного отрезка [math][x_i,x_{i+k}][/math] с шаблоном [math](x_i,x_{i+1},ldots,x_{i+k})[/math], можно записать многочлен Лагранжа в виде

[math]L_k(x)= sumlimits_{m=i}^{i+k} frac{(x-x_i)(x-x_{i+1})ldots (x-x_{m-1})(x-x_{m+1})ldots (x-x_k)}{(x_m-x_i)(x_m-x_{i+1})ldots (x_m-x_{m-1})(x_m-x_{m+1})ldots (x_m-x_k)}f_m= sumlimits_{m=i}^{i+k} P_{km}(x)cdot f_m.[/math]

(4.15)

где [math]P_{km}(x)[/math] — коэффициенты Лагранжа, которые для внутренних точек шаблона записываются следующим образом:

[math]P_{km}(x)= frac{(x-x_i)(x-x_{i+1})ldots (x-x_{m-1})(x-x_{m+1})ldots (x-x_k)}{(x_m-x_i)(x_m-x_{i+1})ldots (x_m-x_{m-1})(x_m-x_{m+1})ldots (x_m-x_k)},.[/math]

(4.16)

Легко проверить, что [math]P_{km}(x)[/math] удовлетворяют условию

[math]P_{km}(x_j)= begin{cases}1,& j=m,\ 0,& jne m,end{cases} i leqslant j leqslant i+k.[/math]

(4.16)

Если положить [math]i=0,, k=n[/math], то приходим к глобальному способу решения задачи. Тогда интерполяционный многочлен Лагранжа n-й степени имеет вид

[math]L_n(x)= sumlimits_{i=0}^{n} frac{(x-x_0)(x-x_{1})ldots (x-x_{i-1})(x-x_{i+1})ldots (x-x_n)}{(x_i-x_0)(x_i-x_{1})ldots (x_i-x_{i-1})(x_i-x_{i+1})ldots (x_i-x_n)}f_i= sumlimits_{i=0}^{n} P_{ni}(x)cdot f_i,[/math]

(4.17)

где коэффициенты Лагранжа [math]P_{ni}(x)[/math] во внутренних точках отрезка записываются в форме

[math]P_{ni}(x)= frac{(x-x_0)(x-x_{1})ldots (x-x_{i-1})(x-x_{i+1})ldots (x-x_n)}{(x_i-x_0)(x_i-x_{1})ldots (x_i-x_{i-1})(x_i-x_{i+1})ldots (x_i-x_n)},.[/math]

Очевидно, многочлен [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.

[math]begin{array}{|c|c|c|c|c|c|c|} multicolumn{7}{r}{mathit{Table~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{matrix}{}\[-8pt]{}end{matrix}Pi_{n+1}(x)= (x-x_0)(x-x_1)(x-x_2)ldots (x-x_n)}& D_i&f_i \hline end{array}[/math]

Здесь [math]D_i[/math] — произведение элементов i-й строки, [math]Pi_{n+1}(x)[/math] — произведение элементов главной диагонали, [math]y_i= f(x_i)= f_i,~ i=overline{0,n}[/math]. Тогда многочлен Лагранжа может быть записан в форме

[math]L_n(x)= Pi_{n+1}(x)cdot sumlimits_{i=0}^{n} frac{f_i}{D_i},.[/math]

(4.18)

Замечания

1. Если заданная сеточная функция такая, что [math]f_i=1~(i=overline{0,n})[/math], то из (4.17) следует, что [math]L_n(x)equiv1[/math], и поэтому справедливо равенство [math]textstyle{sumlimits_{i=0}^{n} P_{ni}(x)=1}[/math] (сумма всех коэффициентов Лагранжа в точке [math]x[/math] равна нулю), которое можно использовать для контроля правильности расчетов.

2. Коэффициенты Лагранжа [math]P_{ni}(x)[/math] для некоторой функции [math]y_i= f(x_i)[/math] определяются лишь узлами сетки и точкой [math]x[/math], в которой необходимо вычислить значение многочлена. Если в некоторой точке [math]x_{ast}[/math] требуется определить значения нескольких интерполируемых функций [math]f_1(x_i),f_2(x_i),ldots~ (i=overline{1,n})[/math], то коэффициенты Лагранжа для всех исходных функций подсчитываются только один раз.

3. При введении дополнительных узлов интерполяции все коэффициенты многочлена Лагранжа необходимо пересчитывать заново, что неудобно на практике. От этого недостатка свободны многочлены Ньютона.

Перейдем к рассмотрению примеров решения задачи интерполяции на основе вышеизложенной методики.

▼ Пример 4.1

Построить многочлен Лагранжа третьей степени для сеточной функции, заданной табл. 4.2. Вычислить значение функции в точке [math]x_{ast}=2,!5[/math]. Для записи многочлена использовать формулу (4.18).

[math]begin{array}{|c|c|c|c|c|} multicolumn{5}{r}{mathit{Table~4.2}}\hline i& 0& 1& 2&3 \hline x_i& 2& 3& 4&5 \hline begin{matrix}{}\[-8pt]{}end{matrix} f(x_i)= f_i & 7& 5& 8&7 \hline end{array}[/math]

Решение. 1. Построим многочлен Лагранжа. Для этого составим табл. 4.3, соответствующую табл. 4.1.

[math]begin{array}{|c|c|c|c|c|c|} multicolumn{6}{r}{mathit{Table~4.3}}\hline x-2&-1&-2&-3&-6cdot(x-2)&7 \hline 1&x-3&-1&-2& 2cdot(x-3)&5 \hline 2& 1& x-4&-1&-2cdot(x-4)&8 \hline 3& 2& 1& x-5& 6cdot(x-5)&7 \hline multicolumn{4}{|c|}{begin{matrix}{}\[-8pt]{}end{matrix} Pi_4(x)= (x-2)(x-3)(x-4)(x-5)}& D_i&f_i \hline end{array}[/math]

По формуле (4.18) получаем:

[math]begin{aligned}L_3(x) &= Pi_4(x) sumlimits_{i=0}^{3} frac{f_i}{D_i}=\ &=-frac{7}{6}(x-3)(x-4)(x-5)+ frac{5}{2}(x-2)(x-4)(x-5)-4(x-2)(x-3)(x-5),+\ &quad +,frac{7}{6} (x-2)(x-3)(x-4)=-frac{3}{2},x^3+16x^2-frac{107}{2},x+62. end{aligned}[/math]

2. Вычислим значение функции в заданной точке: [math]L_3(2,!5)= 4,!8125[/math].

Погрешность интерполяции многочленами Лагранжа

При определении значения [math]f(x),~ xne x_{i}[/math] для функции [math]y_i=f(x_i)~ (i=overline{0,n})[/math] с помощью многочлена Лагранжа возникает погрешность или остаточное слагаемое [math]R_n(x)colon[/math]

[math]f(x)= L_n(x)+ R_n(x).[/math]

(4.19)

Здесь предполагается, что используется глобальный способ интерполяции и что [math]f(x)in C_{n+1}[a,b][/math]. Последнее предположение требуется для применения соответствующих теорем математического анализа, однако, приведенное ниже соотношение для [math]R_n(x)[/math] может использоваться и для сеточных функций.

На основе указанных предположений доказано, что при интерполяции функции [math]y_i=f(x_i)[/math], заданной в общем случае на неравномерной сетке [math]Omega_n[/math], интерполяционным многочленом Лагранжа [math]textstyle{L_n(x)= sumlimits_{i=0}^{n} P_{ni}(x)f_i}[/math] для произвольного значения [math]xin (x_0,x_n)[/math] возникает погрешность

[math]R_n(x)= frac{f^{(n+1)}(xi)}{(n+1)!}cdot omega_n(x),[/math]

(4.20)

где [math]omega_n(x)= (x-x_0)(x-x_1)ldots (x-x_n)[/math] — многочлен (n+1)-й степени, а [math]xiin (a,b)[/math]. Поскольку точно найти [math]R_n(x)[/math] нельзя (из-за неопределенности точки [math]xi[/math]), то при проведении вычислений обычно находятся только приближенные оценки погрешностей интерполяции, которые являются априорными.

Оценка погрешности интерполяции в некоторой произвольной фиксированной точке [math]x_{ast}in[a,b][/math] имеет вид, где [math]M_{n+1}= max_{xin[a,b]}|f^{(n+1)}(x)|[/math]

[math]bigl|f(x_{ast}-L_n(x_{ast})bigr| leqslant frac{M_{n+1}}{(n+1)!}cdot bigl|omega_n(x_{ast})bigr|,[/math]

(4.21)

Оценка максимальной погрешности интерполяции в любой точке [math]xin[a,b][/math], т.е. на всем отрезке [math][a,b]colon[/math]

[math]bigl|f(x-L_n(x)bigr| leqslant frac{M_{n+1}}{(n+1)!}cdot max_{xin[a,b]}bigl|omega_n(x_{ast})bigr|.[/math]

(4.22)

Замечание. Для сеточных функций с фиксированными узлами сетки (узлами интерполяции) также можно проводить оценки погрешности по формулам (4.21), (4.22), однако для этого необходимо численно определять [math]M_{n+1}[/math] с помощью аппарата численного дифференцирования. При этом следует учитывать, что при вычислении производных высокого порядка возникают большие погрешности.

▼ Пример 4.2

С какой точностью можно вычислить значение [math]f(x)= Bigl.{sqrt{x} }Bigr|_{x_{ast}=85}[/math], если вычисления производить на основе интерполяционного многочлена Лагранжа первой и второй степени, а в качестве сеточной функции принять [math](x_i,y_i)= bigl{(16;4),(36;6),(100;10)bigr},~ (i=0,1,2)[/math]?

Решение. Так как требуется вычислить погрешность только в одной точке [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_2colon[/math]

[math]f'(x)= frac{1}{2}x^{-1!!not{phantom{|}},2};qquad f”(x)=-frac{1}{4}x^{-3!!not{phantom{|}},2};qquad M_2= max_{xin[36;100]} bigl|f”(x)bigr|= frac{1}{4 sqrt{36^3}}= frac{1}{864},.[/math]

Тогда [math]omega_1(x_{ast})= (85-36)(85-100)=-735;~~ |R_1(x)| leqslant frac{M_2}{(n+1)!}|omega_1(x_{ast})|= frac{735}{2cdot864}= 0,!425[/math].

Найдем оценку погрешности в точке [math]x_{ast}[/math] для [math]L_2(x)[/math]. В качестве “окна” квадратичной интерполяции [math](n=2)[/math] выбираем отрезок [math][x_0,x_2]= [16;100][/math]. Определим величину [math]M_3colon[/math]

[math]f”'(x)= frac{3}{8}x^{-5!!not{phantom{|}},2};qquad M_3= max_{xin[16;100]} frac{3}{8 sqrt{16^5}}= frac{1}{2730,!7},.[/math]

Тогда [math]omega_2(x_{ast})= (85-16)(85-36)(85-100)=-50715;~~ |R_2(x)| leqslant frac{M_3}{(n+1)!}|omega_2(x_{ast})|=3,!1[/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_{n+1}[a,b][/math].

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

[math]Omega_1= bigl{x_0^{(1)},x_1^{(1)}bigr},quad Omega_2= bigl{x_0^{(2)}, x_1^{(2)},x_2^{(2)}bigr},quad ldots,quad Omega_n= bigl{x_0^{(n)},x_1^{(n)},ldots ,x_n^{(n)}bigr}[/math]

с возрастающим разбиением [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]limlimits_{ktoinfty} L_k(x_{ast})= f(x_{ast})[/math] (поточечная сходимость).

Для отрезка [math][a,b][/math] существует понятие равномерной сходимости в некоторой норме, например max, [math]max_{substack{xin[a,b]\ ntoinfty} bigl|f(x)-L_n(x)bigr|to0}[/math].

Характер сходимости или расходимости интерполяционного процесса зависит как от гладкости и поведения функции [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|,~ xin[-1;1][/math] и [math]f(x)= (1+25x^2)^{-1},~ xin[-1;1][/math]. Кроме того, применение многочленов высоких степеней приводит к так называемым “провалам” между узлами интерполяции, часто называемым осцилляциями. Указанные свойства интерполяционных процессов обусловливают нецелесообразность применения интерполяционных многочленов высоких степеней. В связи с этим в вычислительной практике для сеточных функций степень [math]n[/math] не берут выше [math]5div 8~ (nleqslant 5div 8)[/math] и задание частичного отрезка согласуют с выбранной степенью многочлена.

Рассмотрим часто использующиеся на практике линейную и параболическую интерполяцию.

Линейная и параболическая интерполяция с помощью многочлена Лагранжа

В прикладных расчетах часто применяется простейшая кусочная интерполяция, основанная на многочленах первой степени [math]L_1(x)[/math] или второй степени [math]L_2(x)[/math]. В этом случае функциональная интерполяция называется линейной или параболической (квадратичной) соответственно. Рассмотрим возможные способы их реализации и найдем оценки их погрешностей.

Пусть для сеточной функции [math]y_i=f(x_i)[/math], заданной на сетке [math]Omega_n= {x_0,x_1, ldots, x_n}[/math], требуется выполнить интерполяционный процесс для определения значения [math]f(x_{ast})[/math], где [math]x_{ast}ne x_i~ (i=overline{0,n})[/math], и оценить погрешности. Для обоснованного выбора степени интерполяционного многочлена необходимо указать, какую погрешность имеют значения исходной функции [math]f(x_i)[/math] в узлах. Если эта погрешность составляет величину [math]O(h_{i+1}^2)[/math] или [math]O(h_{i+1}^3[/math], а в широком классе вычислительных задач обеспечиваются именно такие погрешности, следует использовать линейную или параболическую интерполяцию.

Методика решения задачи линейной интерполяции

1. По расположению заданной точки [math]x_{ast}[/math] на оси [math]Ox[/math] выбрать из всех частичных отрезков [math][x_0,x_1], [x_1,x_2],ldots, [x_i,x_{i+1}], ldots, [x_{n-1},x_n][/math], заданных своими крайними значениями и в совокупности образующих сетку [math]Omega_n[/math], “окно” интерполяции [math]Oequiv [x_i,x_{i+1}][/math], такое, что [math]x_i< x_{ast}< x_{i+1}[/math].

2. Для отрезка [math][x_i,x_{i+1}][/math] вычислить значения коэффициентов Лагранжа [math]P_{1i}(x_{ast})= frac{x_{ast}-x_{i+1}}{x_i-x_{i+1}}[/math] и [math]P_{1i+1}(x_{ast})= frac{x_{ast}-x_i}{x_{i+1}-x_i}[/math] входящих в формулу (4.13) для многочлена [math]L_1(x)[/math]. Правильность полученных значений [math]P_{1i}(x_{ast}),, P_{1i+1} (x_{ast})[/math] проверить по условию [math]P_{1i}(x_{ast})+ P_{1i+1}(x_{ast})=1[/math], которое должно выполняться.

3. Вычислить искомое значение [math]f(x_{ast})[/math] согласно (4.13):

[math]f(x_{ast})approx L_1(x_{ast})= P_{1i}(x_{ast}) f_i+ P_{1i+1} (x_{ast}) f_{i+1}.[/math]

Как показано ниже, порядок этой аппроксимации равен двум, т.е. [math]O(h_{i+1}^2)[/math].

Геометрическая интерпретация линейной интерполяции при известной формульной функции [math]y=f(x)[/math] (штриховая линия) изображена на рис. 4.3. Здесь прямая [math]AB[/math] соответствует графику функции [math]y=L_1(x)[/math] на отрезке [math][x_i, x_{i+1}][/math]. Приближенное значение функции равно [math]L_1(x_{ast})[/math] ( точка [math]C[/math]), и оно отстоит от точного значения [math]f(x_{ast})[/math] на величину [math]CD approx O(h_{i+1}^2)[/math] вдоль оси [math]Oy[/math].

Методика решения задачи параболической интерполяции

1. Из всей совокупности спаренных частичных отрезков [math][x_0,x_2], [x_1,x_3], ldots, [x_{i-1}, x_{i+1}], [x_i, x_{i+2}], ldots, [x_{n-1},x_n][/math], образующих сетку [math]Omega_n[/math], по заданной величине [math]x_{ast}[/math] выбрать два пересекающихся “окна” [math]O_1equiv [x_{i-1},x_{i+1}],~ O_2equiv [x_i, x_{i+2}][/math] (предполагается, что [math]x_{ast}[/math] принадлежит внутреннему отрезку [math][x_i,x_{i+1}][/math]).

2. Для “окон” [math]O_1[/math] и [math]O_2[/math] вычислить значения коэффициентов Лагранжа: [math]P_{2,i-1}^{(1)}(x_{ast}),, P_{2,i}^{(1)}(x_{ast}),, P_{2,i+1}^{(1)}(x_{ast})[/math] — для “окна” [math]O_1[/math] и [math]P_{2,i}^{(2)}(x_{ast}),, P_{2,i+1}^{(2)}(x_{ast}),, P_{2,i+2}^{(2)}(x_{ast})[/math] — для “окна” [math]O_2[/math]. Путем суммирования проверить правильность полученных значений коэффициентов:

[math]sumlimits_{k=i-1}^{i+1} P_{2,k}(x_{ast})=1,qquad sumlimits_{k=i}^{i+2} P_{2,k}(x_{ast})=1.[/math]

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]f(x_{ast})approx frac{1}{2}bigl[L_2^{(1)}(x_{ast})+ L_2^{(2)}(x_{ast})bigr]equiv L_{2text{sr}}(x_{ast}).[/math]

При этом порядок интерполяции повышается на единицу, т.е. [math]L_{2text{sr}}(x_{ast})[/math] аппроксимирует точное значение с четвертым порядком, поскольку погрешность составляет величину [math]O(H^4)[/math], где [math]H=max{h_i,h_{i+1},h_{i+2}}[/math].

Замечание. Если [math]x_{ast}in [x_0,x_2][/math] или [math]x_{ast}in [x_{n-1},x_n][/math], то выбирается одно “окно” [math][x_0,x_2][/math] или [math][x_{n-1},x_n][/math] соответственно.

Геометрическая интерпретация параболической интерполяции изображена на рис. 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_{2text{sr}} (x_{ast})[/math].

Приведем оценки погрешностей линейной и параболической интерполяции. Вначале предположим, что сетка [math]Omega_n[/math] равномерная (это имеет значение только для параболической интерполяции). Формулы (4.13), (4.14) и оценки (4.22), записанные для “окон” интерполяции, упрощаются, если ввести в рассмотрение новую переменную — фазу интерполяции и [math]u=frac{x-x_i}{h}colon[/math]

[math]begin{gathered}L_1(u)= (1-u)cdot f_i+ ucdot f_{i+1},\ L_2(u)= frac{(u-1)(u-2)}{2}cdot f_i-u(u-2)cdot f_{i+1}+ frac{u(u-1)}{2}cdot f_{i+2}. end{gathered}[/math]

Здесь учтено, что [math]x-x_{i+1}= x-(x_i+h)= uh-h= h(u-1),~ x-x_{i+2}= h(u-2)[/math]. Очевидно, что для [math]L_1(u)[/math] величина [math]{u}[/math] изменяется в диапазоне [math]0 leqslant u leqslant 1[/math], а для [math]L_2(u)[/math] — в диапазоне [math]0 leqslant u leqslant 2[/math]. Для получения мажорант в оценке (4.22) необходимо найти [math]max|omega_1(x)|[/math] и [math]max|omega_2(x)|[/math]. Преобразуя зависимости [math]omega_1(x)[/math] и [math]omega_2(x)[/math] к новой переменной [math]{u}[/math], получаем [math]omega_1(u)= h^2ucdot (1-u),[/math] [math]omega_2(u)= h^3ucdot (1-u)(2-u)[/math] и находим максимумы:

[math]Omega^1= max_{0 leqslant u leqslant 1} bigl|omega_1(u)bigr|= frac{h^2}{4},,qquad Omega^2= max_{0 leqslant u leqslant 2} bigl|omega_2(u)bigr|= frac{2 sqrt{3}h^2}{9},.[/math]

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

[math]bigl|f(u)-L_1(u)bigr|leqslant frac{h^2}{8},M_{2,i},[/math]

(4.23)

[math]bigl|f(u)-L_2(u)bigr|leqslant frac{sqrt{3}h^3}{27},M_{3,i},[/math]

(4.24)

где [math]M_{2,i}= max_{[x_i, x_{i+1}]}bigl|f”(x)bigr|,~ M_{3,i}= max_{[x_i, x_{i+2}]}bigl|f”'(x)bigr|[/math]. При этом предполагается, что [math]f(x)[/math] принадлежит классам функций [math]f(x)in C_2[a,b],~ f(x)in C_3[a,b][/math] соответственно для линейной и параболической интерполяции. Таким образом, из оценок (4.23), (4.24) следует, что линейная интерполяция обеспечивает на частичном отрезке [math][x_i, x_{i+1}][/math] второй порядок аппроксимации или погрешности по [math]h[/math], а параболическая (без осреднения) на двойном отрезке [math][x_i, x_{i+2}][/math] — третий порядок. Данные пофешности, как отмечалось во введении и предыдущих разделах, сокращенно записываются как [math]O(h^2)[/math] и [math]O(h^3)[/math]. При реализации алгоритма с осреднением порядок параболической интерполяции становится равным [math]O(h^4)[/math].

Замечания

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] приведена оценка

[math]max_{[x_i,x_{i+2}]} bigl|R_2(x)bigr| leqslant frac{sqrt{3}}{27},H^3M_{3,i},[/math]

(4.25)

где [math]H= max(h_{i+1}, h_{i+2}),~ h_{i+2}= x_{i+2}-x_{i+1},~ M_{3,i}= max_{[x_i,x_{i+2}]} bigl|f”'(x)bigr|[/math].

2. Если гладкость функции [math]f(x)[/math] не достигает вышеуказанной и класс гладкости понижен на единицу [math]bigl(f(x)in C_2[a,b]bigr)[/math], то порядок параболической интерполяции также понижается на единицу:

[math]max_{[x_i,x_{i+2}]} bigl|R_2(x)bigr| leqslant 0,!1546cdot H^2M_{2,i}.[/math]

(4.26)

3. Повышение класса гладкости функции [math]f(x)[/math] выше [math]C_3[a,b][/math] не приводит к увеличению порядка параболической интерполяции относительно Н , т.е. происходит как бы его “замораживание” или “насыщение”. Указанные свойства, вытекающие из оценок (4.25), (4.26), носят общий характер и имеют место в других оценках аппроксимации многочленов, сплайнов, производных и интегралов.

4. Для произвольной степени интерполяционного многочлена при [math]h=text{const},~ f(x)in C_{n-1}[a,b][/math] погрешность функциональной интерполяции на отрезке [math][x_0,x_n][/math] выражается следующим образом:

[math]bigl|f(x)-L_n(x)bigr|_{[a,b]}= max_{xin [a,b]} bigl|f(x)-L_n(x)bigr| leqslant h^{n+1} frac{M_{n+1}}{(n+1)!},Omega^n.[/math]

(4.27)

где [math]Omega^n= max_{0 leqslant u leqslant n} bigl|omega_n(u)bigr|,~ M_{n+1}= max_{[x_0,x_n]} bigl|f^{(n+1)}(x)bigr|,~ omega_n(x)= u(u-1)cdot ldotscdot (u-n)[/math]. Для многочленов с [math]n leqslant 5[/math] величины [math]Omega^n[/math] или их оценки являются такими:

[math]Omega^1=frac{1}{4};quad Omega^2=frac{2 sqrt{3}}{9};quad Omega^3=1;quad Omega^4<3,!7;quad Omega^5<17.[/math]

Из (4.27) следует, что на отрезке [math][x_0,x_n][/math] величина [math]|R_n(x)|[/math] есть [math]O(h^{n+1})[/math] и при уменьшении шага [math]h[/math] в два раза мажоранта уменьшается по крайней мере в [math]2^{n+1}[/math] раз. Отсюда вытекает, что для непрерывных функций можно по заданной точности интерполяции выбирать шаг [math]h[/math]. При этом можно путем использования кусочной интерполяции в некоторых пределах изменять степень интерполяционного многочлена. Мажоранту в оценке (4.27) при кусочной интерполяции можно также снизить путем выбора “окна” интерполяции [math][x_i,x_{i+k}][/math] так, чтобы точка [math]x_{ast}[/math] располагалась как можно ближе к его середине. Это обусловлено тем, что колебания функции [math]omega_n(x)[/math] вблизи середины [math][x_i,x_{i+k}][/math] меньше, чем у его концов.

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

Перейдем к рассмотрению примеров решения задач линейной и параболической интерполяции.

▼ Пример 4.3

Дана сеточная функция, являющаяся сеточным представлением формульной функции [math]f(x)=x^3[/math] (табл. 4.4). Найти значение [math]f(x_{ast})[/math] при [math]x_{ast}=2[/math] с помощью линейной и параболической интерполяции.

[math]begin{array}{|c|c|c|c|c|c|} multicolumn{6}{r}{mathit{Table~4.4}}\hline i& 0& 1& 2& 3& 4 \hline x_i& -1& 0& 1& 3& 4 \hline f_i& -1& 0& 1& 27& 64 \hline end{array}[/math]

Решение. Применение линейной интерполяции. Воспользуемся методикой.

1. Выбираем “окно” интерполяции [math][x_i,x_{i+1}]= [1;3]~ (i=2,~ i+1=3)[/math].

2. Вычислим коэффициенты Лагранжа: [math]P_{1,2}= frac{x_{ast}-x_3}{x_2-x_3}= 0,!5;~~ P_{1,3}= frac{x_{ast}-x_2}{x_3-x_2}=0,!5[/math]. Вычисления выполнены верно, так как [math]P_{1,2}+P_{1,3}=1[/math].

3. Определим искомое значение [math]f(2)[/math] по формуле (4.13):

[math]f(2)approx L_1(2)= P_{1,2}cdot f_2+ P_{1,3}cdot f_3=14.[/math]

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

1. Выбираем “окна” интерполяции [math]O_1equiv [x_1;x_3]= [0;3];~ O_2equiv [x_2;x_4]= [1;4][/math].

2. Вычислим коэффициенты Лагранжа [math]P_{21}^{(1)},, P_{22}^{(1)},, P_{23}^{(1)}[/math] для “окна” [math]O_1[/math] и [math]P_{22}^{(2)},, P_{23}^{(2)},, P_{24}^{(2)}[/math] для “окна” [math]O_2[/math]. Получим

[math]P_{21}^{(1)}=-frac{1}{3},quad P_{22}^{(1)}=1,quad P_{23}^{(1)}=frac{1}{3};qquad P_{22}^{(2)}=frac{1}{4},quad P_{23}^{(2)}=1,quad P_{24}^{(2)}=-frac{1}{3},.[/math]

Вычисления выполнены правильно, так как [math]textstyle{sumlimits_{k=1}^{3} P_{2k}^{(1)}=1;~ sumlimits_{k=2}^{4} P_{2k}^{(2)}=1}[/math].

3. Определим значения [math]L_2^{(1)}(x_{ast}),~ L_2^{(2)}(x_{ast})[/math] по формуле (4.14):

[math]L_2^{(1)}(x_{ast})=-frac{1}{3}cdot0+1cdo1+ frac{1}{3}cdot27=10;qquad L_2^{(2)}(x_{ast})= frac{1}{3}cdot1+ 1cdot27-frac{1}{3}cdot64=6.[/math]

Искомое значение получим в результате осреднения: [math]f(2)approx frac{L_2^{(1)}(x_{ast})+ L_2^{(2)}(x_{ast})}{2}=8[/math]. Найдено точное значение [math]f(2)=8[/math], что объясняется тем, что оператор осреднения повышает порядок интерполяции. В этом случае в остаточное слагаемое входит [math]f^{(4)}(xi)equiv0[/math] и оно равно нулю.

Многочлены Ньютона

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

Разделенные разности вводятся для функции [math]y_i= f(x_i)=f_i,~ i=overline{0,n}[/math], заданной на неравномерной сетке [math](h_{i+1}= text{var})[/math], а конечные разности — для функции [math]y_i= f(x_i)=f_i,~ i=overline{0,n}[/math], определенной на равномерной сетке [math](h_{i+1}= text{const})[/math].

Выбрав внутри неравномерной или равномерной сетки соответствующие шаблоны интерполяции [math](x_i,x_{i+1}), (x_i,x_{i+1},x_{i+2}), ldots, (x_i,x_{i+1},ldots, x_{i+k})[/math] введем следующие определения разделенных и конечных разностей:

– разделенная разность первого порядка: [math]f(x_i,x_{i+1})= frac{f_{i+1}-f_i}{x_{i+1}-x_i}[/math];

– разделенная разность второго порядка: [math]f(x_i,x_{i+1},x_{i+2})= frac{f(x_{i+1},x_{i+2})-f(x_i,x_{i+1})}{x_{i+2}-x_i}[/math];

– разделенная разность k-го порядка: [math]f(x_i,x_{i+1},ldots, x_{i+k})= frac{f(x_{i+1},x_{i+2},ldots, x_{i+k})-f(x_i,x_{i+1},ldots, x_{i+k-1})}{x_{i-k}-x_i}[/math];

– конечная разность первого порядка: [math]Delta f_i= f_{i+1}-f_i[/math];

– конечная разность второго порядка: [math]Delta^2f_i= Delta (Delta f_i)= Delta f_{i+1}-Delta f_i= f_{i+2}-2f_{i+1}+f_i[/math];

– конечная разность k-го порядка: [math]textstyle{Delta^kf_i= Delta (Delta^{k-1}f_i)= sumlimits_{j=0}^{k} (-1)^jC_k^j f_{i+j}}[/math], где [math]C_k^j= frac{k!}{(k-j)!,j!}[/math].

Последовательность получения разделенных и конечных разностей при [math]k=3[/math] для произвольной функции наглядно представляют табл. 4.5 и 4.6.

Для гладких функций числовые значения [math]f(x_j, x_{j+1}, ldots, x_{j+k})[/math] и [math]Delta^kf_i[/math] при возрастании [math]k[/math] уменьшаются и стремятся к нулю, т.е. [math]f(x_j, x_{j+1}, ldots, x_{j+k})to0[/math] и [math]Delta^kf_ito0[/math] при [math]ktoinfty[/math].

Связь между разделенными и конечными разностями k-го порядка при [math]h=text{const}[/math] устанавливается следующим соотношением:

[math]f(x_i,x_{i+1},ldots,x_{i+k})= frac{Delta^kf_i}{k!,h^k},.[/math]

(4.28)

Действительно, при [math]k=1[/math] для разделенной разности [math]f(x_i,x_{i+1})[/math] получаем [math]f(x_i,x_{i+1})= frac{Delta f_i}{h}= frac{Delta f_i}{1!,h^1}[/math], то есть (4.28) при [math]k=1[/math] справедливо.

Пусть [math]k=2[/math]. Тогда получаем

[math]f(x_i,x_{i+1},x_{i+2})= frac{f(x_{i+1},x_{i+2})-f(x_{i},x_{i+1})}{}= frac{1}{2h}! left(frac{f_{i+2}-f_{i+1}}{h}-frac{f_{i+1}-f_i}{h}right)= frac{Delta f_{i+1}-Delta f_i}{2h^2}= frac{Delta^2f_i}{2h^2},.[/math]

Таким образом, связь (4.28) выполняется и при [math]k=2[/math]. Справедливость (4.28) при произвольном к можно доказать методом математической индукции.

Интерполяционный многочлен Ньютона для неравномерной сетки

Пусть исходная (интерполируемая) сеточная функция [math]y_i=f(x_i),~ i=overline{0,n}[/math], задана на неравномерной сетке [math]Omega_nequiv {x_0,x_1,x_2,ldots,x_n}[/math], характеризующейся шагами [math]h_{i+1}= x_{i+1}-x_i=text{var}[/math].

Воспользуемся сначала кусочным способом. Из всей совокупности узлов выбираем шаблон [math](x_i,x_{i+1},ldots,x_{i+k})[/math], соответствующий некоторому “окну” интерполяции [math][x_i,x_{i+k}][/math].

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

[math]begin{aligned}N_k(x)= f_i&+ f(x_i,x_{i+1})cdot (x-x_i)+ f(x_i,x_{i+1}, x_{i+2})cdot (x-x_i)cdot (x-x_{i+1})+ ldots+\ &+, f(x_i,x_{i+1},ldots,x_{i+k})cdot (x-x_i)cdot (x-x_{i+1})cdot ldotscdot (x-x_{i+k-1}). end{aligned}[/math]

(4.29)

Действительно, [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_{i+1})[/math]) и [math]k=2[/math] (шаблон [math](x_i,x_{i+1},x_{i+2}})[/math]). Пусть [math]k=1[/math]. Тогда

[math]N_1(x)= f_i+f(x_i,x_{i+1})cdot (x-x_i),[/math]

(4.30)

и поэтому

[math]N_1(x_i)= f_i,;qquad N_1(x_{i+1})= f_i+ frac{f_{i+1}-f_i}{x_{i+1}-x_i}(x_{i+1}-x_i)= f_{i+1}.[/math]

Таким образом, условия интерполяции для [math]N_1(x)[/math] выполнены, следовательно, многочлен (4.30) может быть использован для линейной интерполяции кусочным способом. Пусть [math]k=2[/math]. Тогда

[math]N_2(x)= f_i+ f(x_i,x_{i+1})(x-x_i)+ f(x_i,x_{i+1},x_{i+2}) (x-x_i)(x-x_{i+1}),[/math]

(4.31)

и поэтому

[math]begin{gathered}N_2(x)=f_i;qquad N_2(x_{i+1})= f_i+ frac{f_{i+1}-f_i}{x_{i+1}-x_i}(x_{i+1}-x_i)+0= f_{i+1},,\ N_2(x_{i+2})= f_i+ frac{f_{i+1}-f_i}{x_{i+1}-x_i}(x_{i+2}-x_i)+ frac{1}{x_{i+2}-x_i}! left(frac{f_{i+2}-f_{i+1}}{x_{i+2}-x_{i+1}}-frac{f_{i+1}-f_{i}}{x_{i+1}-x_{i}}right)!(x_{i+2}-x_i)(x_{i+2}-x_{i+1})=f_{i+2}. end{gathered}[/math]

Таким образом, условия интерполяции для многочлена [math]N_2(x)[/math] также выполнены и он может использоваться для параболической интерполяции кусочным способом.

Для произвольного [math]k[/math] справедливость равенств [math]N_k(x_j)=f_i,~ j=overline{i,i+k}[/math], проверяется методом математической индукции.

Полагая [math]i=0,~ k=n[/math], приходим к глобальному способу. Тогда интерполяционный многочлен Ньютона n-й степени имеет вид

[math]begin{aligned}N_n(x)= f_0&+ f(x_0,x_1)(x-x_0)+ f(x_0,x_1,x_2) (x-x_0)(x-x_1)+ ldots+\ &+,f(x_0,x_1,ldots,x_n)(x-x_0)(x-x_1)cdot ldots+(x-x_{n-1}).end{aligned}[/math]

(4.32)

Замечания

1. Согласно теореме 4.1 многочлен Ньютона (4.32) является тождественным многочлену [math]textstyle{widetilde{f}_n(x,overline{a})= sumlimits_{j=0}^{n} a_jx^j}[/math] с коэффициентами, получаемыми из системы (4.10), либо многочлену Лагранжа, т.е. [math]L_n(x)= N_n(x)= widetilde{f}_n(x,overline{a})[/math], если узлы интерполяции и интерполируемая функция одинаковы.

2. Интерполяционный многочлен Ньютона (4.29) или (4.32) (так же, как и многочлен Ньютона, выражаемый ниже через конечные разности) записан не через значения функции, как это имеет место для многочлена Лагранжа, а через разделенные разности. Поэтому при изменении степени [math]k[/math] в процессе интерполирования у многочлена Ньютона [math]N_k(x)[/math] требуется только добавить или отбросить соответствующее число слагаемых. Это иногда упрощает алгоритм интерполирования.

3. При интерполяции на основе (4.29) или (4.32) узлы интерполяции [math]x_i,x_{i+1}, ldots, x_{i+k}[/math] или [math]x_0,x_1,ldots,x_n[/math] определяющие шаблоны интерполяции, целесообразно выбирать так, чтобы точка [math]x_{ast}[/math] была расположена возможно ближе к середине отрезка [math][x_i,x_{i+k}][/math] или [math][x_0,x_n][/math].

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_{i+1}=text{const}[/math] для всех [math]i[/math], то многочлен (4.29), соответствующий шаблону [math](x_i,x_{i+1},ldots, x_{i+k})[/math], путем подстановки в него вместо разделенных разностей их выражений через конечные разности, согласно (4.28), преобразуется к виду

[math]N_{k}^{(I)}(q)= f_i+ frac{Delta f_i}{1!},q+ frac{Delta^2 f_i}{2!},q(q-1)+ ldots+ frac{Delta^k f_i}{k!}q(q-1)cdot ldotscdot (q-k+1),[/math]

(4.33)

где [math]q=frac{x-x_i}{h}[/math] — фаза интерполяции, определенная относительно точки [math]x_i[/math]; [math]Delta^jf_i~(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_{i+1},ldots, x_{i+k})[/math], когда [math]q>0[/math]. Если [math]x_i=x_0[/math], а [math]x_{ast}<x_0[/math], то этот же многочлен используется и для экстраполяции левее точки [math]x_0~(q<0)[/math]. Поэтому многочлен [math]N_k^{(I)}(q)[/math] называется интерполяционным многочленом для интерполяции вперед (в начале таблицы ) или для экстраполяции назад. В (4.33) этот многочлен обозначен цифрой I, указанной в скобках вверху. Для определенности назовем его первым интерполяционным многочленом Ньютона.

Полагая [math]i=0,~k=n[/math], получаем решение задачи глобальной интерполяции на всем отрезке [math][x_0,x_n]colon[/math]

[math]N_{n}^{(I)}(q)= f_0+ frac{Delta f_0}{1!},q+ frac{Delta^2 f_0}{2!},q(q-1)+ ldots+ frac{Delta^n f_i}{n!}q(q-1)cdot ldotscdot (q-n+1),[/math]

(4.34)

где [math]q=frac{x-x_0}{h}[/math]. Остаточное слагаемое этого многочлена имеет вид

[math]R_n(q)= h^{n+1}cdot frac{q(q-1)cdotldotscdot(q-n)}{(n+1)!}cdot f^{(n+1)}(xi),[/math]

где [math]xiin(a,b)[/math] — некоторое промежуточное значение между узлами [math]x_0,x_1,ldots,x_n[/math] и точкой [math]x[/math].

Если фазу интерполяции определить относительно [math]x_{i+k}[/math] некоторой конечной точки шаблона [math](x_i,x_{i+1},ldots, x_{i+k})[/math], то есть [math]widehat{q}= frac{x-x_{i+k}}{h}[/math], то вместо (4.33) получается второй интерполяционный многочлен Ньютона:

[math]N_{k}^{(II)}(widehat{q})= f_{i+k}+ frac{Delta f_{i+k-1}}{1!},widehat{q}+ frac{Delta^2 f_{i+k-2}}{2!},widehat{q}(widehat{q}+1)+ ldots+ frac{Delta^k f_i}{k!}widehat{q}(widehat{q}+1)cdot ldotscdot (widehat{q}+k-1),[/math]

(4.35)

Данный многочлен удобно применять в конце выделенного шаблона или всей таблицы [math]y_i= f(x_i),~ i=overline{0,n}[/math].

Если [math]i+k=n[/math], а [math]x_{ast}>x_n[/math], то (4.35) используется для экстраполяции правее точки [math]x_n~(widehat{h}>0)[/math]. Поэтому многочлен [math]N_{k}^{(II)}(widehat{q})[/math] называется многочленом для интерполяции назад (в конце таблицы) или экстраполяции вперед.

Полагая [math]i=0,~ i+k=n[/math], получаем решение задачи глобальной интерполяции — второй интерполяционный многочлен Ньютона n-й степени (где [math]widehat{q}= frac{x-x_n}{h}[/math]):

[math]N_{n}^{(II)}(widehat{q})= f_{n}+ frac{Delta f_{n-1}}{1!},widehat{q}+ frac{Delta^2 f_{n-2}}{2!},widehat{q}(widehat{q}+1)+ ldots+ frac{Delta^n f_0}{n!}widehat{q}(widehat{q}+1)cdot ldotscdot (widehat{q}+n-1),[/math]

(4.36)

Остаточное слагаемое многочлена (4.36) имеет вид

[math]R_n(widehat{q})= h^{n+1}cdot frac{widehat{q}(widehat{q}+1)cdot ldotscdot (widehat{q}+n)}{(n+1)!}cdot f^{(n+1)}(xi)[/math], где [math]xiin(a,b)[/math].

Схема выбора узлов интерполяции при изменении степеней интерполяционных многочленов [math]N_{k}^{(I)}(q),~ N_{k}^{(II)}(widehat{q})~(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_{n-1},x_n][/math], то — второй (рис. 4.5,а). Если точка [math]x_{ast}[/math] находится вдали от концов отрезка [math][x_0,x_n][/math], то может применяться как первый интерполяционный многочлен, так и второй (рис. 4.5,б).

Замечания

1. Из формул интерполяционных многочленов [math]N_{k}^{(I)}(q)[/math] и [math]N_{k}^{(II)}(widehat{q})[/math] видно, что повышение их степеней в процессе реализации алгоритма не требует пересчета предыдущих слагаемых, входящих в многочлены с меньшими степенями.

2. Для гладких функций при повышении порядка конечных разностей справедливо свойство: [math]Delta^kf_ito0[/math] при [math]ktoinfty[/math], и поэтому, как только очередное слагаемое в рассматриваемых многочленах становится меньше требуемой точности интерполяции, увеличение степени [math]N_{k}^{(I)}(q),~ N_{k}^{(II)}(widehat{q})[/math] следует прекратить. Это замечание справедливо также и для [math]N_k(x)[/math].

▼ Пример 4.4-4.6

Пример 4.4. Построить многочлен Ньютона третьей степени для сеточной функции, заданной табл. 4.7. Вычислить значение функции в точке [math]x_{ast}=2,!5[/math]. Решить задачу интерполяции при включении одного дополнительного значения сеточной функции: [math]f(1)=5[/math].

[math]begin{array}{|c|c|c|c|c|} multicolumn{5}{r}{mathit{Table~4.7}}\hline i& 0& 1& 2& 3\hline x_i& 2&3& 4& 5\hline f(x_i)=f_i& 7&5& 8& 7\hline end{array}[/math]

Решение

Построим многочлен Ньютона (4.32), справедливый для произвольного расположения узлов. Для этого составим табл. 4.8, аналогичную табл. 4.5:

[math]begin{gathered}f(x_0,x_1)= frac{5-7}{3-2}=-2;qquad f(x_1,x_2)= frac{8-5}{4-3}=3;qquad f(x_2,x_3)= frac{7-8}{5-4}=-1;\[2pt] f(x_0,x_1,x_2)= frac{3-(-2)}{4-2}=frac{5}{2};qquad f(x_1,x_2,x_3)= frac{-1-3}{5-3}=-2;qquad f(x_0,x_1,x_2,x_3)= frac{-2-frac{5}{2}}{5-2}=-frac{3}{2},.end{gathered}[/math]

По формуле (4.32) для [math]n=3[/math] имеем

[math]begin{gathered}N_3(x)= f(x_0)+ (x-x_0)f(x_0,x_1)+ (x-x_0)(x-x_1)f(x_0,x_1,x_2)+\ +(x-x_0)(x-x_1)(x-x_2)f(x_0,x_1,x_2,x_3)= 7+(x-2)cdot(-2)+ (x-2)(x-3)cdot frac{5}{2}+\ (x-2)(x-3)(x-4)! left(-frac{3}{2}right)=-frac{3}{2},x^3+ 16x^2-frac{107}{2},x+62. end{gathered}[/math]

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

[math]N_{3}^{(I)}(q)= f(x_0)+ frac{Delta f_0}{1!},q+ frac{Delta^2f_0}{2!},q(q-1)+ frac{Delta^3f_0}{3!},q(q-1)(q-2),[/math]

где [math]q=frac{x-x_0}{h}= frac{x-2}{1}=x-2[/math]. Составим табл. 4.9, аналогичную табл. 4.6. Имеем:.

[math]begin{gathered}Delta f_0= 5-7=-2;quad Delta f_1=8-5=3;quad Delta f_2=7-8=-1;quad Delta^2f_0=3-(-2)=5;\[2pt] Delta^2f_1=-1-3=-4;qquad Delta^3f_0=-4-5=-9. end{gathered}[/math]

Поэтому

[math]begin{aligned}N_{3}^{(I)}(q)&= f(x_0)+ frac{(-2)}{1!},q+ frac{5}{2!},q(q-1)+ frac{(-9)}{3!},q(q-1)(q-2)= left.{left(-frac{3}{2},q^3+7q^2-frac{15}{2},q+7right)}right|_{q=frac{x-2}{1}}\ &=-frac{3}{2}cdot (x^3-6x^2+12x-8)+ 7cdot (x^2-4x+4)-frac{15}{2}cdot (x-2)+7=-frac{3}{2},x^3+ 16x^2-frac{107}{2},x+62. end{aligned}[/math]

Запишем также второй интерполяционный многочлен Ньютона (4.36):

[math]begin{aligned}N_{3}^{(II)}(q)&= left.{left(f_3+ frac{Delta f_2}{1!},widehat{q}+ frac{Delta^2f_1}{2!}, widehat{q} (widehat{q}+1)+ frac{Delta^3f_0}{3!}, widehat{q}(widehat{q}+1)(widehat{q}+2)right)!}right|_{widehat{q}=frac{x-x_0}{1}}=\ &=7+ frac{-1}{1!},widehat{q}+ frac{-4}{2!}, widehat{q} (widehat{q}+1)+ frac{-9}{3!}, widehat{q}(widehat{q}+1)(widehat{q}+2)=\ &= left.{left(-frac{3}{2},widehat{q}^3-frac{13}{2},widehat{q}^2-6widehat{q}+7right)!}right|_{widehat{q}=x-5}=-frac{3}{2},x^3+ 16x^2-frac{107}{2},x+62. end{aligned}[/math]

Сравнивая с результатом примера 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]begin{aligned}N_4(x)&= N_3(x)+(x-x_0)(x-x_1)(x-x_3)(x-x_4)f(x_0,x_1,x_2, x_3,x_4)=\[2pt] &=-frac{3}{2},x^3+ 16x^2-frac{107}{2},x+62+ (x-2)(x-3)(x-4)(x-5)! left(-frac{3}{4}right)=\[2pt] &=-frac{3}{2},x^3+ 16x^2-frac{107}{2},x+62-frac{3}{4}cdot (x^4-14x^3+71x^2-154x+120)=\[2pt] &=-frac{3}{4},x^4+ 9x^3-frac{149}{4},x^2+62x-28. end{aligned}[/math]

Легко проверить, что новый узел можно добавить и в начале таблицы, т.е. будет найден тот же многочлен [math]N_4(x)[/math]. На рис. 4.6 изображены полученные интерполяционные многочлены.

Вычислим значение функции в точке [math]x_{ast}=2,!5colon[/math]

[math]N_3^{(I)}(x_{ast})= N_3^{(II)}(x_{ast})= 4,!8125;qquad N_4^{(I)}(x_{ast})= 5,!5156.[/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_{i+1}),~ f(x_i,x_{i+1},x_{i+2})[/math] (табл. 4.11)

Здесь [math]f(x_1,x_2)= frac{27-1}{3-1}=13,~ f(x_3,x_4)= frac{64-27}{4-3}=37,~ f(x_2,x_3,x_4)= frac{37-13}{4-1}=8.[/math]. Определим [math]N_1(x),~N_2(x)colon[/math]

[math]begin{aligned}N_1(x)&= f_2+f(x_2,x_3)cdot (x-x_2)= 1+13cdot(x-1)=13x-12.\[2pt] N_2(x)&= f_2+ f(x_2,x_3)(x-x_2)+ f(x_2,x_3,x_4)(x-x_2)(x-x_3)=\ &=13x-12+8(x-1)(x-3)= 8x^2-19x+12. end{aligned}[/math]

Сравнивая [math]N_2(x)[/math] с соответствующим многочленом Лагранжа L_2(x) (см. пример 4.3), видим, что они совпадают. Это подтверждает единственность решения задачи о построении интерполяционного многочлена.

2. Вычислим значения [math]Bigl.{N_1(x)}Bigr|_{x=x_{ast}},~ Bigl.{N_2(x)}Bigr|_{x= x_{ast}}colon[/math]

[math]N_1(2,!5)= 13cdot2,!5-12=20,qquad N_2(2,!5)= 8cdot2,!5^2-19cdot 2,!5+12=14,!5.[/math]

Итак, линейная интерполяция приводит к отличию от точного значения на [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;qquad x_{ast 2}=4,!5;qquad x_{ast 3}=3,!5;qquad x_{ast 4}=5,!5.[/math]

Решение

Для подсчета значения функции в точке [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]begin{gathered}N_1^{(I)}(q)= f_0+ frac{Delta f_0}{1!},q= 7+frac{-2}{1},q=7-2q,\[2pt] N_2^{(II)}(q)= N_1^{(I)}(q)+ frac{Delta^2f_0}{2!},q(q-1)= 7-2q+ frac{5}{2},q(q-1)= frac{5}{2},q^2-frac{9}{2},q+7. end{gathered}[/math]

где [math]q=frac{x-2}{1}[/math] (см. табл. 4.9). В результате подстановки получаем

[math]N_1^{(I)}(x)=-2x+11;qquad N_2^{(I)}(x)= frac{5}{2},x^2-frac{29}{2},x+26.[/math]

В заданной точке [math]N_1^{(I)}(2,!5)=6;~ N_2^{(I)}(2,!5)=5,!375[/math].

Для подсчета значения в точке [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]begin{gathered}N_1^{(II)}(widehat{q})= f_3+ frac{Delta f_2}{1!},widehat{q}= 7+frac{-1}{1},widehat{q}= 7-widehat{q},\[2pt] N_2^{(II)}(widehat{q})= N_1^{(II)}(widehat{q})+ frac{Delta^2f_1}{2!}, widehat{q}(widehat{q}+1)= 7-widehat{q}+ frac{-4}{2},widehat{q} (widehat{q}+1)=-2widehat{q}^2-3widehat{q}+7, end{gathered}[/math]

где [math]widehat{q}=frac{x-x_3}{h}=x-5[/math]. В результате подстановки получаем

[math]N_1^{(II)}(x)=-x+12;qquad N_2^{(II)}(x)=-2x^2+17x-28.[/math]

В заданной точке [math]N_1^{(II)}(4,!5)=7,!5;~ N_2^{(II)}(4,!5)=8[/math].

Многочлен [math]N_1^{(II)}(x)[/math] можно использовать для экстраполяции, т.е. для подсчета функции в точке [math]x_{ast 4}=5,!5colon, 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]begin{gathered}N_1^{(I)}(q)= f_1+frac{Delta f_1}{1!},q= 5+frac{3}{1},q=5+3q,\[2pt] N_2^{(I)}(q)= N_1^{(I)}(q)+ frac{Delta^2f_1}{2!},q(q-1)= 5+3q+frac{-4}{2},q(q-1)=-2q^2+5q+5, end{gathered}[/math]

где [math]q=frac{x-x_1}{h}=x-3[/math]. После подстановки получаем

[math]N_1^{(I)}(x)= 3x-4;~~ N_2^{(I)}(x)=-2x^2+17x-28[/math] и [math]N_1^{(I)}(3,!5)= 6,!5;~~ N_2^{(I)}(3,!5)=7[/math].

Подсчитать значения в точке [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,~ k=1)[/math]. По формуле (4.35) получаем

[math]begin{gathered}N_1^{(II)}(widehat{q})= f_2+ frac{Delta f_1}{1!},widehat{q}= 8+frac{3}{1},widehat{q},\[2pt] N_2^{(II)}(widehat{q})= N_1^{(II)}(widehat{q})+ frac{Delta^2 f_0}{2!},widehat{q}= 8+frac{3}{1},widehat{q}(widehat{q}+1)= frac{5}{2},widehat{q}^2+ frac{11}{2},widehat{q}+8. end{gathered}[/math]

где [math]widehat{q}=frac{x-x_2}{h}=x-4[/math]. После подстановки

[math]N_1^{(II)}(x)=3x-4;~~ N_2^{(II)}(x)= frac{5}{2},x^2-frac{29}{2},x+26[/math] и [math]N_1^{(II)}(3,!5)=6,!5;~~ N_2^{(II)}(3,!5)=5,!875[/math].

Математический форум (помощь с решением задач, обсуждение вопросов по математике).

Кнопка "Поделиться"

Если заметили ошибку, опечатку или есть предложения, напишите в комментариях.

Источник

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

Полиномы Чебышева. Интерполяция по Чебышевским узлам.

Материал из MachineLearning.

Содержание

Полином в каноническом виде

Способ вычисления полинома в точке

Анализ метода

Сложность вычислений

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

Выбор узлов интерполяции

Вычислительный эксперимент

Пример: Интерполяция синуса

Рекомендации программисту

Предварительные настойки

Использование программы

Вывод

Прикреплённые файлы

Литература

Смотри также

Ссылки

Методы функциональной интерполяции

Постановка задачи

Многочлен Лагранжа

Погрешность интерполяции многочленами Лагранжа

Сходимости функционального интерполяционного процесса для непрерывных функций

Линейная и параболическая интерполяция с помощью многочлена Лагранжа

Методика решения задачи линейной интерполяции

Методика решения задачи параболической интерполяции

Многочлены Ньютона

Интерполяционный многочлен Ньютона для неравномерной сетки

Интерполяционные многочлены Ньютона для равномерной сетки

Вам также может понравиться

Как в ida найти адрес

Как найти всех культистов глаза космоса

Как найти характеристическое уравнение для системы

Добавить комментарий Отменить ответ