Интерполяция сплайнами

Полиномиальная интерполяция не всегда дает удовлетворитель­ные результаты при аппроксимации зависимостей. Так, например, при представлении полиномами резонансных кривых колебательных си­стем большая погрешность возникает на концах ("крыльях") этих кри­вых. Несмотря на выполнение условий Лагранжа в узлах интерполя­ционная функция может иметь значительное отклонение от аппрокси­мируемой кривой между узлами. При этом повышение степени интер­поляционного полинома приводит не к уменьшению, а к увеличению погрешности. Возникает так называемое явление волнистости [4].

Для проведения гладких кривых через узловые значения функции чертежники используют упругую металлическую линейку, совмещая ее с узловыми точками. Математическая теория подобной аппроксима­ции развита за последние тридцать лет [5, 6] и называется теорией сплайн-функций. Разработано и обширное программное обеспечение для применения сплайнов в различных областях науки и техники [1,7].

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

(2.18)

Функцию и будем использовать для аппроксимации зависимости f(x), заданной в узлах ,...,х_пзначениями .

Если в качестве функции выбрать полином, то в соответствии с уравнением (2.18) степень полинома должна быть не выше третьей. Этот полином называют кубическим сплайном, который на интервале записывают в виде

, (2.19)

где и — коэффициенты сплайна, определяемые из дополни­тельных условий; i= 1,2,..., п— номер сплайна.

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

Коэффициенты сплайнов определяются из условий сшивания со­седних сплайнов в узловых точках:

1) Равенство значений сплайнов и аппроксимируемой функ­ции f(x) в узлах — условия Лагранжа

. (2.20)

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

(2.21)

 

(2.22)

Кроме перечисленных условий необходимо задать условия на кон­цах, т.е. в точкахx₀ иx_n. В общем случае эти условия зависят от кон­кретной задачи. Часто используют условия свободных концов сплай­нов. Если линейка не закреплена в точках вне интервала[x₀, x_n], то там она описывается уравнением прямой, т.е. полиномом первой степени. Следовательно, исходя из условий (2.22), т.е. непрерывности вторых производных сплайнов на концах интервала, запишем соотношения:

(2.23)

(2.24)

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

Получим алгоритм определения коэффициентов кубических сплай­нов из условий (2.20) – (2.24). Условия (2.20) в узлах x_i-1и x_i после подстановки i-гo сплайна принимают вид

a_i = ƒ_{i-1 ,} (2.25)

a_i + b_ih_i +c_ih_i² + d_ih_i³ = ƒ_i, (2.26)

гдеh_i= x_i— x_i-1, 1 ≤ i ≤ n.

Продифференцируем дважды сплайн (2.19) по переменной x:

φ′_i(x) = b_i+ 2c_i (x— x_i-1) + 3d_i (x— x_i-1)², (2.27)

φ″_i(x) = 2c_i+ 6d_i (x— x_i-1). (2.28)

Из условий непрерывности производных (2.21) и (2.22) при пере­ходе в точке x_iот i-го к i₁сплайну с учетом выражений (2.27) и (2.28) получим следующие соотношения:

b_i + 2c_i h_i +3d_ih_i²= b_i+1 , (2.29)

c_i + 3d_ih_i= c_i + 1. (2.30)

И, наконец, из граничных условий (2.23) и (2.24) на основании вы­ражения для второй производной (2.28) получим, что

c₁ = 0, (2.31)

c_n+ 3d_nh_n= 0. (2.32)

Соотношения (2.25), (2.26), (2.29) — (2.32) представляют собой полную систему линейных алгебраических уравнений относительно ко­эффициентов сплайнов a_i, b_i,c_i иd_i. Но прежде чем решать эту систе­му, выгодно преобразовать ее так, чтобы неизвестными была только одна группа коэффициентов c_i.

Из уравнений (2.30) коэффициенты d_iвыразим через коэффициен­ты c_i:

(2.33)

 

Объединяя уравнения (2.25) и (2.26) с уравнением (2.33), предста­вим коэффициенты b_iтакже через коэффициенты c_i:

f _i — f _i-1 (c _i+1 + 2c _i)h

b_i = ———— - ————— . (2.34)

h_i 3

После подстановки выражений (2.33) и (2.34) в соотношение (2.29) получим уравнение, в которое входят только неизвестные коэффици­енты c_i. Для симметричности записи в полученном уравнении умень­шим значение индекса i на единицу:

(2.35)

где 2≤ i ≤n.

При i = n, учитывая условие свободного конца сплайна, в уравне­нии (2.35) следует положить

c_n+1 = 0. (2.36)

Таким образом, уравнение вида (2.35) вместе с условиями (2.31) и (2.36) образует систему линейных алгебраических уравнений для определения коэффициентов c_i. Коэффициенты d_i и b_iвычисляют­ся после нахождения c_i по формулам (2.33) и (2.34), коэффициенты a_iравны значениям аппроксимируемой функции в узлах в соответствии с формулой (2.25).

В каждое из уравнений типа (2.35) входит только три неизвестных с последовательными значениями индексов c_i-1,c_i,c_i+1. Следователь­но, матрица системы линейных алгебраических уравнений относитель­но c_i является трехдиагональной, т.е. имеет отличные от нуля элементы только на главной и двух примыкающих к ней диагоналях. Для реше­ния систем с трехдиагональной матрицей наиболее эффективно при­менять так называемый метод прогонки, являющийся частным случаем метода Гаусса.

Рассмотрим алгоритм метода прогонки применительно к решаемой задаче. Для сокращения записи уравнение (2.35) представим в виде

h_i-1c_i-1+s_ic_i + h_ic_i+1 = r_i, (2.37)

гдеs_i = 2 (h_i-1 - h_i).

r_i= 3 . (2.38)

Так же, как и метод Гаусса, метод прогонки разделяется на два эта­па — прямой и обратный ходы. В процессе прямого хода метода про­гонки вычисляют значения коэффициентов линейной связи каждого предыдущего неизвестного c_iс последующим c_i+1.

Из уравнения (2.37) при i = 2 с учетом граничного условия (2.31) установим связь коэффициента с₂ с коэффициентом с₃:

 

c₂ = k₂- l₂с₃ , (2.39)

где k₂ и l₂ — прогоночные коэффициенты,

 

Затем, подставляя выражение (2.39) в уравнение (2.37) при i = 3, получим линейную связь коэффициента с₃ с коэффициентом с₄:

c₃ = k₃- l₃с₄.

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

c_i = k_i- l_iс_i+1. (2.40)

В процессе выполнения прямого хода метода прогонки необходимо вы­числить значения всех прогоночных коэффициентов k_i и l_i, для кото­рых получим рекуррентные соотношения. Подставим формулу связи (i -1)-го и i-го коэффициентов

c_i-1 = k_i-1- l_i-1с_i

в уравнение (2.37), в результате получим

.

Сравнение последнего соотношения с формулой (2.40) позволяет за­писать рекуррентные соотношения для определения прогоночных ко­эффициентов k_iи l_i:

(2.41)

Учитывая граничное условие (2.31) и соотношение

c₁ = k₁- l₁с₂,

а также полагая c₂≠0, задаем начальные коэффициенты k₁ = 0 и l₁ = 0. Затем по формуле (2.41) вычислим все ппар прогоночных ко­эффициентов k_i и l_i..На основании соотношения c_n = k_n- l_nc_n+1 и граничного условия (2.36) получим, что

c_n = k_n. (2.42)

Далее последовательно применим формулу (2.40) при i = n — 1,n— 2,.. и вычислим значения искомых величин c_n-1, c_n-2,…, c_2.

Эта процедура называется обратным ходом метода прогонки.

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

После определения всех коэффициентов с_iдругие коэффициенты сплайна вычисляются по формулам (2.25), (2.33) и (2.34). Аппрокси­мирующую функцию φ(x)можно рассчитать с помощью соотношения (2.19) в любой точке xна интервале [х₀, х_n].