Алгоритмы и программы для интерполяции табличных функций

Интерполяционный полином Лагранжа L(x) представляет собой тот же самый полином (1), но записанный в другой форме, исключающей необходимость решения СЛАУ. Он представляется в виде линейной комбинации элементарных полиномов n-ной степени li(x), в которой коэффициентами служат заданные значения yi:

L(x) = y0l0(x) + y1l1(x) + y2l2(x) + ... + ynln(x)   (2)

При таком определении интерполирующего полинома L(x) каждый элементарный полином li(x) должен удовлетворять следующим условиям в узловых точках:

li(xi) = 1,               li(xj) = 0 (j-i)

Нетрудно убедиться, что этим условиям удовлетворяет полином

(x-x0)(x-x1)(x-x2)...(x-xj)..(x-xn)

li(x) = ------------------------------------------(xi-x0)(xi-x1)(xi-x2)...(xi-xj)...(xi-xn)

Окончательно формула полинома Лагранжа может быть представлена в виде:

n

--n           | |(x-xj)

---           j=0

L(x) =  >   yi ------------;                                 (3)

---       n

i=0          --| |(xi-xj)

j=0

j-i

5. Сплайн-интерполяция

Кубический сплайн представляет собой полином третий степени вида:

Si(x) = ai + bi(x-xi-1) + ci(x-xi-1)2 + di(x-xi-1)3,   (4)

коэффициенты которого на каждом элементарном интервале  (xi-1,xi)

имеют разные значения. В принятой системе индексации номер сплайна совпадает с номером правой границы элементарного интервала, на котором он работает.

Далее приводится словесно-формульное описание алгоритма определения коэффициентов кубических сплайнов для функции y(x), заданной набором n+1 точек xi, yi (i=0,1,2,3,...,n). При этом получается n отрезков интерполяции и, соответственно, n сплайнов Si (i=1,2,3,...,n), которые имеют 4n коэффициентов. Для определения последних необходимо 4n уравнений.

Обычно для сплайн-интерполяции используют равноотстоящие по x точки, т.е. все элементарные интервалы равны: xi+1 - xi = h. Коэффициенты ai, bi, ci и di выбираются таким образом, чтобы интерполирующие сплайны имели одинаковые производные первого и второго порядка в точках сшивки. Это позволяет получить очень качественный результат интерполяции, который может быть использован для различных целей.

Эти уравнения можно получить из условия равенства значений интерполируемой функции и сплайна в узловых точках:

Si(xi-1) = yi-1;  Si(xi) = yi;              (i=1,2,3,...,n)             (5)

а также сшивки на границах элементарных интервалов первых производных сплайнов:

( dSi-1 )   ( dSi )

2-------2 = 2-----2 ;           (i=2,3,...,n)       (6)

9 dx   0xi 9 dx  0xi

и вторых производных:

( d2Si-1 )   ( d2Si )

2--------2 = 2------2;          (i=2,3,...,n)       (7)

9 dx2 0xi 9 dx2 0xi

Недостающие два уравнения получаются из условия "свободного закрепления концов":

( d2S1 )   ( d2Sn  )

2------2 = 0;    2-------2 = 0                       (8)

9 dx2 0x0 9 dx2 0xn

Подстановка (4) в уравнения (5) - (8) и некоторые преобразования приводят к следующей системе уравнений:

(

2 ai = yi-1;                                        (i=1,2,...,n)

2 ai + bih + cih2 + dih3 = yi;   (i=1,2,...,n)

* bi + 2cih + 3dih = bi+1;                 (i=1,2,...,n-1)                 (9)

2 ci + 3dih = Ci+1;                           (i=1,2,...,n-1)

2 c1 = 0;

2 cn + 3dnh = 0,

9

где h = xi+1 - xi = xi - xi-1 - расстояние между соседними точками (шаг по x).

Коэффициенты ai определяются просто:

ai = yi-1;   (i=1,2,...,n)       (10)

Затем можно определить коэффициенты ci. Для этого из уравнений системы (6) необходимо исключить ai, bi и di, что приводить к выражению:

yi-1 - 2yi + yi+1

ci + 4ci+1 + ci+2 = 37-------------------; (i=1,2,...,n-1) (11)

h2

Эта запись при развертывании индекса i в указанных пределах превращается в СЛАУ относительно ci. Данная система имеет трехдиагональную матрицу коэффициентов, в главной диагонали ее содержится 4, а в прилегающих к ней - 1. Для решения таких систем обычно используется специальный метод - метод прогонки.

После определения коэффициентов ci можно определить bi:

yi - yi-1           h

bi = ----------- - ---7(ci+1 + 2ci)                                       (12)

h               3

и di:

di+1 - di

di = -----------            (i=1,2,3,...,n)                          (13)

3h

6. Используемые подпрограммы

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

Nint - число интервалов интерполяции;

XA, YA - одномерные массивы - таблицы точек X и Y, по которым проводится интерполяция. Тип данных этих массивов - Vec100, определенный в модуле def_type; нумерация элементов в них начинается с нулевого. Число точек на единицу больше числа интервалов.

Выходные параметры SplineCoeff - одномерные массивы коэффициентов сплайнов A, B, C и D (тип Vec100). Переменная Err представляет собой код ошибки: Err=0 при благополучном завершении интерполяции. Остальные значения Err объяснены во входном комментарии подпрограммы.

Подпрограмма SplineCoeff находится в файле spline.pas (spli- ne.for, spline.cpp), который в паскалевском варианте представляет собой модуль