Интерполяция функций: методы приближения

5 Методы приближения функций. Интерполяция табличных
функций.
5.1 Постановка задачи приближения функций.
Аппроксимация функций заключается в приближенной замене заданной
функции f(x) некоторой функцией φ(x) так, чтобы отклонение функции φ(x) от f(x)
в заданной области было наименьшим. Функция φ(х) при этом называется
аппроксимирующей.
Необходимость аппроксимации функций в основном связана с двумя
причинами:
1) Функция f(x) имеет сложное аналитическое описание, вызывающее
определенные
трудности
при
его
использовании
(например,
является
спецфункцией: гамма-функцией, эллиптической функцией и др.).
2) Аналитическое описание функции f(x) неизвестно, т. е. f(x) задана
таблично. При этом необходимо иметь аналитическое описание, приближенно
представляющее f(x) (например, для вычисления значений f(x) в произвольных
точках, определения интегралов и производных от f(x) и т. п.)
Задачи аппроксимации можно разделить на 2 вида:
1) Интерполяция
2) Наилучшее приближение
Также
выделяют
экстраполяцию
являющуюся
частным
случаем
интерполяции или наилучшего приближения, при котором ищется приближенное
значение табличной функции вне заданного интервала.
Рассмотрим аппроксимацию табличной функции:
xi
x0
x1
x2
...
xn
f(xi)
y0
y1
y2
...
yn
Точки xi данной сеточной функции называют узлами аппроксимацию, а их
совокупность – аппроксимационной сеткой.
Пары (xi, yi) называют точками данных или базовыми точками.
Шаг сеточной функции может быть как переменным, так и постоянным.
Если функция F(x), удовлетворяет условию:
F( x 0)= y 0, F (x 1)= y 1, ... , F (x n )= y n
(5.1)
1
То
такая
функция
называется
интерполирующей
функцией
или
интерполянтом, а процесс называется интерполяцией.
Интерполирующая функция должна принадлежащую известному классу.
Наиболее часто в качестве интерполирующей функции используют различного
вида полиномы.
Интерполяцию можно разделить на глобальную – непрерывная функция для
всех точек xi. Кусочная (или локальная) – кусочно-непрерывная функция для всех
точек xi, при этом несколько соседних узлов интерполируются непрерывной
функцией.
Методы глобальной интерполяции обычно применяют для функций,
заданных небольшим количеством точек, т. к. при увеличении количества точек
увеличивается порядок интерполирующего многочлена, что отрицательно
сказывается на гладкости получаемой функции. Для кусочной интерполяции
количество узлов сетки большого значения не имеет.
Наибольшее распространение из методов глобальной интерполяции имеют
полиномы Лагранжа, первая и вторая формулы Ньютона, формула Бесселя,
формула Стирлинга, первая и вторая формулы Гаусса. Для кусочного
интерполирования
чаще
всего
применяются
линейная,
квадратичная
интерполяции, и интерполяции кубическими и B-сплайнами.
5.2 Глобальная интерполяция.
5.2.1 Интерполяция полиномами Лагранжа.
Интерполяционный полином Лагранжа n-ой степени есть линейная
комбинация базисных полиномов Лагранжа:
n
Ln ( x ) =∑ c i⋅li ( x )
(5.2)
i=0
Для того, чтобы удовлетворить условиям интерполяции, в качестве
коэффициентов линейной комбинации можно взять значения табличной функции
в узлах сетки:
ci= yi
А
в
качестве
базовых
(5.3)
полиномов,
полиномы,
в
узлах
сетки
удовлетворяющие условию:
2
li ( x j )= 0, i≠ j
1, i= j
{
(5.4)
Первой части условия (5.4) удовлетворяют полиномы вида:
i−1
n
li ( x j )= A i⋅∏ ( x−x j )⋅ ∏ ( x−x j )
j=0
(5.5)
j=i+1
Из второй части имеем:
1
A i = i−1
n
(5.6)
( x i− x j )⋅ ∏ ( x i− x j )
∏
j=0
j=i+1
Отсюда, базисные полиномы Лагранжа имеют вид:
i−1
n
∏ ( x−x j )⋅ ∏ ( x− x j )
j=0
li ( x j )= i−1
j=i+1
n
(5.7)
∏ ( xi −x j )⋅ ∏ ( xi −x j )
j=0
j=i+1
А интерполяционный полином принимает значение:
i−1
n
∏ ( x−x j )⋅ ∏ ( x−x j )
j=0
Ln ( x ) =∑ y i⋅ i−1
i=0
(5.8)
j=i+1
n
(5.8)
∏ ( xi −x j )⋅ ∏ ( xi −x j )
j=0
Уравнение
n
называется
j=i+1
Лагранжевой
формой
записи
интерполирующего полинома.
Для примера, Лагранжева форма записи линии, проходящей через точки
(x0, y0) и (x1, y1):
L1 ( x ) = y 0⋅
x−x 1
x−x 0
+ y 1⋅
x 0−x 1
x 1−x 0
(5.9)
Параболы, проходящей через точки (x0, y0), (x1, y1) и (x2, y2):
( x−x 1 )⋅( x−x 2 )
( x−x 0 )⋅( x−x 2)
( x−x 0 )⋅( x− x 1)
+ y 1⋅
+ y 2⋅
( x 0 −x 1 )⋅( x 0−x 2 )
( x 1−x 0 )⋅( x 1−x 2 )
( x 2−x 0 )⋅( x 2−x 1 )
L2 ( x ) = y 0⋅
Уравнения (5.9)-(5.1) есть Лагранжева форма записи
(5.10)
линейной и
квадратичной интерполяции.
3
5.2.2 Интерполяционные формулы Ньютона.
Если сетка, на которой задана табличная функция имеет постоянный шаг,
для построения интерполирующего полинома можно воспользоваться введенным
до этого понятием конечных разностей.
Зададим интерполирующий полином в виде:
P n (x)=a 0+ a1⋅( x−x 0) +a 2⋅( x−x 0 )⋅( x−x 1 )+
...+a n⋅( x−x 0 )⋅( x− x 1)⋅...⋅( x−x n)=
n
i−1
(5.11)
=a 0+ ∑ a i⋅∏ ( x−x j )
i=1
j=0
Коэффициенты ai находятся из условий:
P n (x i )= y i
(5.12)
Так, для i = 0 имеем:
a 0= y 0
(5.13)
a 0+ a1⋅( x 1−x 0 )= y 0
(5.14)
y 1− y 0 Δyy 0
=
x 1− x 0
h
(5.15)
Для i = 1:
Откуда:
a 1=
Общая формула для коэффициентов ai будет иметь вид:
i
Δy y 0
( i)
y (x 0 )
a i=
=
i!
i ! hi
(5.16)
И интерполирующий полином преобразуется к виду:
n
( i)
y (x 0 ) i−1
P n (x)= y 0 +∑
⋅∏ ( x− x j )
i!
j=0
i=1
(5.17)
Схема (5.17) называется прямой (первой) интерполяционной формулой
Ньютона.
Заметим, что коэффициенты ai аналогичны коэффициентам формулы
Тейлора при разложении в точке x = x0, только вместо вместо производной
аналитической функции используется ее аппроксимация правой разностной
производной.
Так же как и формула Тейлора, полином (5.17) обладает наибольшей
точностью в окрестностях точки x0, при удалении от которой точность падает.
4
Если рассуждать аналогично, но в качестве базовой точки взять xn, получим
полином:
n
n−i+ 1
P n (x)= y n +∑ ai⋅ ∏ ( x− x j )
i=1
(5.18)
j=n
С коэффициентами:
i
Δy y n
( i)
y (x n )
a i=
=
i!
i ! hi
(5.19)
Только в отличии от (5.16) в (5.19) используются левые разностные
производные.
Полином (5.18) с коэффициентами (5.19) называется обратной (второй)
интерполяционной формулой Ньютона.
Обратная интерполяционная формула Ньютона обладает наибольшей
точностью в окрестностях точки xn.
Формулы Ньютона обладают большим удобством использования по
сравнению с формулой Лагранжа, т. к. степень полинома в них легко
регулируется количеством слагаемых в формулах (5.17) и (5.18).
5.2.3 Интерполяционные формулы Гаусса
Формулы
Ньютона,
использующие
правые
и
левые
разностные
производные наиболее хорошо интерполируют табличную функцию в начале (x0)
и конце таблицы (xn).
Для лучшей интерполяции центральной части таблицы используются
формулы, основанные на центральных разностных производных, каковыми
являются интерполяционные формулы Гаусса.
Базовыми точками для интерполяционных формул Гаусса являются точки,
близкие к центру таблицы.
В отличие от (5.11) в полиноме Гаусса при увеличении степени происходит
не последовательное подключение узлов (0, 1, 2… для первой формулы Ньютона
и n, n - 1, n - 2… для второй формулы Ньютона), а поочередное подключение
узлов как справа, так и слева от базовой точки. При этом, если очередность
подключения узлов сетки имеет вид: k, k + 1, k - 1, k + 2, k - 2… , т. е. сначала
происходит подключение правого узла, потом левого, интерполяционный
полином носит название прямой (первой) интерполяционной формулой Гаусса
5
если же сначала подключается левый узел, потом правый, то имеем обратную
(вторую) интерполяционную формулу Гаусса.
Базовую точку (xk) выбирают следующим образом, если в таблице нечетное
количество узлов, выбирается узел, симметричный, относительно границ
таблицы:
xk = xn/2. Если в таблице четное количество узлов, то для первой
интерполяционной формулы берется узел, ближе к началу таблицы: xk = (xn - 1)/2,
для второй формулы берется узел, ближе к концу таблицы: xk = (xn + 1)/2, так
делается для того, чтобы использовать все узлы сетки для получения
интерполирующего полинома.
Первая интерполяционная формула Гаусса имеет вид:
(2)
y (x k )
Gn (x)= y k + ( x−x k ) y (x k+1 /2)+( x−x k )⋅( x− x k+ 1)
+
2!
y (3) (x k+1 /2)
+( x−x k )⋅( x−x k+ 1)⋅( x− x k−1 )
+...
3!
( 1)
(5.20)
Или в обобщенном виде:
n
{
{
j=
i−1
2
y (i) (x k+ 1/2)
∏ ( x− xk+ j ) i ! , if i mod 2 ≠ 0
1−i
Gn (x)= y k +∑ j= 2
i=1
j=
i
2
∏ ( x−x k+ j )
j=
2−i
2
(5.21)
(i)
y (x k )
, if i mod 2 = 0
i!
Вторая интерполяционная формула Гаусса в общем виде:
n
j=
i−1
2
y (i) (x k−1 /2)
∏ ( x− xk+ j ) i ! , if i mod 2 ≠ 0
1−i
Gn (x)= y k +∑ j= 2
i=1
j=
i−2
2
∏ ( x−x k+ j )
j=−
i
2
(5.22)
(i)
y ( xk )
, if i mod 2 = 0
i!
Формулы (5.21) и (5.22) иногда еще называют формулами нижней и верхней
интерполяции, т. к. в них используются верхние и нижние центральные разности,
т. е. разности на полшага выше и ниже базовой точки. При четных значения
степени полинома используются разностные производные в базовой точке, при
6
нечетных, производные, в точке смещенной на полшага вправо или влево от
базовой точки.
5.2.4 Интерполяция по формуле Стирлинга.
На
интерполяционных
формулах
Гаусса
основываются
формулы,
позволяющие дать более симметричную оценку вблизи базовой точки. Одной из
них является формула Стирлинга, определенная как среднее арифметической
первой и второй формул Гаусса. Если табличная функция имеет нечетное
количество точек, то формула Стирлинга принимает вид:
n
S n (x)= y k + ∑
i=1
{
j=
i−1
2
y ( i) (x k−1 /2)+ y (i) (x k+1 /2)
, if i mod 2 ≠ 0
∏ ( x− xk+ j )
2i !
1−i
j=
(
x−x
2
2−i
k+
2
+ x−x
)(
i−2
k+
2
2
)
(5.23)
i
j= −1
2
y (i)(x k )
∏ ( x−xk+ j ) i ! , if i mod 2 = 0
i
j=1−
2
Интерполяционный полином Стирлинга точно центрирован на середину
таблицы и дает наилучший результат при x = xk ± h/4.
5.2.5 Интерполяция по формуле Бесселя.
На
основе
интерполяционных
формулах
Гаусса
основана
также
интерполяционная формула Бесселя, которая является средним арифметическим
второго интерполяционного полинома Гаусса в базовых точках xk и xk+1:
Интерполяционный полином Бесселя центрирован в точке xk+1/2 и дает
наилучший результат при x = xk + h/2 ± h/4.
5.3 Кусочная интерполяция.
5.3.1 Кусочно-линейная интерполяция.
Из
методов
кусочной
интерполяции
самым
простейшим
является
кусочно-линейная интерполяции, при которой все точки табличной функции
соединяются отрезками прямой (рис. 5.1).
7
f(x)
x2
x1
x0
x
Рисунок 5.1. Кусочно-линейная интерполяция.
Функция f(x) имеет вид:
f 1( x), x∈[x 0, x 1 ]
f 2( x), x∈[x 1, x 2 ]
f (x)= .............................
f i (x), x∈[ x i−1 , x i ]
.............................
f n (x), x∈[ x n−1 , x n]
{
(5.24)
Где функции fi(x) представляют собой отрезки прямой:
f i ( x)=a i⋅x+b i
(5.25)
Коэффициенты уравнений прямых ai и bi легко находятся из условий
прохождения прямой через точки xi-1 и xi. Для удовлетворения этих условий
коэффициенты ai и bi должны удовлетворять системам уравнений:
ai⋅x i−1 +b i= y i−1
ai⋅x i +bi = y i
{
(5.26)
Кусочно-линейная интерполяция является самой простой, и поэтому
довольно часто применяется для расчета значений между узлами интерполяции
«на коленке», для построения интерполирующей зависимости, используемой в
дальнейших научных и инженерных расчетах, обычно используются более
продвинутые методы интерполяции.
5.3.2 Кусочно-квадратичная интерполяция
Путем
повышения
степени
интерполирующего
полинома
получают
кусочно-квадратичную интерполяцию, при применении которой, кусочная
функция состоит из кусков парабол, проходящих через три точки (рис. 5.2).
8
f(x)
x0
x2
x1
x
Рисунок 5.2. Кусочно-квадратичная интерполяция.
Функция f(x) имеет вид:
f 1( x), x∈[ x 0, x 2 ]
f 2 (x), x∈[ x 2, x 4 ]
f (x)= .............................
f i ( x), x∈[x 2 i− 2 , x 2 i ]
.............................
f n (x), x∈[ x n−2 , x n]
{
(5.27)
2
Уравнения парабол:
f i ( x)=a i⋅x 2 +bi⋅x+c i
(5.28)
Коэффициенты парабол находятся аналогично предыдущему случаю, через
условия принадлежности узлов сетки параболе:
ai⋅x 22i−2 +bi⋅x 2 i−2 +c i = y 2 i−2
ai⋅x 22i−1 +bi⋅x 2 i−1 +c i = y 2i−1
ai⋅x 22i +b i⋅x 2 i +c i = y 2 i
{
(5.29)
Рассуждая аналогичным образом можно построить как кусочно-кубические
интерполяции так и интерполяции полиномами более высоких степеней, однако
на практике такое встречается редко, т. к. уменьшается гладкость получаемых
кусочный функций и увеличивается сложность получения коэффициентов
полинома.
5.3.3 Сплайн-интерполяция
Существенным
недостатком
кусочно
полиномиальной
интерполяции
является отсутствие гладкости в глобальном масштабе. Т.е., несмотря на то, что
каждый из отдельных «кусков» обладает непрерывной производной, производная
9
глобальной функции имеет разрывы в узлах перехода от одного «куска» к
другому.
Данный недостаток позволяет решить интерполяция сплайнами.
Введем некоторые обозначения:
l – гладкость сплайна, порядок высшей непрерывной производной сплайна;
m – порядок сплайна, показатель степени интерполирующего полинома;
d – дефект сплайна, разность между порядком сплайна и его гладкостью.
Следуя
данным
терминам,
интерполяционный
сплайн
это
l
раз
непрерывно дифференцируемый на отрезке x∈[ x 0, x n ] полином порядка m. При
этом узлы интерполяционного сплайна не всегда совпадают с узлами табличной
функции.
Таким образом, рассмотренная нами кусочно-линейная интерполяция
является сплайном порядка 1 и дефекта 1, а кусочно-квадратичная – сплайном
порядка 2 и дефекта 2. Все ранее рассмотренные кусочно-полиномиальные
интерполяции обладают нулевой гладкостью.
Попытаемся построить кубический сплайн дефекта 1, с узлами,
совпадающими с узлами табличной функции.
Данный сплайн иногда называют чертежным, или естественным сплайном.
Данное определение пошло из способа построения данного сплайна на бумаге:
бралась достаточно гибкая линейка и прикладывалась к чертежу так, чтобы
проходить через все точки, по профилю выгнутой линейки строился график.
Чертежный сплайн должен быть дважды непрерывно дифференцируем на отрезке
x∈[ x 0, x n ] .
Зададим такой сплайн уравнением:
g1 (x), x∈[ x 0, x 1]
g2 (x), x∈[ x 1, x 2 ]
f (x)= .............................
gi (x), x∈[ x i−1 , x i ]
.............................
g n (x), x∈[x n−1 , x n ]
{
(5.30)
Кусочные функции будут являться полиномами третьего порядка:
3
2
gi (x)=ai⋅( x−x i ) +b i⋅( x−x i ) +c i⋅( x−x i ) +d i
(5.31)
10
Таким
образом
в
функции
(5.30)
присутствует
4n
неизвестных
коэффициентов ai, bi, ci, di.
Функция (5.30) должна удовлетворять условиям интерполяции:
g1 (x 0)= y 0
(5.32)
gi (x i )= y i
(5.33)
Условиям гладкости:
gi−i (x i−1 )=g i ( x i−1 )
( 1)
g(1)
i−i (x i−1 )=g i (x i−1 ) , i = 2..n
( 2)
g(2)
i−i (x i−1 )=g i (x i−1 )
(5.34)
Формулы (5.32-5.34) содержат 4n – 2 уравнений, недостающие два
уравнения, позволяющие однозначно определить коэффициентны ai, bi, ci, di.
возьмем из допущения о отсутствии выпуклости/вогнутости сплайна на концах
отрезка (условие расслабленности концов сплайна):
(2)
g1 (x 0 )=0
(2)
g n (x n )=0
(5.35)
Таким образом система (5.33)-(5.35) однозначно разрешима относительно
коэффициентов ai, bi, ci, di.
Зададим шаг между узлами через:
h i= x i− x i−1
(5.36)
Подставляя функции (5.31) в уравнения (5.33)-(5.35) имеем:
−a1⋅h 31+b 1⋅h21−c 1⋅h1 +d 1= y 0
(5.37)
di = yi
(5.38)
d i−1=−a i h 3i +b i h2i −c i h i +d i
, i = 2..n
c i−1=3⋅ai h2i −2⋅b i h i +c i
b i−1=−3⋅a i⋅h i +b i
(5.39)
−3⋅a 1⋅h1 +b 1=0
(5.40)
b n=0
(5.41)
Система (5.37)-(5.41) представляет собой СЛАУ, которую можно решать
как прямыми, так и итерационными методами.
Задачи построения сплайнов другого прядка и дефектов решаются
аналогично.
11
Список литературы
1. Вержбицкий, В.М. Основы численных методов / В.М. Вержбицкий. – М.:
Высшая школа, 2002. – 840 с.
2. Бахвалов Н.С. Численные методы / Н.С. Бахвалов, Н.П. Жидков, Г.М.
Кобельков. – М.: Бином, 2004. – 634 с.
12