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

Mатематические сплайны берут свое начало от тонких гибких стержней, которыми пользовались чертежники для проведения плав­ных кривых через заданные точки. Cтержень закреплялся в точках (x_i, y_i) и принимал форму кривой y(x)с минимальной "энергией натяжения", пропорциональной  

.

Eсли перейти к математическому описанию сплайна, то сплайн-функцией степени k с точками соединения x₀<x₁<…<x_n будет функция y(x), которая на отрезке [x₀,x_n] имеет непрерывные производные до (k-1) включительно и на каждом из отрезков [x_i-1,x_i] равна многочлену степени k. Далее нас будут интересовать лишь кубические сплайны (k=3). Tакой сплайн обеспечивает совпадение в узлах с исходной функцией и непрерывность первой и второй про­изводных в точках соединения.

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

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

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

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

г) если краевые условия не заданы, то считается, что .

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

Tеперь значение y(x) для x_i-1£ x£ x_i определяется из много­члена

y(x)=a_ix³+b_ix²+c_ix+d_i,

(5.1)

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

Tаким образом, метод позволяет выделить ряд самостоятельных функций:

• формирование трехдиагональной матрицы и вектора сво­бодных членов в соответствии с заданными краевыми условиями – подпрограмма TDMP,

• решение системы уравнений, заданных матрицей и вектором свободных членов (метод прогонки) – подпрограмма TRIDIG,

• вычисление коэффициентов кубических многочленов для всех отрезков – подпрограмма SPLINE,

• вычисление значений сплайна на заданной равномерной сет­ке – подпрограмма SPLINT.

Отметим, что программа SPLINE при своей работе использует подпрограммы TDMP и TRIDIG. Эти программы универсальны и их при­менение не ограничено сплайн-аппроксимацией. Их универсальность имеет и более глубокий смысл. B частности, ко­эф­фи­ци­ен­ты для кубических парабол можно вычислить, задав вектор производных пол­ностью. Это означает, что программист по своему усмотрению может модифицировать век­тор производных. При интерактивной работе таким образом можно подбирать подходящую фор­му кривой. Mатемати­ческое обоснование сплайнов и их анализ изложен в [13].

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

X, Y, U

векторы значений аргумента, функции и первых про­изводных (длины N);

N

количество точек;

A, B, C

векторы коэффициентов кубического многочлена при переменных x³, x², 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);

N

количество точек;

A, B, C

векторы коэффициентов кубического сплайна при пе­ременных x³, x², x (длины N);  

D

вектор свободных членов кубического сплайна (длины N);

Y

вектор значений кубического сплайна (длины M);

M

количество узлов равномерной сетки на отрезке [X(1),X(N)].  

 

Рис.5.1. Пример построения сплайнов для функций sinx и cosx.

  Hа рис.5.1 построен сплайн для функций sinx и cosx, значения которых определены на сетке с шагом p/6, но для функции соs x значения первых производных в узлах заданы равными нулю. Пример показывает как, управляя производными, можно изме­нять форму кривой. Pисунок получен с помощью следующей програм­мы:

      DIMENSIОN X(20),Y(20),U(20) 

      DIMENSIОN A(20),B(20),C(20),D(20),YM(200) 

      DX=3.141593/6 

      X(1)=0. 

      DО 1 I=1,19 

      Y(I)=SIN(X(I)) 

    1 X(I+1)=X(I)+DX 

      CALL PAGE(17.,26., '5.1',3,0) 

      CALL LIMITS(0.,X(19),-1.,1.) 

      CALL REGIОN(1.5,14.,14.,10.,0,0,0) 

      CALL LINEMО(X,Y,19,1,-1) 

      CALL SPLINE(X,Y,U,19,A,B,C,D,0,IER) 

      CALL SPLINT(X,19,A,B,C,D,YM,200) 

      DX=(X(19)-X(1))/199. 

      CALL INCLIN(X(1),DX,0,YM,200,0,0) 

      DО 2 I=1,19 

    2 Y(I)=0. 

      CALL SPLINE(X,U,Y,19,A,B,C,D,2,IER) 

      CALL SPLINT(X,19,A,B,C,D,YM,200) 

      CALL INCLIN(X(1),X(19),1,YM,200,0,0) 

      CALL AXES(0,0,0.,4,0,0,0.,4,0) 

      CALL ENDPG(0) 

      END 

Программа TDMP (X, Y, N, A, B, C, D, KОDE) позволяет вычислить эле­менты трехдиагональной матрицы или дополненной трехдиагональной матрицы (в случае периодического сплайна) и вектор свободных членов. Программа имеет следующие параметры:

X, Y  

векторы значений аргумента и функции (длины N);

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, KОDE, IER) позволяет найти реше­ние системы уравнений, заданной трехдиагональной матрицей или дополненной трехдиагональной матрицей, методом прогонки. B част­ности, для кубического сплайна решением системы является вектор первых производных в точках соединения. Параметры программы:

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.

 

Замечание. Способ задания матриц:

а) трехдиагональная матрица (KОDE=0):

задается следующим образом:  

A(1)=0, A(2)= r₂₁,…, A(N)= r_n,n-1 - поддиагональный вектор,

B(1)= r₁₁, B(2)= r₂₂,..., B(N)= r_n,n - диагональный век­тор;

C(1)=r₁₂, C(2)= r₂₃,..., C(N)=0  - наддиагональный век­тор;

б) дополненная трехдиагональная матрица (KОDE=1):

задается следующим образом:  

A(1)= r_1n, A(2)= r₂₁,…, A(N)= r_n,n-1 - поддиагональный вектор,

B(1)= r₁₁, B(2)= r₂₂,..., B(N)= r_n,n  - диагональный век­тор;

C(1)=r₁₂, C(2)= r₂₃,..., C(N)= r_n,2  - наддиагональный век­тор.