5.1. Cплайн-апроксимація

Mатематичні сплайни беруть свій початок від тонких гнучких стрижнів, якими користувалися креслярі для проведення плавних кривих через задані точки. Стрижень закріплювався в точках (xi, yi) і набував форми кривої y(x) з мінімальною "енергією натягу", пропорційною .

Якщо перейти до математичного опису сплайну, то сплайн-функцією ступеня k з точками з'єднання x 0 1 буде функція y ( x ) , яка на відрізку [ x 0, xn ] має безперервні похідні до ( k -1) включно і на кожному з відрізків [xi-1, xi] дорівнює многочлену ступеня k. Далі нас цікавитимуть лише кубічні сплайни (k=3). Такий сплайн забезпечує збіг у вузлах з вихідною функцією та безперервність першої та другої похідних у точках з'єднання.

Насамперед визначаються перші похідні у всіх точках з'єднання. Вирішується тридіагональна система ( n -1) рівнянь з домінуючою головною діагоналлю. Вона може бути вирішена, якщо задати дві крайові умови. Програми, включені до Графору, дозволяють задавати чотири типи крайових умов:

а) відомі значення перших похідних на кінцях сітки;

б) відомі значення других похідних на кінцях сітки;

в) при побудові періодичного сплайн значення функції в першій і останній точках збігаються, а також і ;

г) якщо крайові умови не задані, вважається, що .

Після того, як побудована тридіагональна матриця (або доповнена тридіагональна матриця у разі в)) і вектор вільних членів, система рівнянь може бути вирішена. Для вирішення такої системи існує ефективний алгоритм, який у вітчизняній літературі називається методом прогону. В результаті розв'язання системи рівнянь виходить вектор першихпохідних у точках з'єднання.

Тепер значення y ( x ) для xi - 1 £ x £ xi визначається з многочлена

коефіцієнти якого легко знайти, оскільки відомі значення функції та перші похідні в точках з'єднання.

Таким чином, метод дозволяє виділити низку самостійних функцій:

формування тридіагональної матриці та вектора вільних членів відповідно до заданих крайових умов – підпрограма TDMP ,
рішення системи рівнянь, заданих матрицею та вектором вільних членів (метод прогонки) – підпрограма TRIDIG ,
обчислення коефіцієнтів кубічних многочленів всім відрізків – підпрограма SPLINE ,
обчислення значень сплайну на заданій рівномірній сітці – підпрограма SPLINT.

Зазначимо, що програма SPLINE при роботі використовує підпрограми TDMP і TRIDIG . Ці програми універсальні та їх застосування не обмежене сплайн аппроксимацією. Їхня універсальність має і глибший зміст. Зокрема, коефіцієнти для кубічних парабол можна обчислити, задаючи вектор похідних повністю. Це означає, що програміст на власний розсуд може модифікувати вектор похідних. При інтерактивній роботі таким чином можна підбирати відповідну форму кривої. Математичне обґрунтування сплайнів та їх аналіз викладено в [13].

Програма SPLINE (X, Y, U, N, A, B, C, D, KОDE, IER) обчислює коефіцієнти кубічного сплайну. Принагідно обчислюються (якщо не задані) значення перших похідних у точках з'єднання. Програма має такі параметри:

X,Y,U вектори значень аргументу, функції та перших похідних (довжини N) ; кількість точок; A,B,C вектори коефіцієнтів кубічного многочлена при змінних x 3 x 2 x (довжини N); D вектор вільних членівкубічного багаточлена (довжини N); KОDE ознака завдання крайових умов:

Значення Сенс

-2	на кінцях сітки задані другі похідні; перед викликом програми і повинні бути занесені відповідно до D (1) і D (N),
-1	на кінцях сітки задані перші похідні; перед викликом програми і повинні бути занесені відповідно до D (1) і D (N),
0	крайові умови не задані, що рівносильно ,
1	заданий періодичний сплайн, тобто і ,
2	вектор похідних значень U заданий повністю;

IER код помилки; цей код для програми SPLINE є транзитним і передається у тому вигляді, як він отриманий у підпрограмі TRIDIG .

Коефіцієнти A (I), B(I), C(I) і D(I) відповідають відрізку від X(I) до X(I+1), I = 1,…,N-1.

Програма SPLINT (X,N,A,B,C,D,Y,M) дозволяє обчислити значення кубічного сплайну на рівномірній сітці із заданим кроком. Параметри програми:

X вектор значень аргументу (довжини N); кількість точок; A,B,C вектори коефіцієнтів кубічного сплайну при змінних x 3 x 2 x (довжини N); D вектор вільних членів кубічного сплайну (довжини N); Y вектор значень кубічного сплайну (довжини M); M кількість вузлів рівномірної сітки на відрізку [X(1), X(N)].

Рис.5.1. Приклад побудови сплайнів для функцій sin x та cos x.

Hа рис.5.1 побудований сплайн для функцій sin x і cos x значення яких визначені на сітці з кроком p/6, але для функції cos x значення перших похідних у вузлах задані рівними нулю. Приклад показує, як, керуючи похідними, можна змінювати форму кривої. Малюнок отримано за допомогою наступної програми:

Програма TDMP (X,Y,N,A,B,C,D,KODE) дозволяє обчислити елементи тридіагональної матриці або доповненої тридіагональної матриці (у разі періодичного сплайну) та вектор вільних членів. Програма має такі параметри:

X,Y вектори значень аргументу та функції (довжини N); кількість точок; A,B,C піддіагональний, діагональний і наддіагональний вектори матриці (довжини N); D вектор вільних членів (довжини N); KОDE ознака завдання крайових умов:

Значення Сенс

-2	на кінцях сітки задані другі похідні; перед викликом програми і повинні бути занесені відповідно до D (1) і D (N),
-1	на кінцях сітки задані перші похідні; перед викликом програми і повинні бути занесені відповідно до D (1) і D (N),
0	крайові умови не задані, що рівносильно ,
1	заданий періодичний сплайн, тобто і .

Програма TRIDIG (U,N,A,B,C,D,KODE,IER) дозволяє знайти рішення системи рівнянь, заданої тридіагональної матрицею або доповненої тридіагональної матриці, методом прогонки. Зокрема, для кубічного сплайну рішенням системи є вектор перших похідних у точках з'єднання. Параметри програми:

U вектор результатів; N кількість рівнянь; A,B,C піддіагональний, діагональний та наддіагональний вектори матриці; D вектор вільних членів; KОDE характеристика матриці:

Значення Сенс

0	тридіагональна матриця,
1	доповнена тридіагональна матриця;

IER код помилки:Значення Сенс

0	помилки немає,
1	якщо B (1) = 0 ,
2	якщо B(J) + C(J) ´A(J-1) = 0 .

Зауваження. Спосіб завдання матриць:

а) тридіагональна матриця (KODE = 0):

задається так:

A (1) = 0 ,A (2) = r 21 ,…,A ( N ) = r n,n-1	- Піддіагональний вектор,
B (1) = r 11 ,B (2) = r 22 ,…,B ( N ) = r n,n	- діагональний вектор;
C (1) = r 12 ,C (2) = r 23 ,…,C ( N ) = 0	- Наддіагональний вектор;

б) доповнена тридіагональна матриця (KODE = 1):