Обобщение одного алгоритма построения рекурентних сплайнов Текст научной статьи по специальности «Математика»

■а о

Запропоновано модифшаци алгоритму Д. О. Силаева побудови згладжуючо-го сплайна рiзних порядтв гладкостi: вид нульового до другого, як спрямоваш на тдвищення стiйкостi цього алгоритму. ОбГрунтовано рекомендаци щодо форми подання полiномiв, як опису-ють ланки сплайтв указаного виду. На обчислювальних прикладах показано можлив^ть узагальнення алгоритму Д. О. Силаева на нерiвномiрнi ытки

Ключовi слова: згладжуючий сплайн, часовий ряд, стштсть алгоритму, обу-мовлетсть матрищ, ресурсоемтсть

алгоритму

□-□

Предложены модификации алгоритма Д. А. Силаева построения сглаживающего сплайна разных порядков гладкости: от нулевого до второго, которые направлены на повышение устойчивости этого алгоритма. Обоснованы рекомендации относительно формы представления полиномов, которые описывают звенья сплайнов указанного вида. На вычислительных примерах показана возможность обобщения алгоритма Д. А. Силаева на неравномерные сетки

Ключевые слова: сглаживающий сплайн, временной ряд, устойчивость алгоритма, обусловленность матрицы, ресурсоемкость алгоритма

■о о

UDC 519.632.4

|DOI: 10.15587/1729-4061.2018.128312

GENERALIZATION OF ONE ALGORITHM FOR CONSTRUCTING RECURRENT SPLINES

G. ^luchen^

Doctor of Technical Sciences, Professor Department of Higher Mathematics and Mathematical Modelling Kherson National Technical University Beryslavske highway, 24, Kherson, Ukraine, 73008 Е-mail: tuluchenko.galina@ukr.net

G. Virchen^

Doctor of Technical Sciences, Associate Professor Department of Descriptive Geometry, Engineering and Computer Graphics National Technical University of Ukraine «Igor Sikorsky Kyiv Polytechnic Institute» Peremohy ave., 37, Kyiv, Ukraine, 03056 E-mail: virchga@gmail.com G. Getu n PhD, Associate Professor* Е-mail: GalinaGetun@ukr.net V. Маrtynov Doctor of Technical Sciences, Associate Professor* Е-mail: arx.martynov@gmail.com M. Tymofieiev PhD, Associate Professor* Е-mail: cotim1951@gmail.com *Department of Architectural Constructions Kyiv National University of instruction and Architecture Povitroflotskyi ave., 31, Kyiv, Ukraine, 03037

1. Introduction

The task on finding an effective analytical description of a time series is commonly in many application areas. The scope of application of interpolation splines is known to have been limited to processing experimental dependences that do not contain measurement errors or contain errors that can be ignored. In the opposite case, the smoothing splines are employed. There are many methods for smoothing. That relates to the fact that the smoothing conditions are chosen to match a particular applied problem. An expert in the subject often struggles to understand the peculiarities of different smoothing methods and to assign the required parameters of smoothing. Therefore, there is a need to develop algorithms of smoothing whose constraints are interpreted so that a user who may lack a proper mathematical training on this subject can understand them.

The task on constructing an optimal smoothing spline implies sorting out all the possible cases of distribution of grid knots among the links of the spline. It is clear that

the resource intensity of algorithms that could fully solve such a problem makes them unsuitable for the application in real-time systems [1]. Therefore, it is a relevant task to build algorithms for constructing the smoothing splines, close to optimal, provided there is the possibility to modify a small number of spline links in proportion to the arrival of new data.

Such algorithms would make it possible to achieve the approximation of an experimental dependence with acceptable quality for much less time than when solving a given problem in classical statement.

2. Literature review and problem statement

The tasks on constructing smoothing splines in various statements have traditionally attracted attention of experts in computational methods and specialists in different subject areas that employ the apparatus of splines [2, 3]. However, these studies mainly concern the development of algorithms

уз

©

for constructing smoothing splines, fulfillment of which is not strictly limited in time. The latter constraint is relevant for many applied contexts of using splines.

One of the important components of the task on constructing smoothing splines is the search for knots, in which spline links are joined. A solution to this auxiliary problem using a genetic algorithm is proposed in paper [4]. The same problem for B-splines is solved in article [5] using the author's two-stage algorithm taking into consideration a curvature of the experimental dependence.

In paper [6], in order to establish coordinates of the knots where spline links are joined, it is proposed first to form their redundant sequence followed by chopping using optimization methods. The framework of this approach includes the results obtained by authors of article [7].

The issue of constructing a robust smoothing B-spline is examined in paper [8]. The above list shows that most of the results were obtained for B-splines and rely on methods that are time-expensive.

We shall consider two algorithms for constructing smoothing splines using which the problems on robustness and splitting a spline into links, close to optimal, can be solved simultaneously. These algorithms have the common feature that implies a partial imposing of the spline's adjacent links one on another. We note that the traditional approach to the construction of splines implies that adjacent links have one common point. Both algorithms are focused on the post-construction of a new link of the spline in proportion to the arrival of new experimental data without changing all the preceding links of a given spline during approximation of a time series.

The first algorithm is described, in particular, in paper

[9]. The algorithm allows the use of non-uniform grids. To describe a spline link, a polynomial with a specific structure is used, the properties of which are examined below. The algorithm is based on the synthesis of ideas about Hermitian interpolation and mean-square approximation. The result of its execution is the smoothing spline of order C0.

The second algorithm is presented, for example, in paper

[10]. The existence of a semi-local spline and its convergence were proven for uniform grids [11]. To describe a link of the spline, the authors used a polynomial with a representation in a power base. The result of the algorithm execution makes it possible to obtain a spline with a preset order of smoothness [12].

Each of the algorithms has both undeniable advantages and certain disadvantages. Thus, a shortcoming of the first algorithm is the low degree of smoothness of the spline obtained. In addition, local estimation of part of the coefficients of a polynomial from the current spline's link for experimental dependences with significant measurement errors may lead to a significant worsening in the approximating properties of the spline in general.

Effectiveness of applying the second algorithm is limited by the requirement to construct links of the spline strictly equal in length. For the plots of experimental sequences that demonstrate the slow nature of changes, it is advisable to apply longer links of the spline. This, in turn, makes it possible to bring down the cost of resources when maintaining a description of the resulting spline and during its application. Therefore, it is expedient to synthesize these algorithms in order to improve their advantages and to eliminate disadvantages.

3. The aim and objectives of the study

The aim of present study is to examine a possibility to generalize the D. A. Silaev algorithm for irregular grids and to strengthen its computational stability. This would make it possible to reduce resource intensity of the examined algorithm when used in monitoring systems.

To accomplish the aim, the following tasks have been set:

- to substantiate the feasibility of transition from representing a polynomial to describe links of the spline in the power base to other types of bases, specifically the bases of Hermite polynomials, the base of N. D. Dicoussar polynomial;

- to show a relation between the N. D. Dicoussar polynomial and other known types of polynomials;

- to examine the impact of the form of polynomial representation to describe links of the spline on the magnitude of the conditionality numbers of stability matrix and the Gram matrix, which are used in the D. A. Silaev algorithm;

- to explore the effectiveness of constraints that are applied to determine the length of the spline's current link in the N. D. Dicoussar algorithm, in order to use them when generalizing on the irregular grids of the D. A. Silaev algorithm;

- using test examples, to explore the resistance of the D. A. Silaev algorithm, and modifications of this algorithm, against the errors in experimental data.

4. Materials and methods of research

4. 1. Relation between the N. D. Dicoussar polynomial and other kinds of polynomials

In cases when the study relates to one link of the spline in order to reduce the number of indexes we shall equate the region for determining this link to section [ xa; xb ].

Paper [9] proposed to describe a link of the cubic smoothing spline by a polynomial of the following form:

PD ( x ) = r d + 6 Q,

(1)

where r=(fa; fj; f0) is the vector of weight coefficients at basis elements; d=(dù d2; d3) is the vector of basis elements in the form of polynomials of the second degree:

d1 =

d2 =

(x-xo)(x-xb) _

( x0- - x a )( xb - xa )'

(x - xo )(x " " x a )

( xb - xo)( xb - xa )

( x - x a )( x - "xb )

( x0 xa )( x0 xb )

Q = ( x - xa Mx - x0 M x - xb )

is the auxiliary "zero" cubic parabola.

The basis elements (in the terminology of paper [9]) (df; d2; d3) are the basis functions in the traditional sense, associated with knots xa, x0 and xb. That is, these functions are equal to 1 in their knot and are equal to 0 in all other knots. Function Q is zero in all three knots, as derived from its formula.

We shall show in advance that the set of functions (di; d2; d3; Q) forms the basis of space of polynomials with degree not higher than the third degree [13]. To this end, we shall construct a matrix from coefficients of the specified functions:

( fa ; fb ; /0; eK1 =

= ( fa; /b; /1; fb ) ' Viî1 .

D =

xn + xu

1

I x„ 11 xu

xu x x x „

I xa x„ 11 x. x„

x xn x x„

a b

xu + x„

ix0 xa)\x0 xb x„xr.xL

I xn x 11 xn xu

( x0 xa )( xb xa 1

(xb -x0 )(xb xa ) 1

(x0 -xa )(x0 -xb )

-(x + x + x, )

Hence, we obtain r 0

' 1 1 1

xa xb x0 0

x2 a x 2 b x 2 b 0

x 3 a x 3 b x 3 b

PD(x)=(/a; /; /0; e).d■ x= = (f; /b; /0; e}yD 1.x.

phi(x) = (/; /b; //; //K1 ■ x,

where

Vh, =

( 1 1 0 0

xa xb 1 1

x2 a x 2 b 2 x 2x

x 3 a x 3 b 3 x2 a 3 xt

(2)

and compute the determinant:

I I 1

(x0 - xa)(x0 - xb)(xb - xa)

It is obvious that functions { d\, d2; d3; Q Q} are linearly independent and form a basis of the space of polynomials with a degree not higher than the third degree.

Find the matrix that is inverse to matrix D and denote it VD:

VD = D-1 =

By using the technology of block matrices for constructing basis functions, which is described in paper [14] (in relation to the construction of basis functions in a finite element method), one can show that polynomial (1) in the matrix form can be recorded as:

V-11 V =

(/; /b; /0; e)=

=( L ; /b ; /a ; / K1 V.

With respect to the structure of matrices V— and VD, we obtain:

2 ■■

Cb - x0 ) (2x0 - 3xa + xb ) 2

1 0

0 1 -

0 0

0 0

2■(xb -xa) (x0 - xa )2 (2x0 + xa - 3xb ) 2 ■( xb - xa )

I x: xn I I xn x„

Hence, we obtain an expression for parameter e:

e = -

( xb - x a )

// +// -

2.( /b -1 )

(3)

Author of paper [9] obtianed an estimate of the magnitude of parameter 8 through the value of function f(x) and its derivatives at the ends of segment [xa; xb]. We show that it can be achieved in another way.

Let us consider the Hermitian polynomial, which hereafter is referred to as the Hermite polynomial of the first kind, defined by their values and the values of the first derivatives at the ends of current segment [xa; xb]. Paper [14] shows that it can be represented in the following matrix form:

which fully coincides with the expression, derived in the N. D. Dicoussar papers.

We note that the value of this parameter does not depend on the localization of knot x0.

It will be shown below that it is appropriate to choose the origin of a local coordinate system for the spline's link at point x0, that is, in this case, x0.=0. The matrix of the N. D. Dicoussar polynomial coefficients (2) then takes the form:

D=

(4)

x

1

0 -

1 0

xb xa ) ■ xa

x

( xb - x a 1

I x, - x„

■ x.

x>, + x„

( xb - x a 1

xu ■ x„

Because formulae (3) and (4) describe the same polynomial in different bases, we can equate the right sides of these formulae:

This allows us to consider the N. D. Dicoussar polynomial as a hierarchical form of the Hermitian polynomial.

4. 2. Approximating properties of the N. D. Dicoussar polynomial

When executing both of the examined algorithms, a system of linear algebraic equations is solved, which is traditional for the least squares method. It is known the matrix of this system coincides with the Gram matrix for basis functions used to represent the approximating

1

polynomial. Therefore, an important characteristic of the basis that is chosen to represent a polynomial is the conditionally number of its Gram matrix. Not reducing the generality of results, we shall compute elements of the Gram matrix for the basis functions of polynomial (1) for segment [xa; xj] = [-1;1j:

(G )=fÎ (d • d ) dx

V-1

(5)

where i, je{1, 2, 3, 4} considering that d4=Q.

We shall tabulate values for a conditionality number of the the Gram matrix from the position of point x0. Paper [9] proposes fixing the positions of point x0 as second knot of approximation at segment [xa; xjj. Our calculations show that this recommendation may be disregared. The dependence of the Gram matrix conditionality number on the position of point x0 (Fig. 1) shows that a shift of the right end of segment xj leads to a rapid growth in the Gram matrix conditionality number (along with the accumulation of calculation errors). That is why knot x0 must be chosen in a spline grid knot, which is the closest to the midpoint of segment [xa; xj]. The results obtained (Table 1) correspond to the point of view on this issue, known from literature, for example, from ref. [15].

Fig. 1. Dependence (5) of conditionality number cond(G) of the Gram matrix G on the localization of point x0 at segment [-1; 1]

Table 1

Estimation of approximating properties of basis functions of the N. D. Dicoussar and Hermitian interpolation polynomials

Characteristic N. D. Dicoussar basis (xo=0) Hermitian basis

Gram matrix determinant 0.0027 0.000042

Gram matrix conditionality number 32.80 292.82

Based on data from Table 1, we conclude that the application of the N. D. Dicoussar interpolation polynomial (1) in approximate calculations demonstrates considerable advantages in comparison with the Hermitian interpolation polynomial because its basis functions possess better approximating properties.

4. 3. Generalization of the D. A. Silaev algorithm

We shall introduce changes to the algorithm described in paper [10]. This relates to the quantity of points, used in the method of least squares, to find coefficients at senior degrees of the polynomial. We also propose changes to the form of a polynomial representation, which describes the current link of the spline.

We shall consistently add knots to the grid at which the problem on the minimization of a functional of least squares is solved until the resulting polynomial deviates from experimental data by less than the specified number 5>0. The final number of knots from the current interval at which this condition is fulfilled will be designated (M+1).

The number of points that finally remains in the composition of the current link of the spline will be designated (m+1).

To assess the impact of the form of representation of a polynomial on the stability of the modified algorithm, we shall use such kinds of polynomials.

When building a spline of smoothness C0, we shall apply:

- polynomials in a power basis;

- the N. D. Dicoussar polynomials of form (1).

When building a spline of smoothness C1, we shall use:

- polynomials in a power basis;

- Hermitian polynomials of the first kind with two knots in which values for the function and its first derivatives are assigned;

- Hermitian polynomials of the second kind with three knots in which values for the function are assigned, an in the first knot - values for the function's derivative.

When building a spline of smoothness C2, we shall use:

- polynomials in a power basis;

- Hermitian polynomials of the second kind with two knots in which values for the function are assigned, an in the first knot - values for the derivatives of the first two orders.

4. 4. Constructing a spline of smoothness Co

The matrices employed in the examined algorithm, when applying a polynomial with a power basis, are described in detail in the paper of its developer [9]. To maintain a unified approach, we shall rewrite basis functions of the N. D. Dicoussar polynomial (1) in a local coordinate system with the origin at point xa=0. Then the other two knots will have the coordinates: xb=Mh, x0=qh, where q=[M/2], [ ] is a whole part from the number, h is the distance between knots. Then the basis functions will be assigned by formulae:

, 1

d. =--

1 qMh

d2 (M - q) Mh2

d (M - q ) qh2

-( x - qh)( x - Mh);

;( x - qh);

•x ( x -Mh);

Q = x (x - Mh)(x - qh).

(6)

We shall write respective matrices when applying the N. D. Dicoussar polynomial (6). Basis function d., which is associated with knot xa, will be used under conditions of smooth joining, and coefficients at functions d2, d3 and Q will be searched for by the method of least squares.

A condition for joining the l-th and (l+1)-th links of the /+1 = B1 ■ A-1 ■ Pl + U ■ fla, (7)

spline with a smoothness of order C0 will be consistently written in the form: where

PD (mh) = PD+1(0);

(d1(mh))( J ) + (d2(mh) d3(mh) Q(mh))-

( fi\ fb

f'

Jo

e'

v у

=( fa'+1 );

A ■( fl )+A

'

fb

f' e'

v у

= ( /«'+1 )-

where

(m - q)(m - M)

A =-'

A =

qM

m (m - q) (m - M) m

M (M - q)' (M - q ) q

m (m - M)(m - q) h3

4 ■ /+4

fb

e

v у

=P'.

Hence, we obtain

fb

f'

o e'

v у

= (Ai)-1 ■(P' -40 ■ /') = (4)-1 ■(P' -40 ■ f-1).

U = 50 - A ■ ■ A0

is the stability matrix.

To examine stability of the algorithm depending on the form of a polynomial representation, we shall tabulate the values:

- of the largest module of eigen number of stability matrix U;

- of the conditionality number of matrix of the least squares method Ai, for different ratios of values of Mto m.

Table 2 will include ratios between M and m, which lead to the lowest of the largest modules of eigenvalues of stability matrices U at a fixed value of M. We consider h=1 during tabulation.

Table 2

Stability parameters of the algorithm for constructing a smoothing spline of order C0

When implementing a least squares method, we shall use a local numbering of knots for index i, assuming xa=X0, X,+i=X,+h.

A matrix of the least squares method takes a traditional form:

( M

X d, (X,) d2 (X,)

¿=0 M

X d, (X,) d3 (X,) ■ f +

¿=0 M

X d, (X, )Q (X,)

0

M

Xd2 (X,)d2 (X,) Xd3 (X,)d2 (X,) XQ(X,)d2 (X,)

S=0 S=0 S=0

M M M

X d 3 (X,) d 2 (X,) X d3 (X,) d3 (X,) XQ (X,) d3 (X,)

S=0 S=0 S=0

M M M

XQ(X,)d2(X,) XQ(X,)d3(X,) XQ(X,)Q(X,)

S=0 S=0 S=0

or according to denotation applied, for example, in paper [12]:

М m min max Я M 1 ' 1 Conditionality number of the matrix of least squares method A1

Polynomial with a power basis N. D. Dicoussar polynomial

4 1 0.05797 32,952.12 18.26

4 3 0.05797 32,952.12 18.26

5 3 0.03306 58,995.13 115.20

6 4 0.02564 1.02105 309.58

7 4 0.05085 1.68105 1,019.02

8 5 0.00000 2.64405 2,157.58

9 2 0.02376 3.99405 5,179.17

10 6 0.01770 5.8M05 9,454.08

f M I yd ( X )

f fb 1 ¿=0 M

/0 = I yd ( X )

leJ ¿=0 M

I yQ( X ) ^ ¿=0

Thus, we obtain a recurrent formula to calculate coefficient /+1 from the next link of the spline through the value of coefficient /la from the current link of the spline, which ensures smoothness of joining the links of the spline of order C0:

Data from Table 2 testify to significant benefits of polynomials in the N. D. Dicoussar form when applying them for the problem on the construction of a smoothing spline. These benefits are ensured by a significant decrease in the conditional-ity number of the matrix of least squares method Ai.

4. 5. Constructing a spline of smoothness С

It is obvious from the geometrical content of coefficients of the N. D. Dicoussar polynomial that the form of this polynomial representation leads to significant complication in the algorithm for constructing an approximating spline with a smoothness order larger than С0. Given the relation between the N. D. Dicoussar polynomials and the Hermitian polynomials, established in chapter 4, we shall explicitly record the matrices used in the D. A. Silaev algorithm for the Hermitian polynomial (4).

We shall also rewrite the basis functions of polynomial (4) in the local coordinate system with origin at point xa=0. In addition, we shall change the numbering order for basis functions in order to simplify the system of equations describing conditions for the smoothness of the spline. Thus, polynomial (4) takes the following form:

PBi( x )=( fa ; //; fb ; // K1 ■ x, where

1

(8) or according to the notation used, for example, in paper [12]:

( fi \

2Mh

( x - Mh)-

/a

i

-A1

The basis functions of polynomial (8) will then be assigned by formulae:

,( /a

Hence, we obtain

/

( /1 )

= Pl.

(10)

1 2 1

U —__o,q (2x + Mh)(x-Mh) ; u2 — , ., x(x

1 M3 h3' 1

u = —

m 2h2

- Mh)2;

( / )

=( 4 r-

pl - A

.( //

(11)

M h

-(2x - 3Mh)- x2; u4 — x2 (x - Mh).

Recurrent formula (7) takes the following form in this case:

Coefficients at basis functions u and u2 will be found from the conditions for the smoothness of joining links of the spline, and the coefficients at functions u3 and u4 will be searched for by the method of least squares.

/l+ J a

( //:

— B1 - A1-1 - Pl + U

f J a

.( /l\

(12)

and ensures smoothness of joining links of the spline of

Condition for joining the i-th and (i+1)-th links of the order C1.

spline with a smoothness of order C1 will be consistently written in the form:

PH(mh) — P^O);

x))X r = mh = №(x))'

We shall consider another form of representation of the Hermitian polynomial with three knots in which values for the function are assigned, in the first knot - also the value of the function's derivative:

x = 0'

u1(mh)

u2(mh)

ph2( x )=( /a ; f; /0; f k1 - x,

where

(13)

u1(x ))

x = mh u3(mh)

2 ( x))

x = mh u4(mh)

( //

3( x ))■

, (u4(x))

; = mh

; — mh

f ■I \ b f /i+1 N J a

( / / )i = ( f )l+1

yv ! \V a !

yH„ —

f 1 0 1 1 N

0 1 qh Mh

0 0 q2h2 M 2h2

y 0 0 q3 h3 M 3h3 ,

M(m+2m)(m-m)2 M(M -m)2

—(M - m) -K(M - 3m)(M - m) hM 31 ; M21 n

^ (3M - 2m) -M-(M - m)

(M - m) —(2M - 3m)

. hM31 ' M21

/ fi \

/a .( //

xb = Mh, x0 = qh, q = [M/2].

The basis functions of polynomial (13) will then be assigned by formulae:

1

1 „2^,21,3

q2M h

( x - Mh)( x - qh)((M + q) x + Mqh);

f fl > Jb f /1+1 ^ a 1

l( // )l J = l( // )l+1J ; —-TT 2 qMh2

Be - f il ^ a ( f / )l +B - f / ( l > b / )l — f /1+1 ^ a ( f / )l+1

l(fa ) J ) J l(/a ) J

x(x - Mh)(x - qh); 1

( x - Mh)- x2;

q2 (M - q ) h3 1

' M2 (M - q) h3

; - qh) -

A matrix of the least squares method takes a traditional form:

/ M

X u ( X ) u3 ( X, ) X u (X, )u3 ( X, )

i=0 i=0 M M

X u ( Xi ) u4 ( Xi ) X u2 (Xi )u, ( X, )

0 i=0 M M

Xu3 (Xi)u3 (Xi) Xus (Xi)u (Xi)

i=0 i=0 M M

Xu3 (Xi)u, (Xi) Xu, (X)u4 (Xi)

( /

/l ( / )

M

X yu ( Xi )

i=0 M

X yiu4 ( Xi )

Coefficients at basis functions Vj and v2 will also be found from the conditions for smoothness of joining links of the spline, and coefficients at functions v3 and v4 will be searched for by the method of least squares.

Condition for joining the i-th and the (i+1)-th links of the spline with a smoothness of order C1 will be consistently written in the following form:

(9)

PH2(mh) — P£(0);

(PH2(x))X x = mh = (PH+1(x))

x = 0'

3

3

2

x.

4

v1(mh)

/

v2(mh)

v1( x))

v

v3( x ))'

V

/

x = mh v3(mh)

x = mh

V2( X))

v

(v4( x ))■

x = mh Vi(mh)

x = mh

fl Ja

f / )l +

fa )

( fl ^ f0 ' fl+1 N J a

V fbj l( f/ )l+1]

PB3(x) = (fa; f/; f//; f„ K1 ■ X, where

'10 0 1 ^ 0 1 0 Mh 0 0 2 M2h2 . 0 0 0 M3h3;

(14)

K =

1 mh

—(m - q)(M - m)(Mm + Mq + mq) ~ q)(M - m)

-1

2 ■( M2 + Mm + m2 )-3m ( M + q)) — (2m ( M + q)-( Mq + 3m

hq2M2

fl Ja

.( f/:

(M - q) q-m

h(M - q)q'

-(M - m) -(2M - 3m)

(M - q) M2

h (M - q) M:

-(3m - 2q )

( fl \ fb ' fl+1 N Ja

i( f )lJ l( fl )l*

l

A

.( f/

-A

/ /l+1 \

V-'b 7

( f/

The basis functions of polynomial (14) will then be assigned by formulae:

1

M 3h3 1

~M2h2 1

( x3 - M3 h3 ); x(x2 - M2h2);

2Mh 1

M3 h3 "

(x - Mh)- x2;

Other calculation formulae are similar to formulae (9)-(12), which is why we shall not write them.

Results of studying the stability of an algorithm depend-

Coefficients at basis functions wi, w2 and w3 will be found from conditions fpr smoothness of joining links of the spline, and coefficient at function w4 will be searched for by

ing on the form of a polynomial representation are given in the method of least squares.

Table 3.

Condition for joining the l-th and the (l+1)-th links of the spline with a smoothness of order C2 will be consistently Tab|e 3 written in the form:

Stability parameters of the algorithm for constructing a smoothing spline of order C

Conditionality number of matrix of the least squares method A1

M m min max X. M 1 ' 1 Polynomial with a power basis Hermitian polynomial (8) Hermitian polynomial (13)

4 2 0.2433 667.74 8.57 2.44

5 4 0.1315 915.45 7.60 2.17

6 5 0.1217 1,219.06 7.59 3.42

7 6 0.1156 1,576.11 8.13 2.98

8 7 0.1129 1,985.63 9.07 4.11

9 7 0.1332 2,447.17 10.31 3.61

10 8 0.1270 2,960.52 11.84 4.63

M3 - m3 mh

M3 M2

3m2 1

M3 h M2 ■

6m

M3 h2

PH 3(mh) = PiH+1(0);

(PH3( x ))/| x = mh = (PH+1( x ))

K(xf x = mh = №(x))/

x = 0' x = 0'

w1(mh)

/

w2(mh)

w3(mh)

(w1( x ))' (w1(x ))//

x = mh x = mh

(w2(x ))/ (w2( x ))//

It is obvious that it is appropriate to use the Hermitian polynomials in the form (13), which, all other conditions being equal, ensure the best conditionality of matrix A1, which is applied in the implementation of the examined algorithm.

4. 6. Constructing a spline of smoothness C2

To construct the smoothing spline of order C2, we shall use the Hermitian polynomial with two knots, in which values for the function and its two derivatives are assigned in the first knot, and in the second knot - only a value for the function:

l

( f ( f/

w4(mh)

(w4( x))/ (w4( x ))//

x= mh x= mh

x = mh x = mh

( fl ) =

(w3(x ))/ (w3( x ))//

x = mh x = mh

W+1

\

( f/ ( fl'

\

2M

■(M -

m)

mh

-V (M2 - 3m2 ) — (2M - 3m)

2

6m

~Mh

2M -1 (M - 3m)

3

m3

fl Ja M3

( f / )l 3m2

( fa ) + M3 h

( f/ )l j 6m

M3 h2

( fl ) =

( f .( fl

3

m

/ \ 1 /1 I Ja / \ /1+1 a

Bo ■ ( // )1 + Bi ■( /1 ) = ( // )1+1

l( // )1 J i^ )l+1,

It is obvious from data in Table 4 that the application of Hermitian polynomials in the form (14) has advantages in comparison with polynomials with a power basis because it ensures better stability of the algorithm to the impact of measurement errors in experimental data.

A matrix for the least squares method takes a traditional form:

X« (xh(x) X« (x)«(x) x« (xu (x)

0 i=0 i=0

M I / M

X «4 ( X ) «4 ( x ) ( /l HI ( x, )

J a ( //: ( //

or according to the notation used, for example, in paper [12]:

/ \

f J a

Ao

( //

- A ■( /1 )=P.

K //

Hence, we obtain

/ c w 1 /1 II Ja

( /1 )=( A r p1 - a0 ■ ( // )1

l( // )1 JJ

Recurrent formula (7) takes the form in this case:

( // ( //'

= B1 ■ A-1 ■ P1 + U ■

( // ( //'

5. Discussion of results on a possibility of the generalization of the examined algorithm for irregular grids

It is clear from the results of theoretical research that when constructing the spline of smoothness C0 it is appropriate to apply polynomials in the N. D. Dicoussar form; that of smoothness C1 - the Hermitian polynomials of second kind (13). To construct the spline of the examined kind of smoothness C2, it is expedient to apply polynomials of higher power, since when applying the third-power polynomials, in order to determine by the method of least squares, there is only one coefficient, the model turns out to be too rigid and the stability of the algorithm is significantly worse than in the two previous cases.

To generate an experimental sequence, we shall construct a function that is a combination of the three Lorentz functions with parameters that are given in Table 5.

y( x ) = X Vj( x ) = X-

Aw2

j j

j=i

j=l Wj +(x - Xj 0 )

(15)

where Aj is the amplitude; Wj is the semi-width; Xj0 is the frequency position of a maximum.

Table 5

Parameters of a test function

j Aj Wj Xj0

1 1 0.3 -2

2 3 0.5 0

3 2 0.2 3

and ensures smoothness for joining the links of spline of order C2.

Results of studying the stability of an algorithm, depending on the form of a polynomial representation, will be given in Table 4. As matrix A1 is composed in this case of one element, its conditionality number for all the bases is unity.

We shall add to the values of function y(x) an error that has a normal distribution law with parameters a=0; o=0.075.

We shall generate a sequence of N=120 in the interval [-3; 5]. Using rule 2o, we shall add points to each spline's link until the polynomial starts to deviate from experimental points by larger than 2o. The final number of points at which this requirement is fulfilled will be designated (M+1).

In the first interval, we shall solve an auxiliary task on approximation with additional constraints that take into consideration the physical essence of the problem. Specifically: the approximating polynomial at the beginning of the segment accepts a positive value and has a positive first derivative at the same point. Under such conditions, the algorithm selected 12 links of the spline with smoothness C0 with the polynomials from the current link of the spline represented in the N. D. Dicoussar form (1) (Fig. 2). As follows from the conditions for the construction of such a spline, its derivative at the points of connection has gaps (Fig. 3).

Table 4

Stability parameters of the algorithm for constructing a smoothing spline of order C2

M m min max X M 1

Polynomial with a power basis Hermitian polynomial (14)

4 2 0.81040 0.65880

5 3 0.89022 0.62324

6 3 0.80712 0.65650

7 4 0.77350 0.63001

8 5 0.80456 0.62865

9 5 0.77790 0.63368

10 6 0.80270 0.62822

Computational experiments to construct the spline of smoothness G\ with the representation of polynomials from current links in the second Hermitian form (13) revealed the following. The best approximation results, in terms of their robustness, are obtained when the imposing depth of the spline's links increases by one point relative to the number of points that are determined based on the estimates of eigenvalues modules of stability matrices. In this case, the algorithm selects 14 links (Fig. 4).

Fig. 2. Results of the approximation with a spline of smoothness Co by the modified algorithm: 1 - experimental dependence; 2 - chart of function (15); 3 - chart of the spline

Fig. 3. Results of the approximation with a spline of smoothness c0 by the modified algorithm: 1 - experimental dependence; 2 — chart of the derivative from spline with gaps at points

Fig. 4. Results of the approximation with a spline of smoothness Ci by the modified algorithm: 1 - experimental dependence; 2 - chart of function (15); 3 - chart of the spline

It is clear that the spline derivative chart in this case is continuous.

Fig. 2-4 show that at segments with a slow change in the examined parameter the algorithm selects longer spline links than at sections with a fast ascending or descending of

experimental dependence. As expected, this makes it possible to reduce the total number of the spline's links compared to their number at uniform splitting. Thus, while ensuring an acceptable quality of approximation, the resulting spline requires less memory for its storing and using.

The disadvantages that remain inherent to the algorithm after its generalization include subjective determining of the acceptable magnitude of deviation in the chart of a polynomial from the current link of the spline from experimental data. Note that in the case of an unknown probabilistic law of distribution (or its parameters), by which measurement errors abide, such a disadvantage is common to all the algorithms for constructing the smoothing splines of such kind.

Prospects for the further research include the extension of algorithm to use, in its composition, polynomials of higher degrees to ensure the second order of smoothness of joining links of the spline with the simultaneous avoidance of the emergence of parasitic oscillations at an increase in the length of experimental dependence. It is also expedient to devise a formal criterion for determining the length of the current link of the spline in order to ensure the property of robustness for a smoothing spline. Further research will also relate to studying the question of alignment of the specified criterion with the criterion of algorithm stability.

6. Conclusions

1. It is shown that in a transition from representing a polynomial in a power basis to the representation in the basis of the N. D. Dicoussar polynomial (when constructing a spline of smoothness С0) or to the representation in the bases of the appropriate Hermitian polynomials (in the construction of a spline of smoothness С1 and С2), there is a significant reduction in the conditionality number of the Gram matrix for the basis functions of the recommended forms of a polynomial representation. This ensures better approximating properties of the examined spline.

2. It is shown that the N. D. Dicoussar polynomial is a special case of the Hermitian polynomial, which is represented in a hierarchical form.

3. It was established that the form of a polynomial representation affects the magnitude of eigenvalues of the stability matrix, starting at the second order of smoothness while joining links of the spline. Such an impact is missing at the lower orders of smoothness of joining the links of the spline.

4. In all the test examples, investigated by Authors, the use of constraints, suggested by N. D. Dicoussar for preliminary determining the length of a link of the spline, made it possible to obtain a smoothing spline, based on the generalized D. A. Silaev algorithm, of acceptable quality, with links of different lengths.

5. An analysis of the quality of the obtained solutions to test problems revealed that in the processing of experimental dependences with significant measurement errors, there may occur the need to reduce the final lengths of spline's links compared to the lengths that are determined by the D. A. Silaev rule. Therefore, in the presence of significant errors in experimental data, it is appropriate to run a preliminary analysis employing the methods of mathematical statistics. This would make it possible, when executing the algorithm, to pay special attention to ensuring its stability.

References

1. Wang Y. Smoothing Splines: Methods and Applications. New York, 2011. 384 p. doi: 10.1201/b10954

2. Smoothing spline analysis of variance models: A new tool for the analysis of cyclic biomechanical data / Helwig N. E., Shorter K. A., Ma P., Hsiao-Wecksler E. T. // Journal of Biomechanics. 2016. Vol. 49, Issue 14. P. 3216-3222. doi: 10.1016/j.jbiomech.2016.07.035

3. Liu Z., Guo W. Data Driven Adaptive Spline Smoothing. Statistica Sinica. 2010. Issue 20. P. 1143-1163.

4. Knot selection for least-squares and penalized splines / Spiriti S., Eubank R., Smith P. W., Young D. // Journal of Statistical Computation and Simulation. 2013. Vol. 83, Issue 6. P. 1020-1036. doi: 10.1080/00949655.2011.647317

5. Yuan Y., Chen N., Zhou S. Adaptive B-spline knot selection using multi-resolution basis set // IIE Transactions. 2013. Vol. 45, Issue 12. P. 1263-1277. doi: 10.1080/0740817x.2012.726758

6. Knot calculation for spline fitting via sparse optimization / Kang H., Chen F., Li Y., Deng J., Yang Z. // Computer-Aided Design. 2015. Vol. 58. P. 179-188. doi: 10.1016/j.cad.2014.08.022

7. Brandt C., Seidel H.-P., Hildebrandt K. Optimal Spline Approximation via 10-Minimization // Computer Graphics Forum. 2015. Vol. 34, Issue 2. P. 617-626. doi: 10.1111/cgf.12589

8. Bureick J., Alkhatib H., Neumann I. Robust Spatial Approximation of Laser Scanner Point Clouds by Means of Free-form Curve Approaches in Deformation Analysis // Journal of Applied Geodesy. 2016. Vol. 10, Issue 1. doi: 10.1515/jag-2015-0020

9. Dicoussar N. D. Optimizaciya resheniya v zadachah kusochno-polinomial'noy approksimacii. Dubna, 2016. 14 p.

10. Silaev D. A. Polulokal'niye sglazhivayushchie splayny // Trudy seminara im. I. G. Petrovskogo. 2013. Issue 29. P. 443-454.

11. Silaev D. A. Polulokal'nye sglazhivayushchie S-splayny // Komp'yuternye issledovaniya i modelirovanie. 2010. Vol. 2, Issue 4. P. 349-357.

12. Polulokal'nye sglazhivayushchie splayny klassa S1 / Silaev D. A. et. al. // Trudy seminara imeni I. G. Petrovskogo. 2007. Issue 26. P. 347-367.

13. Liu X., Yang Q., Jin H. New Representations of the Group Inverse of 2x2 Block Matrices // Journal of Applied Mathematics. 2013. Vol. 2013. P. 1-10. doi: 10.1155/2013/247028

14. Pinezhaninov F., Pinezhaninov P. Bazisnye funkcii dlya konechnyh elementov. URL: http://old.exponenta.ru/soft/mathemat/ pinega/a1/a1.asp

15. Van Manh P. Hermite interpolation with symmetric polynomials // Numerical Algorithms. 2017. Vol. 76, Issue 3. P. 709-725. doi: 10.1007/s11075-017-0278-0