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

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

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

Формулировка основной системы линейных уравнений

При описании задачи интерполяции данных квадратичным сплайном S ₂( t ) дефекта один можно воспользоваться следующей схемой формулировки основных уравнений. Пусть носителем такого сплайна является интервал [0, n ]. Так как на каждом из интервалов разбиения множества [0, n ] на интервалы единичной длины сплайн, являющийся многочленом второго порядка, определяется тремя коэффициентами, то его будут описывать на всем интервале 3 n параметров, которые и подлежат определению.

Если учесть, что должны быть выполнены условия непрерывности сплайна S ₂( t ) и его первой производной в промежуточных целочисленных точках основного интервала [0, n ], то это накладывает 2( n –1) условий на параметры сплайна. В результате остаются неопределенными n +2 параметра. Далее, если для каждой средней точки интервалов единичной длины заданы значения анализируемой функции, то неопределенными останутся два параметра, которые обычно находятся из так называемых краевых условий (см. [1]). Подчеркнем, что в рассматриваемой схеме значения интерполируемой функции считаются заданными в полуцелых точках.

Рассмотрим способ описания сплайна S ₂( t ), аналогичный тому, который приведен в монографии [1] для кубических сплайнов. Для этого положим, что первая производная сплайна S 2 ( t ) на ин терва ле [ j , j +1] может быть представлена в виде ( j = 0, n - 1)

5 2° ( t ) = m j ( j + 1 - t ) + m j _{+ l} ( t - j ),                 (1)

где величины m _j определяют значения производных в целочисленных точках интервала [0, n ]. После интегрирования формулы (1) получаем, что

5 2 ( t ) = - m j

( j + 1 - t )²         ( t - j )²

-------------+ m....--------- + c .

2           j 2 j

Воспользовавшись условием 5 ₂ ( j + 0.5) = f j , из

(2) определяем постоянную интегрирования c j , что позволяет представить сплайн S 2 ( t ) в виде

5 2 ( t ) = - m j

⁺ m j + 1

( j + 1 - 1 ) 2

-

+

( t - j ) 2 1

+ f j .

Заметим, что представление сплайна S 2 ( t ) в виде (3) обеспечивает непрерывность как самого сплайна, так и его первых производных в узлах рассматриваемой сетки.

Если использовать условия непрерывности сплайна S ₂( t ) в целых точках, то получается следующая система уравнений для наклонов сплайна m j :

m j - 1 ^{+ 6} m j ⁺ m j + 1 = ⁸⁽ f j ^- f j - 1 ) = d j ,             ⁽⁴⁾

где j = 1, n - 1, поскольку такие уравнения имеют место только для внутренних точек рассматриваемого целочисленного интервала [0, n ]. Величины d j с точностью до множителя дают наклон в данных, относящихся к полуцелым точкам. Систему уравнений (4) следует дополнить краевыми условиями, которые позволят определить два свободных параметра m 0 и m_n .

Граничные условия будут выбраны с использованием дополнительных уравнений, записываемых в следующей форме:

6 m ₀ + m ₁ = d ₀ , m _{n -1} + 6 m _n = d n . (5)

Здесь подлежат определению величины d ₀, d_n . Эти величины будут выбраны после построения решения полученной системы уравнений (4) и (5). Формально же считаем, что величины d 0 , d n нам заданы. Подчеркнем, что выбор дополнительных уравнений в форме (5) удобен тем, что получающаяся при этом матрица системы линейных уравнений будет симметричной теплицевой матрицей, т. е. такой, на каждой диагонали которой стоят одинаковые значения.

Метод решения теплицевой системы линейных уравнений

Решение системы уравнений (4), (5) будет строиться следующим образом. Дополним произвольно набор искомых наклонов m _j и значений правых

to                     А то               to                        p

E mpzp=^ E zp Г E Xj—pdj + E Xp"d p =—TO            A p = —TO     [ j = p + 1             j=—TO

частей уравнений d j до множеств из бесконечного

Заметим, что в результате приравнивания коэффициентов при степенях z получаем решение в виде

числа

элементов { ш , } , { d } . Относительно

вновь вводимых величин будем предполагать, что вновь вводимые наклоны m j имеют нулевые значения. Дополнительные величины d j выбираются рав-

Г m =-XГ V Xj—pd, + У Xp—jd, p   jj

^A [ j = p + 1              j =—TO

Отметим также, что формула (13) фактически может быть получена интегрированием по углу ф ком-

ными левым частям дополнительно вводимых уравнений типа уравнений (4). С помощью полученных

наборов { m j }    и { d j }    формально определим

плексных тригонометрических рядов (12) при условии z = exp( i ф ). При сделанных же ранее предположениях

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

относительно свойств бесконечной системы уравнений справедливость приведенных действий очевидна.

Полагая вводимые дополнительные параметры

{ d j }

— 2

—то

и { d j } 2 равными нулю, рассмотрим фор-

мулу (13) с целью выяснения возникающих равенств для определения величин m 0 , m n . С учетом нулевых

m - 1 + 6 m 0 + m 1 = d 0 , m n -1 ^{+ 6} m n ⁺ m n + 1 = d n ,

значений для наборов

{ d j }

и d

1 ^jln + 2

для индек-

в которой использовано равенство нулю дополнительно вводимых величин наклонов. В итоге получаем формально бесконечную систему уравнений, которая имеет вид mj–1+6mj+mj+1=dj, (7)

где j = —to , to .

Для решения полученной системы (7) применим следующий прием. После умножения каждого из уравнений (7) на переменные z^j (| x |=1) соответственно и суммирования полученных равенств получаем функциональное уравнение

са p = –1 потребуем выполнения условия

⁽ d — 1 = m 0 , d n + 1 = m n ^):

—X I                               I m—1 = —

A [j=0

которое дает одно из уравнений для определения величин m0, mn. Для остальных отрицательных индексов p<–1, соответственно, имеем

1 I n mp = T lE Xj—pdj+X—1-pm 0 +Xn+1—pmn

^A[ j =0

co                     co                                                                    co

E djZ1 = E zj (mi—i+6 mi+mj+i)=Q2( z) E mpzp, j=—TO           j=—TO                                             p=—TO в котором использовано обозначение

Q2 ( z) = z "* + 6 + z = — 1X (1 — X z )(1 — Xjz),      (9)

x—p—1 Гn|

=-----

A

Заметим, что равенства (15) являются следствием формулы (14).

Для положительных значений индексов формулы получаются аналогично предыдущему случаю. Так, при p = n+1 из формулы (13) имеем условие, которое должно быть выполнено:

где корни X± = —3±2V2 многочлена (z2+6z +1) об

ладают свойствами Z++Z_-=-6, X₊X_-=1, 1-X²=A,

|X₊|=|X|<1. С учетом этих свойств имеем

I                                                                                             .I

X } X "1 n+1—j J 1       1 "1 n+2 I A mn+1 = TiE X   dj+mn+X m0 r = 0.

^AIj=0                             J

Q2(z)

—X Г 1      Xz

-------J--1--

A [1 — X / z 1 — Xz

Для индексов p >n +1, соответственно, имеем

Г                                                    1

mp =-X i E X p^—d +X^p—n^—1 mn +X^p+1m, Г =

A

¹-¹ [ j =—TO                                                   J

Полученные соотношения позволяют вместо формулы (8) написать решение системы (7) в следующем виде:

p—n Г p

—— Г У Xn+1—jd, + m„ +Xn+2mn Г = 0.(17)

jn0

A I j = — TO

Q2(z)

co                     co

E dX = E mpzp = j=—TO          p=—TO

—x Г у (x ч a [^E1 z J

co                    co

+E (xz) k\E zdr

k =0         j=—~

Далее, преобразовав уравнение (11), получаем для нахождения решения системы формулу

Тем самым для обеспечения рассматриваемых условий тождественности обращения в ноль вспомогательных переменных mj должны быть выполнены условия:

n

У Xj +1 d + m,. + Xn+2 m = 0, j0              n j=0

n

У X^{n +1—j}d + m + Xⁿ⁺²m„ = 0, jn    0

j=0

И.И. Орлов которые позволяют определить величины m0, mn в зависимости от величин dj.

Рассмотрим теперь формулы (13) при p = 0:

___ I n                                                                     I m0 = — < V Xj+1 d j + Xn+2mn + X2m0 к

A Ij=0                          J

С учетом первого уравнения (18) и свойств собственных значений формула (19) преобразуется к тождеству

m mn = — (1 _ X2 ) = mn.                       (20)

A

Аналогично предыдущим соотношениям, для p = n получаем формулу mn =— lX2mn + VXn j+1 dj + Xn+2m0 [,        (21)

A I      i0

которая, с использованием второго уравнения (18), приводится к тождеству mn = m~ (1 _X2) = mn .

A

Теперь из формулы (13) при p = 1, n - 1 получаем mp

-X I1

— J У Xj^_pd. +УХp^_jd. +X^n+1_pm +Xp⁺¹m J.

j                   j                 n°

Ij=p+1              j=0

Напомним, что величины m0, mn определяются из системы (18). Эти решения задаются формулами:

II m0 = J _У Xj+1 d, + Xn+2 У Xn+1_j d, J, 0 1 a 2 n+4                   jj

• 1    _X I j=0                j=0

• 1    I nn

m„ = J^Xn+1_idi +Xn+2У Xj+1 d, J. n           2n+4                    jj

• 1    X ^ j=0                     j=0

  
  

Заметим, что для больших n формулы (24) можно упростить с учетом малости вторых слагаемых в них вследствие неравенства |λ|<1. Эти формулы будут иметь вид

n m0 ^ _VXj+1 dj, j=0

n mn = _E X"+^'dj .

j=0

Для значений индекса p, принадлежащих интервалу [0, n] и удаленных от его концов, величины mp можно приближенно представить в виде разложения по степеням λ. В этом случае получим приближенное разложение mp -{_Xdp +X2dp+1 + dp_ 1 +X3dp+2 +Xdp_2 + ...}.(26)

Отсюда следует, что имеет место свойство квазилокальности, заключающееся в том, что основной вклад в значение величины mp вносят те dj, индексы которых мало отличаются от индекса p. Это важное общее свойство и оправдывает в ряде случаев использование вместо интерполяционных аппроксимирующих сплайнов.

Набор величин dj и формулы (24) вместе с (23) полностью определяют решение системы (4), и оста- лось задать d0, dn, которые ранее не были определены. Это может быть сделано различными способами, которые обусловлены требованиями к выбранному варианту интерполяционного сплайна.

Возможен такой способ задания дополнительных условий. Продолжив формально сплайн на один интервал влево и вправо, считаем, что величины m–1, mn+1 имеют нулевые значения. Значения же интерполируемой функции в точках t = –0.5 и t = n+0.5 определим с использованием формул линейной экстраполяции, положив, что f–1 = 2f0 – f1, f_n+1 = 2f_n–1 – f_n–2. В таком случае добавляемые уравнения (5) будут иметь в правых частях следующие значения: d0 = 8(f1–f0), dn = 8(fn–1 – fn–2), которые получаются из значений интерполируемой функции по тем же формулам, что и остальные величины правой части системы линейных уравнений (4).

Другой способ задания величин d0, dn может быть осуществлен следующим образом. После того как определены величины {mj}_q через параметры {dj}0, мы можем считать, что значения сплайна в крайних точках интервала [0, n] получаются линейной интерполяцией по паре ближайших к этим точкам значений. Эти условия дают S2(0)=1.5f0–0.5f1, а S2(n)=1.5fn–0.5fn–1. Используя эти условия, из формул типа (3) получаем, что должны быть выполнены условия

S₂ (0) = —m₀ —m_v + f₀ = 1.5f₀ _Q.5f, 2          ₈0   ₈1     0         0         1

„ , , 3         1            .

s2(n) = - mn + - mn_ 1 + fn_ 1 = 1.5 fn _0.5fn_ 1.

Так как все mp являются линейными комбинациями величин {dj }0, то из полученных условий (27) может быть получена система двух линейных уравнений на величины d₀, d_n. Очевидно, что такой способ менее удобен, чем приведенный выше метод задания величин d0, dn непосредственно с использованием методики продолжения сплайна тем или иным способом вне основного рассматриваемого интервала. Возможны также и иные принципы определения граничных условий. Примеры такого рода дополнительных условий можно найти, например, в монографии [1].

Заключение

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

Другой общий метод обращения конечных теплицевых матриц, отличный от рассмотренного в этой работе, содержится в монографии И.Ц. Гохбер-га и И.А. Фельдмана [2], посвященной методам анализа уравнений в свертках. При получении результатов по обращению конечных теплицевых матриц существенно то свойство, что многочлены типа (10)

не имеют корней на единичной окружности, что в рассмотренном нами случае легко проверяется.

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

Использованный в рамках рассмотренной выше схемы решения трехдиагональной теплицевой системы уравнений метод основан на свойствах корней многочленов Q2(z), которые связаны с многочленами Эйлера–Фробениуса. Эти многочлены изучались в работе [3] в связи с задачами функциональной интерполяции. Эти же многочлены, в рамках подхода, отличного от использованного в работе [3], рассматривались в задаче построения периодических сплайнов на равномерной сетке [4].

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

В связи с этой аналогией заметим, что изложенный выше метод решения систем линейных уравнений дает дискретный аналог не только для решений обыкновенного дифференциального уравнения, но и методику нахождения дискретного аналога функции Грина. Действительно, если задавать правые части рассматриваемой системы в виде dj = δj, k (здесь δj, k – символ Кронекера), то из формул (23) с соответствующими значениями для m0, mn можно получить явные выражения для дискретного аналога функции Грина.