Построение интерполяционных сплайнов

3.2.1. Общая постановка задачи. Решение методом неопределённых коэффициентов

 

Рассмотрим решение одномерной задачи. Обозначим величины параметра t в узловых точках `P<sub>0</sub>, `P₁,...,`Р_n через t₀, t₁,...,t_n. При этом t₀ <t₁<...<t_n. На каждом i -ом сегменте интерполирования [t_i-1, t_i ] (i=1,...,n) искомый кусочно-полиномиальный сплайн задается отдельным многочленом S_i(t) таким образом, чтобы:

1) в узловых точках t_i-1, t_i многочлен S_i(t) совпадал со зна-чениями y_i-1, y_i;

2) в тех точках t_i-1, t_i, которые являются внутренними, у многочлена S _i (t) совпадали производные от 1 -го до m -го

порядков с аналогичными производными соседних много-членов.

В математической форме выражение искомой кусоч-но-полиномиальной функции ` P(t) можно представить в виде:

`P(t)=`S_i(t) при t_i-1 £ t £ t_i; i=1,2,…,n – условие кусочности сплайна` P(t);

`S<sub>i</sub>(t<sub>i-1</sub>)=`P_i-1;`S<sub>i</sub>(t<sub>i</sub>)=`P_i, i=1,2,…,n – условия прохождения через заданные точки;

`S<sub>i</sub><sup>-j</sup>(t<sub>i</sub>)=`S_i+1^-j(t_i), j=1,…,m, i=1,…,n-1 - условия на глад-кость во внутренних точках. (3.6)

Поскольку методика интерполирования одинакова по всем координатам (x,y,z), то рассмотрим интерполирование только по одной координате х. Соответствующие сплайны обозначим нижним индексом х. Постановка задачи имеет вид:

x(t)=S_{х i}(t); t_i-1 £ t £ t_i, i=1,…,n;

2. S_{х i} (t_i-1)=x_i-1; S _{х i} (t_i)=x_i , i=1,…,n;

3. S_{х i} ^j(t_i)=S_{х (i+1)} ^j(t_i), j=1,…,m, i=1,…,n-1. (3.7)

Число геометрических условий связи, налагаемых на искомую кусочно-полиномиальную функцию x(t) системой (3.7), складывается из 2n (во второй строке) и m(n-1) (в третьей). Суммарное их число равно:

S_s = 2n + m(n - 1) = n(m + 2) - m.

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

S _K = n(k+1).

Для того чтобы решение существовало и было един-ственным, необходимо выполнение условия S _S = S _K. Для этого:

1) принимают степень многочленов сплайна k = m+1;

2) дополнительно вводят m условий связи, которые называют краевыми условиями, поскольку обычно их задают в конечных точках ` P<sub>0,</sub> ` P_n, исходя из специфики решаемой задачи.

Наиболее распространенным случаем является ин-терполяция при m = 2 (гладкость 2 -ой степени). В этом случае необходимо использовать кубические сплайны (k = 3) и добавить 2 дополнительных краевых условия на производные. Чаще всего налагают следующие виды крае-вых условий.

1) S₁ ¹(t₀)=f ¹(t₀), S_n ¹(t_n)=f ¹(t_n) - условия на первые производные. Применяют в том случае, когда сплайн в концевых точках необходимо гладко соединить с некоторой кривой f(t).

2) S₁ ¹¹(t₀)=f ¹¹(t₀), S_n ¹¹(t_n)=f ¹¹(t_n) - условия на вторые производные. Обычно они используются для уменьшения осцилляции сплайна.

3) S₁ ¹(t₀) = S_n ¹(t_n); S₁ ¹¹(t₀) = S_n ¹¹(t_n) - условия для периодических функций, заключающиеся в том, что в конечных точках интервала равны не только значения, но и две первых производных сплайна.

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

Допустим, каждый i –тый участок сплайна имеет вид

S_xi(t)=C_i0+C_i1t+...+C_ikt^k .

Каждое из условий второй строки системы (3.7) даёт по одному уравнению:

C_i0 + C_i1 t_i-1 +... + C_ik t_i-1 ^k = x_i-1 ;

C_i0 + C_i1 t_i +...+ C_ik t_i ^k = x_i .

Рассмотрим равенство первых производных во внутренних точках S_хi ¹(t_i)=S_х(i+1)¹(t_i) в третьей строке системы (3.7). После взятия первых производных по t в S_xi и S_x(i+1) и подстановки t = t_i получим уравнение

C_i1+2C_i2t_i+…+kC_ikt_i^k-1 = C_(i+1)1+2C_(i+1)2t_i+…+kC_(i+1)kt_i^k-1.

Перенося все слагаемые, содержащие неизвестные коэффициенты, в левую часть, получим:

C_i1+2C_i2t_i+…+kC_ikt_i^k-1-C_(i+1)1-2C_(i+1)2t_i-…-kC_(i+1)kt_i^k-1=0.

При j = 2 из условия S_xi¹¹(t_i) = S_x(i+1)¹¹(t_i) после

а) определения вторых производных,

б) подстановки в них узловых точек и

в) переноса всех слагаемых, содержащих неизвестные (коэффициенты полдиномов сплайна),получим:

2C_i2+…+k (k-1) C_ikt_i^k-2- 2C_(i+1)2-…-k(k-1)C_(i+1)k t_i^k-2 = 0.

Для более высоких степеней производных действия аналогичны.

После объединения всех полученных уравнений полу-чаем линейную систему вида

A`C=`Х,

в которой:

- неизвестный вектор `C задаёт полный набор коэффи-циентов, входящих в многочлены сплайна,

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

- свободный вектор Х содержит заданные узловые значения многочленов сплайна, а также значения производных из краевых условий.

В общем случае решение системы (вектор `C коэф-фициентов, входящих в многочлены сплайна) может быть найдено любыми общими методами, например, по методу Гаусса. Однако с учётом диагонального вида матрицы А системы на практике используют более оптимальные специальные методы решения.

Пример. Построить сплайн P(t), проходящий через три точки (х(t₀)=x₀ , x(t₁)=x₁, x(t₂)=x₂ ), который во внутренней точке x(t₁)=x₁ имел бы гладкость 2 -ой степени (Рис.3.1).

x

S₁(t) S₂(t)

x₁

x₂

x₀

0

t₀ t₁ t₂ t

Рис.3.1

Решение. Интервала два, поэтому искомые многочлены на них обозначим S₁(t), S₂(t). Зададим их следующим образом:

S₁(t)=C₀₁+C₁₁t+C₂₁t²+C₃₁t³; S₂(t)=C₀₂+C₁₂t+C₂₂t²+C₃₂t³.

Условия I рода имеют вид:

S₁(t₀)=x₀; S₁(t₁)=x₁; S₂(t₁)=x₁; S₂(t₂)=x₂.

Условия II рода при m=2 и n=3 будут следующими:

S¢₁ (t₁)= S¢₂ (t₁); S¢¢₁ (t₁)= S¢¢₂ (t₁).

Поскольку n(k+1) - n(m+2) = m = 2, то для одно-значности решения задачи необходимо задать еще два крае-вых условия. Зададим оба условия в начальной точке:

S¢₁(t₀)=0; S¢¢₁ (t₀)=0.

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

C₀₁+C₁₁t₀+C₂₁t²₀+C₃₁t³₀=x₀;

C₀₁+C₁₁t₁+C₂₁t²₁+C₃₁t³₁=x₁ ;

C₀₂+C₁₂t₁+C₂₂t²₁+C₃₂t³₁=x₁ ;

C₀₂+C₁₂t₂+C₂₂t²₂+C₃₂t³₂=x₂ .

Дифференцируя многочлены по параметру и подстав-ляя в первые и вторые производные t = t₁, из условий II рода получаем:

C₁₁ + 2C₂₁t₁ + 3C₃₁t₁² - C₁₂ - 2C₂₂t₁ - 3C₃₂t₁² = 0;

2C₂₁ + 6C₃₁t₁ - 2C₂₂ - 6C₃₂t₁ = 0.

Из краевых условий получаем:

C₁₁ + 2C₂₁t₀ + 3C₃₁t₀ ² = 0;

2C₂₁ + 6C₃₁t₀ = 0.

Объединим все полученные уравнения в систему ви-да A`C=`Х, где

Коэффициенты обоих многочленов сплайна находят-ся одновременно при решении данной системы:

`C = A <sup>–1</sup>` Y.

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

 

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

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

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

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

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

Зафиксируем внутренний узел 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.а) полиномы Эрмита j_1i (t) - j_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 и введя коэф-фициент m_i = h _i/h_i+1, представим уравнение связи величин у¢_i-1 , у¢_i, у¢_i+1 (i=1,..., n-1) в виде:

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

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

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

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

m _n-1 у¢_{n -2} + 2(1+m _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,5m₁)у¢₁ + у¢₂ = b₁ +1,5m₁ (у₁ - у₀ ) / h₁ ,

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

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

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

m _n-1 у¢_{n -2} + С _n-1 у¢_n-1 = b¢_n-1 ,

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

С₁ =2(1+m₁), b¢₁ = b₁ - m₁ f ¹(х₀),

С _n-1 =2(1+m_{n -1}), b¢_n-1 = b_n-1 - f ¹(х_n),

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

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

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

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

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

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

m _n-1 у¢_{n -2} + С _n-1 у¢_n-1 = b¢_n-1 ,

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

А`Y¢ =`B, (3.12 б)

где

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

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

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

у¢_i = a_i у¢_i+1 + b_i. (3.13)

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

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

у¢_n-1 =(b¢_n-1 - m_n-1 b_n-1) / (c_n-1 + m_n-1 a_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 при фиксированных зна-чениях по х.

⇐ Предыдущая123Следующая ⇒

Поделиться:

Воспользуйтесь поиском по сайту: