3.2.2. Кубические интерполяционные сплайны

Как уже отмечалось, на практике наибольшее приме-нение нашли сплайны нечётной степени. При k=1 интерпо-ляционные сплайны совпадают с локальными, поскольку отсутствуют условия на производные во внутренних точках. В случае k=3 локальные сплайны более точно приближают исходную функцию (с точностью до первых производных), однако имеют гладкость в узлах только первой степени. Интерполяционные кубические сплайны менее точно при-ближают функцию (только её значения), но имеют глад-кость в узлах степени 2. Также они обладают следующим важным свойством. Если наложить на сплайн краевые усло-вия S (х₀)= S (х_n)=0, то он будет минимизировать функ-ционал

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

Рассмотрим построение интерполяционного сплайна. В предыдущем разделе дано общее решение по методу не-определённых коэффициентов, однако в данном случае задачу можно решить намного проще с использованием кубических сплайнов Эрмита (локальных). Данные сплайны имеют гладкость степени 1 и у них помимо значений у_i во всех узлах x_i заданы также величины первых производных у_i . Поскольку у интерполяционного сплайна на первые производные никаких исходных условий не накладывается,

80

то непрерывность второй производной во внутренних узлах

и выполнение краевых условий можно обеспечить за счёт варьирования именно этих величин.

Зафиксируем внутренний узел x_i (i=1, . . . , n-1). Усло-вие непрерывности второй производной в нём принимает вид: S_i (x_i) = S_i+1 (x_i). Переходя к безразмерным норми-рованным переменным, представим его в форме: h²_i S_i (1) = h²_i+1 S_i+1 (0), h_i = x_i - x_{i -1} ; h_i+1 = x_i+1 - x_i . Найдем выражения для производных. Подставляя в формулу (3.5.а) полиномы Эрмита _1i (t) - _4i (t), , получим:

S_(i+1)(t) = у_i (1-3t²+2t³)+ у_i+1 ( 3t² –2t³)+ h_i+1 у_i (t - 2t²+ t³)+ h_i+1 у_i+1 ( - t ² + t ³ ).

Отсюда следует:

S_(i+1)(t)= у_i (-6+12t) /h²_i+1 + у_i+1 ( 6 –12t) / h²_i+1+ у_i ( - 4+ 6t) / h_i+1+ у_i+1 ( - 2 + 6t ) / h_i+1,

S_(i+1)(0)= -6у_i / h²_i+1 + 6 у_i+1 / h²_i+1 - 4 у_i / h_i+1 - 2 у_i+1 / h_i+1.

Аналогично рассматриваем участок [x_i-1, x_i] :

S_i (t)= у_i-1 (1-3t²+2t³)+ у_i ( 3t² –2t³)+ h_i у_i-1 (t - 2t²+ t³)+

h_i у_i (- t ² + t ³ ),

S_i(t)= у_i-1 (-6+12t) / h²_i+ у_i ( 6 –12t) / h²_i + у_i-1 ( - 4+ 6t) / h_i + у_i ( - 2 + 6t ) / h_i,

S_i (1)= 6у_i-1 / h²_i - 6 у_i / h²_i + 2 у_i-1 / h_i+ 4 у_i / h_i.

Приравнивая вторые производные, группируя слага-емые и сократив на 2, получим:

у_i-1 /h_i + 2(1/h_i +1/h_i+1)у_i + у_i+1 /h_i+1 = 3[- у_i-1 /h²_i + (1/h²_i -1/h²_i+1 )у_i+ у_i+1 /h²_i+1].

Домножая обе части выражения на h_i+1 и введя коэф-фициент _i = h _i/h_i+1 , представим уравнение связи величин у_i-1 , у_i , у_i+1 (i=1, . . . , n-1) в виде:

_i у_{i -1} +2(1+_i )у_i + у_i+1 = b_i ; (i=1, . . . , n-1) (3.8)

где b_i =3[_i ( у_i - у_i-1 ) / h_i +( у_i+1 - у_i ) / h_i+1 ] - постоянный коэффициент.

81

Рассмотрим краевые условия 1) (п.3.2.1), при кото-рых заданы условия на первые производные в начальной и конечной точках сплайна: S₁ ¹(t₀)=f ¹(х₀), S_n ¹(t_n)=f ¹(х_n ). Их учёт приводит к тому, что в уравнениях для узлов x₁ и x_n вместо у₀ , у_n необходимо подставить, соответственно, заданные величины f ¹(х₀) и f ¹(х_n ). После переноса соответствующих слагаемых в правую часть эти уравнения примут вид:

2(1+₁ )у₁ + у₂ = b₁ - ₁ f ¹(х₀) ,

 _n-1 у_{n -2} + 2(1+ _n-1 )у_n-1 = b_n-1 - f ¹(х_n ). (3.9)

Краевые условия 2) (п.3.2.1) обычно для уменьшения осцилляции сплайна применяют в виде: S₁¹¹(х₀)=S_n ¹¹(х_n)=0. Раскрывая их, после сокращения на 2 получим следующие дополнительные уравнения:

-3у₀ / h²₁ + 3 у₁ / h²₁ -2 у₀ / h₁ - у₁ / h₁ = 0,

3у_{n -1} / h²_n - 3 у_n / h²_n + у_n-1 / h_n + 2 у_n / h_n = 0.

Для приведения системы уравнений к однотипному виду выразим из этих уравнений, соответственно, у₀ че-рез у₁, а у_n - через у_n-1 :

у₀ = -1,5 у₀ / h₁ + 1,5 у₁ / h₁ - 0,5 у₁ ,

у_n = -1,5 у_{n -1} / h_n + 1,5 у_n / h_n -0,5 у_n-1 . (3.10)

Подставляя найденные выражения, соответственно, в уравнения для узлов х₁ и х_n-1, получим их в следующей форме:

(2+1,5₁ )у₁ + у₂ = b₁ +1,5₁ (у₁ - у₀ ) / h₁ ,

 _n-1 у_{n -2} +(1,5+2 _n-1 )у_n-1 =b_n-1+1,5 (у_n-1 -у_n )/h_n . (3.11)

Очевидно, в обоих рассмотренных видах краевых условий уравнения (3.9),(3.11) для узлов х₁ и х_n-1 качест-венно одинаковы. Их можно представить как:

С₁ у₁ + у₂ = b₁ ,

 _n-1 у_{n -2} + С _n-1 у_n-1 = b_n-1 ,

где для условий 1-го вида:

С₁ =2(1+₁ ), b₁ = b₁ - ₁ f ¹(х₀) ,

С _n-1 =2(1+_{n -1} ), b_n-1 = b_n-1 - f ¹(х_n ) ,

82

для условий 2-го вида:

С₁ =2+1,5₁ , b₁ = b₁ + 1,5₁ (у₁ - у₀ ) / h₁ ,

С _n-1 =1,5+2_{n -1} , b_n-1 = b_n-1 +1,5 (у_n-1 - у_n ) / h_n .

В итоге система уравнений для определения значе-ний (у₁ , у₂ , . . . , у_n-1) примет вид:

С₁ у₁ + у₂ = b₁ ,

_i у_{i -1} + 2(1+_i )у_i + у_i+1 = b_i ; (i=2, . . . , n-2) (3.12 а)

 _n-1 у_{n -2} + С _n-1 у_n-1 = b_n-1 ,

В векторной форме:

АY =B , (3.12 б)

где

Матрица системы трёхдиагональна, поэтому для упро-щения её решения применяют специальную модификацию

83

метода Гаусса – метод прогонки. Он также имеет прямой и обратный ход.

При прямой прогонке для всех величин у_i (i=1, . . . , n-2) путём последовательного исключения находят их ли-нейные выражения через у_i+1 вида:

у_i = _i у_i+1 + _i . (3.13)

При обратной прогонке:

1) из совместного решения линейного выражения вида (3.13) для у_n-2 и последнего уравнения системы (3.12) находят величину у_n-1 :

у_n-1 =(b_n-1 - _n-1 _n-1 ) / (c_n-1 + _n-1 _n-1 ) ,

2) последовательно подставляя у_n-1 , у_n-2 , . . . в уравнения (3.13), определяют все искомые величины у_i .

Значения у₀, у_n на краях в случае условий 1-го рода являются нулевыми: у₀=у_n =0. При условиях 2-го рода их рассчитывают по формулам (3.10).

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

Рассмотренные выше кубические интерполяцион-ные сплайны, имеющие гладкость 2-й степени, называют сплайнами Фергюсона или сплайнами Шёнберга. Расчёт их значений производится так же, как и у сплайнов Эрмита по формулам (3.4 а,б).

Двухмерные сплайны Фергюсона, предназначенные для интерполирования поверхностей на прямоугольной сетке (x_i, у_j) (i = 0,1 ,..., n; j = 0,1 ,..., m), представляют собой такие же квадратичные формы, как и двухмерные сплайны Эрмита (3.5). Требуемые значения первых производных z^x_ij , z^у_ij в узлах сетки можно найти последовательно мето-дом одномерной прогонки – сначала z^x_ij при фиксирован-ных значениях по у , затем - z^у_ij при фиксированных зна-чениях по х.

84