Stable numerical evaluation of multi-degree B-splines
aa r X i v : . [ m a t h . NA ] F e b Stable numerical evaluation of multi-degree B-splines
Carolina Vittoria Beccari a, ∗ , Giulio Casciola a a Department of Mathematics, University of Bologna, Piazza di Porta San Donato 5, 40126 Bologna, Italy
Abstract
Multi-degree splines are piecewise polynomial functions having sections of di ff erent degrees. They o ff er significantadvantages over the classical uniform-degree framework, as they allow for modeling complex geometries with fewerdegrees of freedom and, at the same time, for a more e ffi cient engineering analysis. Moreover they possess a setof basis functions with similar properties to standard B-splines. In this paper we develop an algorithm for e ffi cientevaluation of multi-degree B-splines, which, unlike previous approaches, is numerically stable. The proposed methodconsists in explicitly constructing a mapping between a known basis and the multi degree B-spline basis of the space ofinterest, exploiting the fact that the two bases are related by a sequence of knot insertion and / or degree elevation stepsand performing only numerically stable operations. In addition to theoretically justifying the stability of the algo-rithm, we will illustrate its performance through numerical experiments that will serve us to demonstrate its excellentbehavior in comparison with existing methods, which, in some cases, su ff er from apparent numerical problems. Keywords:
Multi-degree spline, B-spline basis, matrix representation, stable evaluation, algorithmic computation,Greville abscissæ
1. Introduction
Spline functions are the foundation of numerous results and methods of approximation theory and, nowadays,are an integral part of geometric modeling and computational engineering analysis systems. Classically, a univariatespline is a piecewise function defined on a partition of a real interval [ a , b ], where each piece belongs to the space ofalgebraic polynomials of degree less than or equal to d > C d − . The success of splines is largely due to the fact that they possess a B-spline basis , namely a normalized, totallypositive basis of compactly supported functions [1, 2]. Besides the elegance of the theoretical framework, this basis hasexcellent properties from the computational point of view, both for its good conditioning, and because its evaluationcan be carried out via an e ffi cient and numerically stable algorithm, the well-known Cox-de Boor recurrence scheme[3, 4]. These classical splines will be hereinafter referred to also as conventional splines .As the name suggests, multi-degree splines ( MD-splines , for short) are a generalization of conventional splineswhere each piece can have a di ff erent degree. They are a natural and extremely powerful extension of the classicalframework, which allows for modeling complex geometries with fewer control points and at the same time leadsto more e ffi cient engineering analysis [5]. Although the concept of multi-degree splines is long-standing [6, 7], forseveral years the interest in these spaces has been mostly theoretical. Only recently, in fact, has it been understoodhow to build a set of functions analogue to the B-spline basis, dubbed MDB-spline basis (or simply
MDB-splines ).This has opened up the possibility of easily integrating multi-degree splines into current computing systems and hasmade them a real full-fledged extension of conventional splines. Multivariate versions of the multi-degree concepthave been as well devised [5, 8, 9].The first approaches for constructing an MDB-spline basis can be traced back to the work by Shen and Wang[10, 11] and are subject to constraints on the continuity attained at the joins. More importantly, they rely on integral ∗ Corresponding author.
Email addresses: [email protected] (Carolina Vittoria Beccari), [email protected] (Giulio Casciola) ff erent degree are joined with continuity at most C [13, 14].Methods capable of dealing with multi-degree spaces with arbitrary structure stem from a common basic idea,which is to map a set of known (or easily computable) functions into the basis of interest. They can be traced backto two di ff erent approaches. The first consists in determining the basis functions by interpolation, exploiting thefact that, under suitable assumptions, Hermite interpolation problems are unisolvent in MD-spline spaces [15]. Thisinvolves solving a number of (small) linear systems that represent the continuity conditions at the joins. Followingthis approach, in [16] normalized MDB-splines are expressed as combinations of transition functions (a notion earlierintroduced in the context of local spline interpolation [17, 18]), which allows to e ffi ciently compute their expansionwith respect to the collection of the Bernstein bases relative to the breakpoint intervals. In [19] the same idea wasused to deal with splines whose pieces are drawn from Extended Chebyshev spaces [2, 20], a powerful and versatileextension of algebraic polynomials.The second approach stems from the idea of calculating in an explicit way, namely without having to solve anylinear system, a matrix operator M that specifies the mapping between a known basis (or collection of bases) andthe set of MDB-splines. The resulting matrix representation N = M N provides a way to evaluate the MDB-splinebasis functions, which are the elements of vector N , as a combination of easily computable functions, which are theelements of vector N . This avenue was firstly pursued in [5], which underpinning idea is to gather the continuityconstraints between spline pieces in a matrix and then calculate its null space by a recursive procedure. The follow-uppaper [21] proves that the output of the algorithm is exactly the entire set of MDB-splines, whereas implementationdetails are given in [22] and a Chebyshevian extension of the construction is presented in [23]. In these series ofpapers, the vector N is composed of a collection of local bases, which can be, in particular, either the Bernstein basesrelative to the breakpoint intervals or conventional B-spline bases, each one relative to a sequence of intervals of equaldegree. In both cases the functions in N are discontinuous.In [24] it is observed that the functions in N and N are related through a sequence of successive knot insertions(the same observation was made in the context of Chebyshevian splines in [23]) and that matrix M can be computed byinverting these steps, giving rise to a process called reverse knot insertion (RKI) . This realization allowed the authorsof [24] not only to generate the matrix representation [5] in a more direct and intuitive way, but more generallyto derive a matrix representation with respect to any set of functions N which are related to the basis N throughthe aforementioned knot-insertion structure. It is shown, in particular, that to minimize the number of performedoperations, it is convenient to start from a basis N composed of conventional B-spline functions connected with C continuity. Such functions are the MDB-spline basis of a piecewise conventional spline space, referred to as a C MDB-spline space, and, as such can be evaluated by known techniques. The same paper also deals with how togenerate the matrix representation when N is the conventional B-spline basis having maximum (over all intervals)degree and same continuities as the basis N . In this case the procedure consists in inverting the sequence of localdegree elevations connecting the target space and the maximum-degree conventional space and is therefore called reverse degree elevation (RDE) . It is also discussed how it is possible to combine successive reverse knot insertion andreverse degree elevation steps in a unique algorithm. This algorithm hence allows to construct a matrix representationin the most general case, that is under the sole assumption that the space spanned by N contains the space spannedby N .However, both methods [5] and [24] have a weakness, which is that they require to calculate higher-order deriva-tives of B-spline (or Bernstein) basis functions (the order of the derivatives to be calculated corresponds to the maxi-mum continuity or maximum degree to be handled). The evaluation of these derivatives can be carried out in a stableway [25]. However, not only is it a price to pay in terms of computational cost, but also, and above all, it can leadto the numerical instability of the algorithm using them. In fact, B-spline derivatives may easily become very largenumbers for high di ff erentiation order and / or very nonuniform partitions, hence the arithmetic operations carried outwith them are potentially risky. The actual occurrence of instability phenomena was observed in [24], where it issuggested to work with a compensated version of the algorithm in order to improve on accuracy (compensation is astandard technique, see, e.g., [26]). 2he main contribution of this paper is a new, stable algorithm that provides the matrix representation of any MDB-spline basis. The algorithm exploits a suitable reformulation of the reverse knot insertion and reverse degree elevationprocesses, thanks to which only numerically stable operations are performed and, in particular, no derivative needs tobe evaluated. We will show how a natural way to arrive at this reformulation is to pass through the concept of Grevilleabscissæ. The Greville abscissæ are defined, similarly as in the case of conventional splines, as the coe ffi cients ofthe identity function in the MDB-spline basis. As is well known, they are essential in various applications rangingfrom approximation and interpolation to isogeometric analysis. We will show that these abscissæ can be obtained byintegrating the MDB-spline basis of the corresponding derivative space.All in all, the resulting stable algorithm has the form of a triangular scheme. As such it has quadratic computationalcomplexity (vs. the linear complexity of previous algorithms), nevertheless the computation time is negligible inpractical situations. An appealing feature of the approach lies in its deep relationship with the usual spline toolsof knot insertion and degree elevation, which use automatically ensures the correctness of the set of MDB-splinesprovided as output. This is an important di ff erence with respect to the approach in [5], where it is required to prove aposteriori that the generated functions are MDB-splines.For ease of presentation and in the interest of clarity, we will develop in all details the so-called RKI Algorithm ,corresponding to the case where N and N are related by iterated reverse knot insertions. The circumstance in whichthe two basis vectors are related by degree elevation (giving rise to the RDE Algorithm ) can be addressed by similargeneral principles and will be discussed more briefly in the last part of the paper. Finally, we will illustrate how it ispossible to mix RDE and RKI steps, like in [24], in such a way as to be able to choose the initial basis vector N thatwill entail the least number of operations and thereby improve the e ffi ciency of the computation. The latter algorithmbuilds on the previous two and due to space constraints we will limit ourselves to providing a quick sketch of theprocedure.The remainder of the paper is organized a follows. Section 2 collects the necessary notions and results on multi-degree splines and their matrix representation. Section 3 presents the new algorithm, in particular showing how tocompute the Greville abscissæ, using these abscissæ to reformulate the reverse knot insertion process without resortingto MDB-spline derivatives and finally introducing the triangular scheme. Section 4 is concerned with the computa-tional / inplementation aspects of the procedure and presents a practical example of its application. The numericalstability of the method is discussed theoretically in section 5, while subsection 5.1 proposes a series of numerical ex-periments which, in addition to confirming the theoretical predictions, highlight the potential inaccuracy of previousmethods. Section 6 illustrates how to derive a matrix representation in terms of the conventional B-spline basis ofmaximum degree. Conclusions are drawn in section 7.
2. Background and basic notions
In this section we gather the notions and results on multi-degree splines of interest for this paper.
Throughout the paper we will deal with piecewise functions, with pieces drawn from polynomial spaces whosedimensions are allowed to change from interval to interval. Spaces of such functions are defined as follows.
Definition 1 (Multi degree spline space) . Let [ a , b ] be a closed bounded real interval, X = { x i } qi = be a partition of[ a , b ] s.t. a ≕ x < x < . . . < x q < x q + ≔ b and d = ( d , . . . , d q ) be a vector of nonnegative integers. Letalso K = ( k , . . . , k q ) be a vector of nonnegative integers such that k i min { d i − , d i } . The corresponding space of multi-degree splines (MD-splines, for short) is the set of functions S ( P d , X , K ) ≔ n f (cid:12)(cid:12)(cid:12) there exist p i ∈ P d i , i = , . . . , q , such that:i) f ( x ) = p i ( x ) for x ∈ [ x i , x i + ] , i = , . . . , q ;ii) D ℓ p i − ( x i ) = D ℓ p i ( x i ) for ℓ = , . . . , k i , i = , . . . , q o , where P d is the space of algebraic polynomials of degree at most d .3ote that the above definition returns a conventional spline space in the particular case where d is a vector withconstant entries, that is d = · · · = d q , and therefore conventional splines can be seen as a subclass of MD-splines.A space S ≔ S ( P d , K , ∆ ) has dimension K ≔ dim( S ) = d + + P qi = ( d i − k i ) . Moreover, as shown in [16], it possesses a B-spline-type basis, dubbed
MDB-spline basis (or
MDB-splines ),sharing many properties with conventional B-splines. For defining this basis we shall introduce two partitions s and t as follows: s ≔ { s j } Kj = ≔ { a , . . . , a | {z } d + , x , . . . , x | {z } d − k times , . . . , x q , . . . , x q | {z } d q − k q times } , (1)and t ≔ { t j } Kj = ≔ { x , . . . , x | {z } d − k times , . . . , x q , . . . , x q | {z } d q − − k q times , b , . . . , b | {z } d q + } . (2)We call s and t the left and right extended partition , respectively, associated with the MD-spline space.Denoted m ≔ max i { d i } , the MDB-spline basis functions N , m , . . . , N K , m are defined recursively over s and t . Therecurrence consists in constructing, for n = , . . . , m , a sequence of functions N i , n , i = m + − n , . . . , K , where N i , n issupported on [ s i , t i − m + n ] and is determined on each nontrivial interval [ x j , x j + ) ⊂ [ s i , t i − m + n ] by the following integralrelation [16]: N i , n ( x ) ≔ , x j x < x j + n = m − d j , R x −∞ (cid:2) δ i , n − N i , n − ( u ) − δ i + , n − N i + , n − ( u ) (cid:3) du , n > m − d j , , otherwise , (3)where δ i , n ≔ Z + ∞−∞ N i , n ( x ) dx ! − . (4)In (3) we assume undefined N i , n functions to be zero and, in this case, we set Z x −∞ δ i , n N i , n ( u ) du ≔ ( , x < s i , , x > s i . The K functions generated by the above integral formulation possess analogous characterizing properties as con-ventional B-splines, that is:i) Compact support : N i , m ( x ) = x < [ s i , t i ];ii) Positivity : N i , m ( x ) > x ∈ ( s i , t i );iii) End point property : N i , m vanishes exactly d ps i − max { j > | s i = s i + j } times at s i and d pt i − − max { j > | t i − j = t i } times at t i , where ps i and pt i are s.t. x ps i = s i and x pt i = t i ;iv) Partition of unity : X i N i , m ( x ) = ∀ x ∈ [ a , b ].Moreover, the above properties i), iii) and iv) also warrant the uniqueness of the MDB-spline basis (see [16] forfurther details). We refer to previous papers for illustrations of MDB-spline basis functions, see e.g. [5, 16, 22, 24]. Multi-degree splines
In the remainder of the paper we will often rely on MD-spline spaces whose elements are conventional splinefunctions connected with C continuity, which we refer to as C MD-splines .These spaces are convenient tools to work with, in that well-established methods can easily be adapted to deal withthem. In particular a generalization of Cox de-Boor recurrence formula [3, 4], the main method for evaluating con-ventional B-splines, is given in [24, Proposition 4] along with a recurrence relation for the computation of derivatives([24, Proposition 5]). 4he integrals of C MDB-splines can as well be e ffi ciently evaluated resorting to existing results. To this aim, itsu ffi ces to observe that a C MD-spline N i , m can be of only one of the following two types: either it is a conventionalB-spline, which means that within its support each piece has same degree, or it consists of only two non-trivial piecesof di ff erent degrees connected with C continuity. In both cases, recalling that the integral of a conventional degree- d B-spline basis function N i , d is equal to the ratio between the support width of N i , d and d +
1, its integral is easilycomputed as: Z x pti x psi N i , m dx = pt i − X j = ps i x j + − x j d j + , (5)with ps i and pt i defined as in iii). Definition 2 (Associated C MD-spline space) . The C MD-spline space associated with a multi-degree spline space S ( P d , X , K ) as in Definition 1 is the unique MD-spline space S ( P d , X , K ) having same breakpoint sequence anddegree vector as S ( P d , X , K ) and vector of continuities K = ( k , . . . , k q ) such that k i = d i − , d i and k i = k i otherwise. If we take an MD-spline space (including, possibly, a conventional spline space) S such that S ⊂ S , we canwrite the relationship between the respective MDB-spline bases in the form N = M N , (6)where N ≔ ( N , . . . , N K ) and N = (cid:16) N , . . . , N K (cid:17) are the vectors containing the basis functions of S and S , respec-tively, and M is a linear operator of size K × K . We refer to the above as the matrix representation of N relative to N . Knowing M and N , one can thus use (6) to evaluate N .One way to compute matrix M is to choose S to be the C MD-spline space associated with S and rely on iteratedapplication of the following result, referred to as Reverse Knot Insertion , RKI for short [24, Proposition 6].
Proposition 1 (Reverse knot insertion) . Let
S ≡ S ( P d , X , K ) and b S ≡ S ( P d , X , b K ) be MD-spline spaces withsame breakpoint sequence and degrees. Let also the respective continuity vectors be K = ( k , . . . , k j , . . . , k q ) and b K = ( k , . . . , k j − , . . . , k q ) , for j ∈ { , . . . , q } . Then, S ⊂ b S and the corresponding MDB-spline bases { N i , m } Ki = and { ˆ N i , m } K + i = are related through N i , m = α i ˆ N i , m + (1 − α i + ) ˆ N i + , m , i = , . . . , K , (7) where, being ℓ the index of the element of s such that s ℓ x j < min( s ℓ + , b ) , the coe ffi cients α i are such that α i = , i = , . . . , ℓ − d j , ∈ ]0 , , i = ℓ − d j + , . . . , ℓ − d j + k j , = , i = ℓ − d j + k j + , . . . , K + . (8) Moreover, the coe ffi cients α i , i = ℓ − d j + , . . . , ℓ − d j + k j , can be computed from the MDB-splines { ˆ N i , m } K + i = by therelation α i = + α i − D k j − ˆ N i − , m | x j − D k j + ˆ N i − , m | x j D k j − ˆ N i , m | x j − D k j + ˆ N i , m | x j . (9) Remark 1.
Knot insertion is a well-established tool for conventional splines and was generalized to multi-degreesplines in [16]. Classically, using knot insertion we pass from the representation in a space S to that in a space b S inwhich the continuity at a breakpoint is decreased. The word reverse refers to the fact that we go from space b S to S byincreasing the continuity at a breakpoint. Moreover, under the assumptions of Proposition 1, classical knot insertionwould mean to compute the coe ffi cients (8) knowing the MDB-spline bases of both spaces S and b S . Conversely, inreverse knot insertion the coe ffi cients (9) only depend on the basis { ˆ N i , m } of b S and can hence be used to compute thebasis { N i , m } of S .
5o derive the matrix representation (6), take a sequence of MD-spline spaces, all defined on [ a , b ], sharing samebreakpoint sequence X and degree vector d , such that S ≔ S g ⊂ · · · ⊂ S ⊂ S . (10)Furthermore suppose that each space S r , r = , . . . , g , is obtained from S r − increasing by one the continuity at abreakpoint in such a way that K r ≔ dim( S r ) = dim( S r − ) −
1. Relation (7) to pass from S r − to S r can be written inthe matrix form N r = A r N r − , with a bidiagonal matrix A r of size K r × ( K r +
1) having on each row the coe ffi cients α ri and (1 − α ri + ) relating the bases N r and N r − as in (7), (8). Iterating the procedure for as many multiplicities andknots as needed, one obtains the matrix M in (6), which is the product of all matrices A r , that is M = A g · · · A .Being based on repeated reverse knot insertions, the above procedure for the construction of the representationmatrix in (6) is called RKI Algorithm [24]. In order for the RKI Algorithm – that is the matrix representation itproduces – to be an e ffi cient tool for evaluating the target MDB-spline basis N , the vector N should be known (ormore precisely computable through established methods). This is the reason why we have chosen it to contain thebasis functions of the C MD-spline space associated with S . The described procedure, however, simply requires that S can be generated from S by repeated reverse knot insertions. It is therefore possible, and sometimes desired, tochoose S in a di ff erent way, see Remark 2 for further details.It shall be noted that, according to equation (9), the described method requires calculating the derivatives ofMDB-spline basis functions at each breakpoint x j to be processed, up to the target continuity k j . As pointed out in[24], B-spline derivatives, and MDB-spline derivatives likewise, may become very large numbers when the order ofdi ff erentitation increases, causing potential numerical issues for very large degrees and / or nonuniform partitions. Tocircumvent this criticality, in that paper it is proposed to work with a compensated version of the algorithm [26].Although this strategy is able to improve the accuracy of the results, the presence of high order derivatives remains anintrinsic feature of the existing algorithms which would be far preferable to avoid. Remark 2.
Our assumption that S be a C MD-spline space is merely dictated by simplicity and conciseness ofpresentation. However, with a view to choosing S in such a way that its basis can be evaluated by known methods,one can as well consider alternative initial spaces, since the general algorithm proposed in [24] is based on the solerequirement that S be defined on [ a , b ] , have same breakpoint sequence as S and S ⊂ S . In particular we may wantto take S to be a piecewise space and its basis vector N to be the collection of the Bernstein bases relative to thebreakpoint intervals, or alternatively a piecewise conventional B-spline basis on a sequence of abutting intervals withequal degree (in the latter situation, the generated matrix M is the H-operator proposed in [5]). The ideas presentedin this work can easily be adapted to both of these situations.Lastly, we may want S to be a conventional spline space of degree m ≔ max i { d i } , whose basis N will be aconventional degree-m B-spline basis. This situation cannot be addressed by reverse knot insertions, but by a similarapproach, as well described in [24], based on reverse degree elevation and will briefly be discussed in section 6.
3. The novel RKI algorithm
In the following we present the general ideas our new algorithm is based on. To formalize our method we willneed to use spaces spanned by derivatives of MD-splines, that are defined as follows.
Definition 3 (Spline space of derivatives) . We denote by D r S ≔ { D r + f | f ∈ S} the function space whose elementsare r th right derivatives of functions in a multi degree spline space S . There follows that D r S has dimension K − r ,where K is the dimension of S .Throughout the paper, for brevity, we simply write D S in place of D S . We will be concerned with spaces D r S , with r ranging from 0 up to max i { d i } , the latter corresponding to the number of levels in the recurrence (3).Therefore D r S may contain functions discontinuous at breakpoints (the right derivative in the definition accountingfor possible discontinuities) and, rather than a “single” MD-spline space, it should be regarded as “piecewise” space,whose functions are MD-splines defined on abutting intervals and possibly discontinuous at breakpoints. In particular,we associate with D r S the vectors of degrees ( d − r , . . . , d q − r ) and continuities ( k − r , . . . , k q − r ), which may benegative integers, in such a way that K ( r ) ≔ dim( D r S ) = K − r = d − r + + P qi = ( d i − k i ), with the convention that the6estriction of D r S to [ x i , x i + ] be the zero function in case d i − r <
0. We can as well define a piecewise MDB-splinebasis N ( r )1 , . . . , N ( r ) K − r spanning D r S , which will be relative to the left and right extended partitions s ( r ) and t ( r ) obtainedby replacing d i and k i with d i − r and k i − r in (1), (2). These correspond to the functions generated by the the integralrecurrence relation (3) for n = m − r and defined on the partitions s and t relative to S in such a way that: N ( r )1 ( x ) = N r + , m − r ( x ) , . . . , N ( r ) K − r ( x ) = N K , m − r ( x ) , ∀ x ∈ [ a , b ] . (11)Furthermore, we can generalize to MD-spline spaces the classical notion of Greville abscissæ , that are the coe ffi -cients of the expansion of the function f ( x ) = x in the B-spline basis of a conventional spline space containing linearfunctions. For a multi-degree spline space D r S the Greville abscissæ relative to the basis { N ( r ) i , m } are defined to be thecoe ffi cients ξ ( r ) i , i = , . . . , K ( r ) , such that P K ( r ) i = ξ ( r ) i N ( r ) i , m ( x ) = x , ∀ x ∈ [ x j , x j + ], j = , . . . , q , such that d j > D S . Asimilar result was proved in [27, Theorem 15] for Chebyshevian splines with all section spaces of the same dimension. Proposition 2 (Computation of Greville abscissae) . Let S be a K-dimensio- nal MD-spline space. Then the Grevilleabscissaæ ξ , . . . , ξ K with respect to the MDB-spline basis N . . . , N K of S are given by: ξ i = a + i − X j = Z ba N (1) j ( x ) dx , i = , . . . , K , (12) where N (1) j , j = , . . . , K − , is the MDB-spline basis of the derivative space D S .Proof. The properties of the MDB-spline basis { N j } entail that ξ = a and ξ K = b . Moreover, from the partition ofunity property and Abel’s lemma we obtain: K X i = ξ i N i , m ( x ) = ξ K X i = N i , m ( x ) + K − X i = ( ξ i + − ξ i ) K X ℓ = i + N ℓ, m ( x ) = a + K − X i = ( ξ i + − ξ i ) K X ℓ = i + N ℓ, m ( x ) . Using relations (4) with n = m −
1, (11) with r = δ − i + , m − = Z + ∞−∞ N i + , m − ( u ) du = Z ba N (1) i ( u ) du = ξ i + − ξ i , By the recurrence definition (3): K X ℓ = i + N ℓ, m ( x ) = Z x −∞ δ i + , m − N i + , m − ( u ) du = δ i + , m − Z xa N (1) i ( u ) du . From the partition of unity property of the functions N (1) i ( x ), i = , . . . , K − K X i = ξ i N i , m ( x ) = a + K − X i = Z xa N (1) i ( u ) du = a + Z xa du = x , which implies that ξ , . . . , ξ K are the Greville abscissæ with respect to the basis N . . . , N K and concludes the proof.Note that (12) guarantees that the Greville abscissæ form an increasing sequence.The following result provides an alternative way to calculate the coe ffi cients of reverse knot insertion, which,unlike previous methods [24, Proposition 6], does not involve di ff erentiating the MDB-spline basis.7 roposition 3 (Reverse knot insertion by Greville abscissæ) . Under the same setting and assumptions of Proposition1, denoted by ξ i and ˆ ξ i the Greville abscissae of S and b S , respectively, the coe ffi cients α i in (8) can be computed asfollows: α i = ˆ ξ i − ξ i − ξ i − ξ i − , i = ℓ − d j + , . . . , ℓ − d j + k j . (13) Proof.
Let f be a function in S ⊂ b S with expansions in the MDB-spline bases of S and b S : f ( x ) = K X i = c i N i , m ( x ) = K + X i = ˆ c i ˆ N i , m ( x ) . According to (8) the above coe ffi cients c i and ˆ c i satisfy the relationship:ˆ c i = c i , i ℓ − d j ,α i c i + (1 − α i ) c i − , ℓ − d j + i ℓ − d j + k j , c i − , i > ℓ − d j + k j + , (14)with s ℓ x j < min( s ℓ + , b ). By taking f ( x ) = x , equation (13) is hence obtained from the middle line of (14).By virtue of the above results, the reverse knot insertion step leading from b S to S can be outlined by a triangularscheme of the form (15) D S b S S
GRKI
The “G” arrow and the “RKI” arrow are both needed to evaluate S . More precisely, to evaluate S we need informationabout two spaces. One is the space b S , obtained from S by inserting a knot, the other is the space D S , spanned bythe derivatives of functions in S . Concerning D S we need the MDB-spline basis in order to compute the Grevilleabscissæ of S (Proposition 2), whereas about b S we need the Greville abscissæ in order to compute the RKI coe ffi cients(Proposition 3). Hence we need to put ourselves in a position where the necessary quantities relative to b S and D S are known. To this end, the main idea is to concatenate several blocks of type (15) giving rise to a triangular schemesuch as the one in (17), as will be detailed in the following. This triangular scheme will be constructed in such away that the information at the starting level (the first column of the scheme) is known (or easy to calculate). Theremaining elements in the triangle can be derived by recursive application of the basis block (15), moving from left toright and progressively generating the information relative to each column, up to reaching the vertex of the triangle,which finally contains the information relative to the target space. Remark 3.
As already noted, the reverse knot-insertion formula (9) entails computing higher order derivatives of theB-spline basis functions. On the contrary, the above triangular scheme does not involve the calculation of derivativesof any order. In fact, one might be misled by the fact that the B-spline basis of D S appears in the formula for theGreville abscissa, however, this basis can be evaluated directly (i.e. without involving any di ff erentiation), as we shallsee shortly. Definition 4 ( C r join of two multi-degree spline spaces) . Let S L and S R be MD-spline spaces on adjacent intervals,[ a , b ] and [ b , c ] ( a < b < c ) respectively. We define the C r join of S L and S R to be the MD-spline space on [ a , c ]whose restriction to [ a , b ] and [ b , c ] coincides with S L and S R , respectively, and whose elements are C r continuous at b . It shall be noted that, according to the definition above, the C join of two conventional spline spaces is a C MD-spline space, whereas the C join of two MD-spline spaces is not, in general, a C MD-spline space.One can join with C r continuity two multi-degree spline spaces S L and S R defined on abutting intervals [ a , b ]and [ b , c ] by performing iterated reverse knot insertions at point b , each of which increases the continuity at the8oin. More precisely, denoted by S and S , respectively, the C r and C join of the two spaces S L and S R , we shallstart from S and, by RKI, generate the C join of S L and S R , which we denote by S . We shall then iterate theprocedure generating a sequence of nested spaces S k as in (10), where S k is the C k join of S L and S R , S k + ⊂ S k anddim( S k + ) = dim( S k ) − k = , . . . , r −
1. In this way, the last space S r will be the target space S . These joins will beperformed by concatenating several basic blocks of type (15). To this end, we will need to use the MDB-spline bases N0 L and N0 R of the C MD-spline spaces associated with S L and S R and the representation matrices M L and M R suchthat N L = M L N0 L and N R = M R N0 R , where N L and N R are the MDB-spline bases of S L and S R , respectively. Moreover, on account of the previousdiscussion, we will as well need such information for all the derivative spaces D r − n S L and D r − n S R , for n = , . . . , r ,as detailed in the following.Being K L the dimension of S L , G L = ( ξ L , , . . . , ξ L , K L ) the vector of its Greville abscissæ and N L = ( N L , , . . . , N L , K L )the vector of the MDB-spline basis functions, with a similar notation for S R , the C join of S L and S R , which weindicate by [ S L , S R ], is an MD-spline space of dimension K L + K R − N = [ N L , N R ] ≔ (cid:0) N L , , . . . , N L , K L − , N L , K L + N R , , N R , , . . . , N R , K R (cid:1) , and Greville abscissæ G = [ G L , G R ] ≔ (cid:0) ξ L , , . . . , ξ L , K L ≡ ξ R , , ξ R , , . . . , ξ R , K R (cid:1) . Furthermore, the matrix representation of N relative to the basis N0 = [ N0 L , N0 R ] of the associated C MDB-splinespace is: N = [M L , M R ] N0 , with [M L , M R ] ≔ M L ∗ M R (16)where the above entry ∗ has the same value in M, M L and M R . Note that, if S L is a conventional spline space, thenM L is trivially known, being the identity matrix of size K L , and similarly for M R .We can regard the operation [ , ] as the C join of the representation matrices, of the vectors of MDB-spline basisfunctions or of those of Greville abscissæ. Note that this join operation acts di ff erently depending on the entity towhich it is applied.To describe the process of concatenating several blocks of type (15) we will need to indicate the aforementionedspaces S k , by a double index, that is S k ≕ S r , k , k = , . . . , r . Suppose for example that S L and S R are to be joinedwith C continuity at b to generate space S . Iterated application of (15) will lead to the following triangular schemeof “size” 4: (17)[ D S L , D S R ] = S , ≔ D S , [ D S L , D S R ] = S , ≔ D S , S , ≔ D S , [ D S L , D S R ] = S , ≔ D S , S , ≔ D S , S , ≔ D S , S = [ S L , S R ] = S , S , S , S , = S GG RKI GGRKI GRKI GRKI RKI RKI
More generally, the C r join of S L and S R will involve a similar triangular scheme having size r +
1, where theelement in row n and column k is S n , k , for n = , . . . , r , k = , . . . , n . To implement the scheme all we need to know arethe Greville abscissae and MDB-spline bases of all spaces S n , , n = , . . . , r (that is all spaces in the leftmost column),9ince the information on all other columns and rows of the triangle can be derived by progressive application of thebasic block (15). Recalling that S n , is the C join of D r − n S L and D r − n S R , these Greville abscissæ and MDB-splinebases are easily obtained from the described C join operation [ , ].We shall now use the above triangular scheme to compute the matrix representation (6) where N is the MDB-splinebasis of S and N0 the MDB-spline basis of the associated C MD-spline space. To this aim, observe that all spaces S n , , . . . , S n , n on a row of the triangle are associated with the same C MD-spline space according to Definition 2,which basis will be indicated by N0 n . We may then rewrite the triangular scheme replacing S n , k with its MDB-splinebasis N n , k and Greville abscissæ G n , k and defining M n , k to be the matrix such that N n , k = M n , k N0 n . In particular,matrices M n , are obtained by joining the representation matrices of D r − n S L and D r − n S R as in (16), whereas, foreach k >
0, matrix M n , k is obtained from the preceding one by the relation M n , k = A n , k M n , k − , where A n , k is the thebidiagonal matrix whose nontrivial entries are the RKI coe ffi cients to pass from S n , k − to S n , k .Relying on the matrix representation it is also easy to compute the vectors IN n , k containing integrals of functionsin N n , k that are needed for calculating the Greville abscissæ. In particular, let N nL and N nR be the vectors containing theMDB-spline basis functions in D r − n S L and D r − n S R , respectively, and IN nL and IN nR be the vectors of their integrals.Then the vectors IN n , , corresponding to spaces S n , in the first column of the triangular scheme, contain the integralsof functions in [ N nL , N nR ], n = , . . . , r , and we will indicate this by IN n , = [ IN nL , IN nR ]. For all spaces appearing inthe subsequent columns, instead, the integrals of basis functions are derived from the relation N n , k = M n , k N0 n , whichyields IN n , k = M n , k IN0 n .The described procedure for computing the C r join of spaces S L and S R is outlined in Algorithm 1. The algorithmtakes as input the Greville absissæ, the integrals and the representation matrices of the MDB-spline bases of the twospaces to be joined and of their derivative spaces up to suitable order. On account of their ease of computation, theintegrals of the C MDB-spline bases N0 n (Algorithm 1, line 2) are evaluated at runtime using (5), but they could aswell be provided as input. The algorithm returns as output the representation matrix M r , r , relative to the C MD-splinespace associated with the join space S . Moreover, in anticipation of having to further join the generated MD-splinespace, it computes and returns all the matrices M n , n , n = , . . . , r −
1, related to the derivative spaces S n , n = D r − n S .Note that the overall procedure does never use the MDB-spline bases of the initial spaces S n , , for n = , . . . , r − IN nL , IN nR and IN0 n . Algorithm 1:
Matrix representation of the C r join of two MD-spline spaces S L and S R Data: G nL , IN0 nL , M nL , G nR , IN0 nR , M nR , n = , . . . , r . Result: M n , n , n = , . . . , r . for n ← to r do IN0 n ← [ IN0 nL , IN0 nR ] ; M n , ← [M nL , M nR ] ; G n , ← [ G nL , G nR ] ; end for n ← to r do for k ← to n do IN n − , k − ← M n − , k − · IN0 n − ; Compute G n , k by IN n − , k − using (12) ; Compute the RKI coe ffi cients from G n , k and G n , k − using (13) ; M n , k ← A n , k M n , k − ; end end At this point, by repeatedly joining MD-spline spaces on abutting intervals, we can generate the matrix represen-tation of an MDB-spline basis vector N , spanning an arbitrary space S ≡ S ( P d , X , K ), with respect to the basis vector N0 of the associated C MD-spline space S ( P d , X , K ). Essentially what we need to do is “break” the target spaceinto a sequence of conventional spline spaces and join these spaces in pairs with the required continuities.In particular, with reference to Definition 2, let J be the vector containing the indices, in ascending order, of the10reakpoints separating intervals with di ff erent degrees, including the first and last breakpoint, that isJ = ( j , j , . . . , j p + ) , with 0 = j < j < · · · < j p + = q . Since all breakpoint intervals contained in each [ x j h , x j h + ] have same degree, that is d j h = · · · = d j h + − , we canbreak the target space S into a sequence of conventional spline spaces, each one defined on [ x j h , x j h + ], and thenjoin these spaces two by two. For example, joining the two sections of S relative to [ x j h − , x j h ] and [ x j h , x j h + ] withcontinuity k j h at x j h will produce a space on the whole interval [ x j h − , x j h + ], which is the restriction of the target space S to that interval. This join will generate the representation matrix of the MDB-spline basis of S relative to theMDB-spline basis of S restricted to [ x j h − , x j h + ]. The resulting space can in turn be connected with the neighboringsections of S at x j h − and / or x j h + with continuities k j h − and k j h + , respectively. Not that these joins must be performedin a specific order, namely from higher continuity to lower continuity, in such a way to guarantee that, before eachrepetition of Algorithm 1, all the necessary information (representation matrices and integrals of the MDB-splines)relative to the derivative spaces (up to the required order of di ff erentiation) to the right and left of the join have beengenerated as the output of the previous joins.In this paper, Algorithm 1 mostly serves as a step-up for the derivation of the actual algorithm (see Algorithm2) which will be presented in the next section. Algorithm 2, in fact, is conceptually similar to Algorithm 1 and willbe designed starting from it. In particular, it represents a reformulation which, although less intuitive, allows forimproving the method from a computational point of view.
4. Stable implementation of the RKI Algorithm
In this section we will introduce some observations that will lead us to reformulate Algorithm 1 in an alternativeway, which, although less intuitive, is numerically stable and more e ffi cient from the point of view of the calculationsto be performed.To this end, we start by observing that (13) may raise some concern about the possibile occurrence of cancel-lation errors, due to the di ff erences at the numerators and denominators. The following result shows that the RKIcoe ffi cients can indeed be determined without resorting to the di ff erences of Greville abscissæ, thus it overcomes theaforementioned stability issues. In addition it also improves on the computational cost of the procedure (intended asthe number of operations to be performed) with respect to using (12) and (13). Proposition 4.
The setting and assumptions being the same as in Proposition 3, the RKI coe ffi cients in (13) can becalculated as follows: α i = α (1) i − R ba ˆ N (1) i − ( x ) dx R ba N (1) i − ( x ) dx , i = ℓ − d j + , . . . , ℓ − d j + k j , (18) where α (1) i are the coe ffi cients of reverse knot insertion from D b S to D S .Proof. Under the above assumptions, the MDB-spline basis functions of D b S and D S are such that b N (1) j = N (1) j , for j = , . . . , ℓ (1) − d (1) j − = , . . . , ℓ − d j − d (1) j the degrees in D S and ℓ (1) the index of the largest knot in s (1) smaller or equal to x j ). This observation and relation (7) between the MDB-spline bases of the derivative spacesyield: ˆ ξ i − ξ i − = a + i − X j = Z ba ˆ N (1) j ( x ) dx − a − i − X j = Z ba N (1) j ( x ) dx = i − X j = ℓ − d j Z ba ˆ N (1) j ( x ) dx − i − X j = ℓ − d j Z ba N (1) j ( x ) dx = i − X j = ℓ − d j Z ba ˆ N (1) j ( x ) dx − i − X j = ℓ − d j α (1) j Z ba ˆ N (1) j ( x ) dx + (1 − α (1) j + ) Z ba ˆ N (1) j + ( x ) dx ! . (19)11ence the numerator of (18) comes from the above identity and the fact that α (1) ℓ − d j =
1, whereas the denominatorstraightforwardly follows from (12).
Remark 4.
Using again relation (7) between the MDB-spline bases of the derivative spaces and (18) , we can obtainthe following formula: − α i = (1 − α (1) i ) R ba ˆ N (1) i ( x ) dx R ba N (1) i − ( x ) dx . (20) This result avoids us to actually perform any di ff erences of type − α i or − α (1) i . In particular, when raising thecontinuity from C to C , that is passing from spaces S n in the first column to spaces S n in the second column of thetriangular scheme, the coe ffi cients α (1) i will all be trivial (that is either zero or one) and thus so will be the di ff erences − α (1) i . Hence, at each subsequent iteration, the evaluation of the right-hand side of (20) will just involve thecalculation of the ratio of two integrals and the product by the value − α (1) i inherited from the previous step andtherefore no subtraction will need to be performed. Remark 5.
Relation (18) , which elegantly emerges passing through Greville abscissæ, could alternatively be provenby induction resorting to the integral definition (3) . The latter approach was pursued in a less general context in [11]to determine the coe ffi cients of knot insertion between two MDB-spline bases. Algorithm 2:
Stable matrix representation of the C r join of two MD-spline spaces S L and S R . Data: M nL and M nR , n = , . . . , r ; IN0 nL and IN0 nR , n = , . . . , r − Result: M nn , n = , . . . , r . for n ← to r − do IN nL ← M nL · IN0 nL ; IN nR ← M nR · IN0 nR ; IN n , ← [ IN nL , IN nR ] ; IN0 n ← [ IN0 nL , IN0 nR ] ; M n , ← [M nL , M nR ] ; end M r , ← [M rL , M rR ] ; ibstart ← ℓ , − d , j + for n ← to r do ib ← ibstart ; IN n − , − ib − ← last element of IN n − L ; IN n − , − ib ← first element of IN n − R ; for k ← to n do α n − , k − ib − ← β n − , k − ib + k − ← for i ← ib to ib + k − do α n , ki ← α n − , k − i − IN n − , k − i − / IN n − , k − i − ; β n , ki ← β n − , k − i IN n − , k − i / IN n − , k − i − ; //β n , ki store 1 − α n , ki end M n , k ← A n , k · M n , k − ; if n < r then Compute IN n , k ← M n , k · IN0 n ; end ib ← ib − end ibstart ← ibstart + end C r join of two MDB-spline spaces can be revisited on account of the above discussion,leading to Algorithm 2. In the algorithm, as well as in the examples presented below, we indicate by α n , ki the RKIcoe ffi cients to pass from S n , k − to S n , k and by β n , ki the di ff erences 1 − α n , ki . Furthermore, for each row n of the triangularscheme, we need to identify the index ibstart of the first nontrivial RKI coe ffi cient to be determined. Its initial value(line 9) is derived from Proposition 1 applied to spaces S , and S , , being d , j the degree in the interval to the leftof breakpoint x j (that is, in the algorithm j is the index of the breakpoint where S L and S R are joined and r standsfor k j ) and being ℓ , computed with respect to the left extended partition of S , . Subsequently, the indices of thefirst non-zero RKI coe ffi cients are determined incrementing ibstart while n increases and decrementing it while k increases.The following example not only illustrates the application of Algorithm 2 on a practical case, but also demonstrateshow to generate the matrix representation of an arbitrary MD-spline space following the genaral outline discussed atthe end of section 3, that is by “breaking” the target space into a sequence of conventional spline spaces and joiningthese spaces two by two with the required continuities. Example 1 (Matrix representation via RKI) . Let us consider the target space S ( P d , X , K ) defined on [0 , X = { , , } , degrees d = (2 , , ,
3) and continuities K = (1 , , C MD-spline spacewill be likewise defined on [0 , X and degrees d and will have continuities K = (1 , , S can be seen as the join of three spaces, and more precisely of a degree-2 conventional spline space S A on[0 , S B on [2 ,
3] and a degree-3 polynomial space S C on [3 , C join at point 2 and a C join at point 3. As we will see, these joins should beprocessed starting from the one of higher continuity, since this guarantees that all the information necessary to performa join is either trivially known or has been computed during the previous ones. Therefore we will first calculate the C join of spaces S B and S C at 3, and then calculate the C join of the resulting space with S A at 2.With reference to Algorithm 2, in which S L and S R will be the spaces S B and S C of this example, the first join isdescribed by the triangular scheme (17), which we rewrite below indicating the degrees and continuities in each space S n , k in the form (degree continuity degree), along with the nontrivial RKI coe ffi cients necessary to pass from one spaceto another: S , : (1 S , : (2 S , : (2 S , : (3 S , : (3 S , : (3 S , : (4 S , : (4 S , : (4 S , : (4 α , α , α , i , i = , α , α , i , i = , α , i , i = , , Being the C join of two polynomial spaces, each S n , , n = , . . . ,
3, is a C MD-spline space having dimension2( n + nL and M nR are identity matrices and IN nL ≡ IN0 nL , IN nR ≡ IN0 nR . For n = , . . . ,
2, the integral vectors
IN0 nL , resp. IN0 nR , can be calculated by observing that the basis functions in D r − n S L , resp. D r − n S R , are conventionalB-splines of degree n +
1, resp. n . Hence, according to (5), IN0 nL has n + / ( n +
2) and
IN0 nR has n + / ( n + IN n , , IN0 n and matrices M n , are obtained by the previously discussed C join operation [ , ]; inparticular in this example IN n , ≡ IN0 n = n + , . . . , n + + n + , . . . , n + ! . n = k = α , = α , IN , − IN , = , and β , = β , IN , − IN , = , from which A , =
00 0 , M , = A , , IN , = M , IN0 = , , ! . For n = , k =
1, we obtain: α , = α , IN , − IN , = , and β , = β , IN , − IN , = , from which A , = , M , = A , , and IN , = M , IN0 = , , , , ! . We shall then proceed to n = , k =
2, obtaining: α , = α , IN , IN , = , β , = β , IN , IN , = ,α , = α , IN , IN , = , β , = β , IN , IN , = , which yields A , =
38 914
00 0 0 , M , = A , M , =
58 38
38 2756 914
00 0 0
17 514 and IN , = (1 / , / , / , / . This completes the second row of the triangular scheme. Proceeding in this way for n = k = , ,
3, eventuallyyields the matrix M , = A , M , =
35 720 15
25 2755 2455 411
744 49165 238495 2845
00 0 0 0
115 745 1745 . Recall that Algorithm 2 returns as output the representation matrices for all derivative spaces up to di ff erentiationorder three of the C join of S B and S C . We shall use this information to apply the algorithm again, this time for14oining with C continuity the conventional spline space S A on [0 ,
2] and the MD-spline space on [2 ,
4] obtained asoutput of the previous join. As for this second round, the triangular scheme will be: S , : (0 − S , : (1 S , : (1 S , : (2 S , : (2 S , : (2 α , α , α , i , i = , Denoted as usual by S L and S R the two spaces to be joined, the corresponding representation matrix M nL will be theidentity of size n +
2, whereas M nR = M n + , n + , being M n + , n + , n = , ,
2, the output of the previous C join. Theintegrals of the C MDB-spline functions required by Algorithm 2 can be evaluated by (5) for the left-hand side spaces S nL , n = ,
1, which gives:
IN0 L = (1 ,
1) and
IN0 L = , , ! . For spaces S nR , n = ,
1, instead, the integrals may be stored as output of the first join or e ffi ciently calculated atruntime by (5), obtaining: IN0 R = , , , ! and IN0 R = , , , , , ! , hence, from lines 2 and 3 of the algorithm, we will obtain IN nL = IN0 nL , n = ,
1, and IN R = , , ! and IN R = , , , ! . Finally, vectors IN n , and IN0 n , n = ,
1, (lines 4 and 5) and matrices M n , , n = , ,
2, (line 6) are the C join of theabove quantities. The output of this second and last join is a matrix M , such that N = M , N0 . The above is the matrix representation of the MDB-spline basis of the target space S relative to the basis of theassociated C MD-spline space N0 = N0 .As previously mentioned, note that processing the joins from higher to lower continuity makes so that, each time,all the information to address the next join is available or has been computed during the previous steps.With the previous example in mind, we can further discuss some details of our implementation. In order tosave on memory allocation, only one matrix M n , k should be stored for each row of the triangular scheme, that iseach n = , . . . , r , since such matrices can be overwritten when moving from one column to the other. In addition,it is unnecessary to create matrices A n , k , which we merely introduced for ease of presentation, as the coe ffi cients α n , ki and β n , ki can be stored in temporary one-dimensional arrays, to be destroyed after been used for the coe ffi cientcomputations at lines 15 and 16 and the RKI steps at lines 18 and 19. The integrals IN n , k can as well be stored intemporary one-dimensional arrays. Moreover, only one array can be used to store all vectors IN0 n , n = , . . . , r −
5. Stability analysis
Unlike how it usually happens, namely that we propose an algorithm and then we analyze its stability, we designedan algorithm that would possess all the characteristics to be numerically stable. This feature becomes clear if we15reak Algorithm 2 into a sequence of basic steps, each involving numerically stable operations only. The results ofthis analysis will be confirmed and highlighted by the numerical experimentation presented in subsection 5.1.Our discussion may benefit from some preliminary considerations. First, it is easy to count how many RKIcoe ffi cients will be calculated over the course of the algorithm. In particular, the “for” loops at lines 14 and 17 showthat we will have to calculate one coe ffi cient α n , i for n = , . . . r , two coe ffi cients α n , i for n = , . . . , r and so on upto r coe ffi cients α n , ri for n = r . Also note that, for k =
1, at lines 18 and 19 reference is made to elements of vectors IN n − , − , never formally initialized, but whose values are trivially known from IN nL and IN nR (lines 12 and 13). Again,for each n and k , the first α n − , k − i and the last β n − , k − i considered are trivially equal to one (lines 15 and 16). Finally, asrecalled earlier, storing the quantities β n − , k − i allows us to avoid the evaluation of the quantities 1 − α n , ki and thereforethe whole algorithm does not contain any floating point subtraction.Bearing in mind these observations, we can break the algorithm in the following basic steps.A) Calculation of the input vectors IN0 nR and IN0 nL . Since functions in N0 nR and N0 nL are C MDB-splines, theevaluation of their integrals involves computing and adding the integrals of conventional B-splines according to(5), all of which are positive quantities. Likewise, the C join of the integral vectors at lines 4 and 5 involvessummations between positive quantities.B) Products between matrices M n , k (as well as M nL and M nR ), all of which entries belong to [0 , IN0 n (as well as IN0 nL and IN0 nR ) (lines 2,3 and 23). Due to the fact that only some elements of the vectors atthe right-hand side of these assignments are used, these products are reduced to dot products between single rowsof matrices M n , k and vectors IN0 n . Note that each iteration involves as many such dot products as the integralsat lines 18 and 19, that is 3 dot products at most (since some of those integrals are used twice, so they could bestored and reused).C) Evaluation of the right-hand sides of the assignments at lines 18 and 19. This amounts to calculating first theproduct, which produces a value in [0 , ,
1] as can be seen from thefact that N n − , k − i − = α n − , k − i − N n − , k − i − + (1 − α n − , k − i ) N n − , k − i .D) Product at the right-hand side of the assignment at line 21. Rather than a matrix product, it is convenient toperform this calculation as a repeated combination of two rows of M n , k − (all of which entries are in [0 , α n , kj m n , k − j − + β n , kj + m n , k − j , where α n , kj and β n , kj + = (1 − α n , kj + ) are entries on the bidiagonal of A n , k and m n , k − j isthe j th row of M n , k − .The above analysis emphasizes that the proposed algorithm consists of summations, ratios and products betweenpositive quantities (most of which belonging to [0 , C r join of two MD-spline spaces, that is: • r operations of type A); • ( r − r operations of type B), or 3(2( r − + r − + . . . ( r − + r ) dot products; • r + r − + r − + . . . ( r −
1) 2 + r operations of type C); • r ( r + operations of type D) or 2 r + r − + r − + . . . + r + r + n , k − , thatis as many as the overall number of nontrivial RKI coe ffi cients α n , ki plus one;and thus to estimate the computational complexity of the algorithm, which amounts to O ( r ) operations. Besides supporting the conclusions of the above stability analysis, the following numerical experiments providea comparison between the new proposal and previous ones. For the sake of brevity, we will refer to the presentmethod and to those in [24] and [22] as RKI / Greville, RKI / Derivative and H-Operator, respectively. Recall thatboth the RKI / Derivative and H-Operator algorithms make use of derivatives (of order up to the target continuity) ofMDB-splines and thus su ff er in a similar way from the fact that those quantities may be very large numbers.16ur analysis is based on calculating and comparing the algorithmic errors on the evaluation of MDB-spline basisfunctions and / or on the representation matrix. To this end, the “exact” values are obtained by symbolic computation,using MATLAB’s Symbolic Math Toolbox, whereas the numerical results rely on MATLAB’s standard precision(rounding unit U ≈ − ). In all the examples, the symbolic implementation of the RKI / Greville algorithm wasable to produce an output within reasonable time, due to the fact that the method performs operations between smallquantities all of which can be stored in rational form. As would be expected, the response times of the symbolicprocedure become impractical for more complex tests.This section contains three experiments. The first (Example 2) is aimed at evaluating how our analysis approach,based the algorithmic error, relates to the a posteriori error bound in Cox’s seminal paper on the evaluation of B-splines [3]. Like the referenced paper, this example is concerned with conventional B-splines and as a consequencethe representation matrix is the identity matrix. In the successive two experiments (Examples 3 and 4) we compare theRKI / Greville Algorithm with previous proposals on a variety of test spaces featured by both uniform and nonuniformdistributions of breakpoints as well as largely inhomogeneous degrees. The parameters of the di ff erent test spaces thatwill be considered are summarized in Table 1. [ a , b ] X d K Test 1 [ − , {− , , } (5 , , ,
5) (3 , , − , {− , , } (3 , , ,
3) (3 , , , { j } , j = , . . . , , , , , , , , , ,
9) (8 , , , , , , , , − , {− − j } , j = , . . . , , , , , , , , , ,
9) (8 , , , , , , , , , { j } , j = , . . . , d i = i = , . . . , k i = i = , . . . , d i = i = , . . . , , , . . . , k i = i = , . . . , , , . . . , d i = i = , . . . , , , . . . , k i = i = , . . . , , , . . . , − , {− , , } (21 , , ,
21) (15 , , Table 1: Test spaces for Examples 3 and 4.
Example 2 (A comparison with conventional B-splines) . This experiment replicates [3, Example 2], which is themost challenging test in the referenced paper. The setting is a conventional spline space of degree 21, defined in theinterval [0 , x j placed at the integers and C continuity at each breakpoint. Notethat choosing both the breakpoints and the evaluation points to be exactly represented in the floating point standardallows for avoiding roundo ff errors in the initial data. Table 2 shows the algorithmic error on the evaluation of the“central” B-spline N , at the breakpoints x j . For the same experiment, [3, Table 2] reports the values of N , alongwith the a posteriori error bounds established in that paper. In particular, the values of N , found by Cox referto non-normalized basis functions and are the same as in the second column of Table 2, whereas the values in thethird column of Table 2 are obtained with the recurrence relation for normalized C MDB-splines in [24], which is asimple generalization of the more established scheme in [4]. The values in the two columns, however, only di ff er bya normalization constant equal to the width of the support.A running error analysis was also integrated in our implementation and returned a posteriori error bounds inaccordance with those reported by Cox (considering that we work in double precision with 16 digits, while Cox with11 digits). It shall be noted, in particular, how the results in the column of absolute algorithmic errors are consistentwith the corresponding error bounds and the corresponding relative errors that will be used to assess the numericalstability of our proposal.In the conclusions of [3], on the basis of the a posteriori error bound, it is expected that the maximum relativeerror attained with a t -digits mantissa cannot exceed (70)2 − t for degree 10 or less, whereas it cannot exceed (700)2 − t for degree 100 or less. It is also observed that the bound on the relative error grows linearly with the degree of aspline. Our experimentation shows that the actual error is even lower. In fact, for degree 100 or less the relativealgorithmic error for most experiments is about 10 − , with only a few values of the order of 10 − , whereas the boundestimated by Cox would be of the order of 10 − . We believe that this may be attributable to cancellation of roundingerror, which may cause the final computed answer to be much more accurate than the intermediate quantities. Thisphenomenon has been described, e.g., in [26, p.19]. 17 j N , ( x j ) Non-Normalized N , ( x j ) Normalized Error Bound Absolute Alg. Error Relative Alg. Error1 8.896791392450574e-22 1.957294106339126e-20 3.2378e-34 1.3644e-36 6.9706e-172 1.865772813284987e-15 4.104700189226971e-14 6.7901e-28 4.0230e-30 9.8009e-173 9.265310806863227e-12 2.038368377509910e-10 3.3719e-24 3.3787e-26 1.6575e-164 3.708541354285271e-09 8.158790979427597e-08 1.3497e-21 7.0656e-24 8.6601e-175 3.402962627063746e-07 7.486517779540241e-06 1.2384e-19 9.0997e-23 1.2155e-176 1.107329203006056e-05 2.436124246613324e-04 4.0299e-18 2.4981e-20 1.0254e-167 1.595958078468785e-04 3.511107772631326e-03 5.8082e-17 9.0643e-19 2.5816e-168 1.156908330166488e-03 2.545198326366273e-02 4.2103e-16 3.6835e-18 1.4472e-169 4.554285942496692e-03 1.001942907349272e-01 1.6574e-15 7.1213e-18 7.1075e-1710 1.019454972176512e-02 2.242800938788327e-01 3.7101e-15 2.1568e-17 9.6165e-1711 1.330103123779249e-02 2.926226872314347e-01 4.8407e-15 8.2012e-17 2.8026e-1612 1.019454972176512e-02 2.242800938788327e-01 3.7101e-15 2.1568e-17 9.6165e-1713 4.554285942496692e-03 1.001942907349272e-01 1.6574e-15 7.1213e-18 7.1075e-1714 1.156908330166488e-03 2.545198326366273e-02 4.2103e-16 3.6835e-18 1.4472e-1615 1.595958078468785e-04 3.511107772631326e-03 5.8082e-17 9.0643e-19 2.5816e-1616 1.107329203006056e-05 2.436124246613324e-04 4.0299e-18 2.4981e-20 1.0254e-1617 3.402962627063746e-07 7.486517779540241e-06 1.2384e-19 9.0997e-23 1.2155e-1718 3.708541354285271e-09 8.158790979427597e-08 1.3497e-21 7.0656e-24 8.6601e-1719 9.265310806863227e-12 2.038368377509910e-10 3.3719e-24 3.3787e-26 1.6575e-1620 1.865772813284987e-15 4.104700189226971e-14 6.7901e-28 4.0230e-30 9.8009e-1721 8.896791392450574e-22 1.957294106339126e-20 3.2378e-34 1.3644e-36 6.9706e-17 Table 2: Numerical experiments reported in Example 2
We conclude by mentioning that a similar study of algorithmic errors was carried out on the evaluation of deriva-tives. Also in this case for splines of degree less than or equal to 50 and order of di ff erentiation up to ten we neverencountered algorithmic errors exceeding ≈ − . Example 3 (Algorithmic error on the evaluation of MDB-splines) . This experiment illustrates how erroneous theresults of the RKI / Derivative method can be for degrees as low as three and five if the knot spacing is highly nonuni-form. Such a case is important in practice since it is often of interest to investigate the case of near-coincident knots.From Table 3 one can observe that at x = − x = / Greville agree for sym-metry, while this is not the case for the corresponding results obtained by RKI / Derivative. Moreover, the values of thealgorithmic errors show that the accuracy of the RKI / Derivative method is limited to the first 6 / x =
0. Similar results are reported in Table 4, from whichone can again see that the RKI / Derivative method returns strongly asymmetric results despite the expected symmetryof the evaluated MDB-spline. In both experiments the results obtained by RKI / Greville agree for symmetry and areextremely accurate, which is consistent with the conclusions of the theoretical analysis.The experiment reported in Table 5 concerns a space with a less challenging uneven distribution of breakpoints,but higher degrees. In this case, the RKI / Derivative method appears adequate up to 9 figures only. Analogous resultswere also obtained for the spaces “Test 4” , “5” and “6”. Overall, the large errors for the RKI / Derivative algorithmshow that the method is potentially unstable. Conversely the small algorithmic errors of the RKI / Greville methodconfirm its stability. The same conclusions are supported by the results illustrated in Example 4, concerned with thealgorithmic errors with respect to the entries of the representation matrices.
Example 4 (Algorithmic error on the representation matrix) . In this second type of test, we consider the algorithmicerror on the representation matrix, calculated as k M
16 digits − M exact k , where M
16 digits is the numerically calculated matrix, whereas M exact is the one obtained by symbolic computation.This is both an absolute and a relative error on account of the fact that k M exact k = KI / Greville RKI / Derivativex N ( x ) Normalized Relative Alg. Error N ( x ) Normalized Relative Alg. Error-9.999000e +
03 4.500275008083014e-09 1.8381e-16 4.500275772672185e-09 1.6990e-070.000000e +
00 5.000083333610773e-01 0.0000e +
00 5.000084045999867e-01 1.4248e-079.999000e +
03 4.500275008083015e-09 0.0000e +
00 4.500275649258610e-09 1.4247e-07
Table 3: Numerical results discussed in Example 3 for space “Test 1” in Table 1.
RKI / Greville RKI / Derivativex N ( x ) Normalized Relative Alg. Error N ( x ) Normalized Relative Alg. Error-9.999000e +
03 2.499250262410031e-12 0.0000e +
00 2.499214206146373e-12 1.4427e-050.000000e +
00 3.750749868799358e-01 0.0000e +
00 3.750749863390447e-01 1.4421e-099.999000e +
03 2.499250262410030e-12 1.6161e-16 2.499250262410031e-12 1.6161e-16
Table 4: Numerical results discussed in Example 3 for space “Test 2” in Table 1.
RKI / Greville RKI / Derivativex N ( x ) Normalized Relative Alg. Error N ( x ) Normalized Relative Alg. Error2.000000e +
00 2.912087112938504e-13 3.4674e-16 2.912087106308203e-13 2.2768e-094.000000e +
00 1.275774160308294e-09 1.6209e-16 1.275774157784237e-09 1.9785e-098.000000e +
00 4.806036147184862e-07 2.2030e-16 4.806036141267946e-07 1.2311e-091.600000e +
01 5.258129295850228e-05 3.8662e-16 5.258129293072319e-05 5.2831e-103.200000e +
01 2.147713272383253e-03 8.0771e-16 2.147713271996800e-03 1.7994e-106.400000e +
01 3.541058939374863e-02 5.8787e-16 3.541058939374988e-02 3.4684e-141.280000e +
02 2.206016671195212e-01 3.7745e-16 2.206016671340502e-01 6.5860e-112.560000e +
02 3.592347216925473e-01 0.0000e +
00 3.592347217235125e-01 8.6198e-115.120000e +
02 4.466585515804859e-02 1.5535e-16 4.466585516215183e-02 9.1866e-11
Table 5: Numerical results discussed in Example 3 for space “Test 3” in Table 1.
The algorithmic error obtained with RKI / Greville is compared with those relative to both RKI / Derivative andH-operator (the code for the latter is taken from [22]. ).Table 6 contains the algorithmic errors obtained for all the test spaces in Table 1. In particular, space “Test 5”is the multi-degree counterpart of the aforementioned experiment [3, Example 2]. “Test 6”, instead, is aimed atcomparing the considered algorithms in case of a very nonuniform partition and high degrees. Finally, Table 7 showsthe algorithmic errors obtained in a test case presented in our previous paper [24]. All the results confirm the adequacyof the new proposal, by contrast with previous methods, which, in some cases, su ff er from serious loss of accuracy.
6. Matrix representation in terms of the conventional B-spline basis of maximum degree
For a given target space S ( P d , X , K ), another way to compute a matrix representation (6) is to choose as initialspace a conventional spline space S ≡ S ( P d , X , K ), with d j = m for all j , being m ≔ max j { d j } . In this setting, itstill holds that S ⊂ S , but this time matrix M needs to be computed performing successive steps of reverse degreeelevation (RDE). As the name suggests, reverse degree elevation is the reverse operation of degree elevation and wecan understand it by referring to Remark 1, where instead of decreasing / increasing the continuity at a breakpoint, oneincreases / decreases the degree on a breakpoint interval. Therefore, each round of reverse degree elevation diminishesby one the degree in an interval, until each interval [ x j , x j + ] reaches the target degree d j . Overall, the number of steps g required to pass from S to S amounts to the total number of RDE steps to be performed, that is g ≔ P qj = ( m − d j ) . The H-operator algorithm, as implemented in [22], returns a representation matrix with respect to a sequence conventional B-spline basesconnected with C − continuity, therefore, with respect to ours, it has replicated columns that have been removed for a fair comparison. xamples Algorithmic Errors Examples Algorithmic ErrorsRKI / Greville RKI / Derivative H-Operator RKI / Greville RKI / Derivative H-OperatorTest 1 1.0 × − × − × − Test 4 6.0 × − × − × − Test 2 6.7 × − × − × − Test 5 1.0 × − × − × − Test 3 3.7 × − × − × − Test 6 1.7 × − × + × + Table 6: Algorithmic errors on the representation matrix for the test spaces in Table 1. k K Algorithmic Errors k K Algorithmic ErrorsRKI / Greville RKI / Derivative H-Operator RKI / Greville RKI / Derivative H-Operator5 17 2.5 × − × − × −
13 25 2.7 × − × − × − × − × − × −
15 27 4.4 × − × − × − × − × − × −
17 29 3.1 × − × − × −
11 23 2.5 × − × − × −
19 31 4.5 × − × − × − Table 7: Target space S ( P d , X , K ) with [ a , b ] = [0 , X = (1) and d = (19 , S is K = The process must be accomplished in such a way to generate a sequence of MD-spline spaces S n ≡ S ( P d n , X , K n ), n = , . . . , g , such that S ≔ S g ⊂ S g − ⊂ · · · ⊂ S ⊂ S , (21)where each space S n is defined on [ a , b ], has same breakpoints X and continuities K as the target space S and hasdimension K n ≔ K + ( g − n ), being K the dimension of S . In general, there may be more than one sequence (21)leading from S to S g and therefore, while S and S g are fixed, the intermediate spaces S , . . . , S g − will depend onthe specific ordering of RDE steps performed.[24, Proposition 7] provides a result akin to Proposition 3, where space b S is obtained from S through (local) degreeelevation. In this case, the respective MDB-spline bases satisfy a relationship analogous to (8), with coe ffi cients α i given by (9), the only di ff erence being that the nontrivial coe ffi cients α i ∈ ]0 ,
1[ correspond to i = ℓ − d j + , . . . , ℓ . Thesecoe ffi cients can still be determined through (13), where ˆ ξ j and ξ j are the Greville abscissæ of b S and S , respectively.On account of Proposition 2, the computation of the Greville abscissæ of b S and S entails integrating the MDB-spline bases of the respective derivative spaces. Hence, an RDE step can be described by the following triangularblock, akin to (15): D S b S S
GRDE
Repeated applications of the above basic block, give rise to the following rhomboid scheme, which is the RDEcounterpart of (17), and in which S k , n , for k = , . . . , r , n = , . . . , g , indicate the derivative spaces D r − k S n with r ≔ max { , max { k i ∈ K}} : S , S , · · · S , g − S , g · · · · · · · · · · · · · · ·S r − , S r − , · · · S r − , g − S r − , g S r , S r , · · · S r , g − S r , gGRDE GRDE GRDE GRDEG GRDE RDEG GGRDE GRDE GRDE GRDE RDE RDE S , n are C MD-spline spaces and may feature breakpoints with negative continuities, as well asintervals with negative degrees. These correspond to the degenerate spaces involved in the generation of the MDB-spline basis of S r , n by the integral recurrence relation (3). Spaces S k , are instead conventional spline spaces ofdegree m − ( r − k ). For all spaces S , n and S k , the MDB-spline basis functions, as well as their integrals, can bestraightforwardly computed by standard approaches, as discussed in section 2.2.Using the rhomboid scheme and the corresponding matrix representations leads to Algorithm 3, where K ( k , n )indicates the dimension of space S k , n , that is K ( k , n ) = K − ( r − k ) + ( g − n ), being K the dimension of the target space S ≡ S r , g . The algorithm requires as input the vectors IN k , , k = , . . . , r , containing the integrals of the conventionalB-spline bases of spaces S k , and the vectors IN , n , n = , . . . , g , of the integrals of the C MDB-splines of the spaces S , n . It returns as output the matrices M k , g such that N k , g = M k , g N k , , k = , . . . , r , where N k , is a conventionalB-spline basis of degree m − r + k . In particular M r , g is the matrix representation of the MDB-spline basis N r , g of thetarget space S g ≡ S r , g with respect to the B-spline basis N r , of the conventional spline space S ≡ S r , of degree m ≔ max i { d i } . We remark that, while the RKI Algorithm joins two spaces at a time, the RDE works globally, i.e. bycarrying out a sequence of reverse degree elevations on all the intervals involved. Algorithm 3:
Matrix representation relative to the conventional B-spline basis of degree m ≔ max i { d i } , with r ≔ max { , max { k i ∈ K}} and g ≔ P qj = ( m − d j ). Data: IN k , , k = , . . . , r ; IN , n , n = , . . . , g . Result: M k , g , k = , . . . , r . for k ← to r do M k , ← I K ( k , ; n ← for j ← to q do for h ← m − to d j do n ← n + ie ← d k , n + + P jh = d k , nh − k k , nh ; ib ← ie − d k , nj + α k − , nib − ← β k − , nie ← for i ← ib to ie do if k == then α , ni ← (cid:16)P i − h = ib − IN , n − h − P i − h = ib − IN , nh (cid:17) / IN , ni − ; β , ni ← (cid:16)P i − h = ib − IN , nh − P i − h = ib − IN , n − h (cid:17) / IN , ni − ; else α k , ni ← α k − , ni − IN k − , n − i − / IN k − , ni − ; β k , ni ← β k − , ni IN k − , n − i / IN k − , ni − ; // β k , ni store 1 − α k , ni end end M k , n ← A k , n M k , n − ; if k < r then IN k , n ← M k , n · IN k , ; end end end end Example 5 (Matrix representation via reverse degree elevation) . In the interval [0 , S ( P d , X , K ) with X = { , } , d = (4 , ,
3) and K = (2 , S , S and S defined on the same interval and having same breakpoint sequence and continuities as S ≡ S and such that S ⊂ S ⊂ S ⊂ S . In particular we take S to be the MD-spline space having d = (4 , , S having d = (4 , , S having d = (4 , ,
4) and S having d = (4 , , S is a conventional spline space and hence its B-spline21asis and corresponding integrals can e ffi ciently be computed by known methods. In this way one can pass from theMDB-spline basis of S to that of S , from that of S to that of S and finally from the MDB-spline basis of S tothat of S performing three successive rounds of RDE.The rhomboid scheme of spaces in this example is as follows, where S n ≡ S , n , n = , . . . , S , : (2 − S , : (2 − S , : (2 − S , : (2 − S , : (3 , S , : (3 S , : (3 S , : (3 S , : (4 S , : (4 S , : (4 S , : (4 α , α , α , i i = , α , α , i i = , α , i , i = , , α , i , i = , α , i , i = , , The nontrivial coe ffi cients α k , ni necessary for each RDE step can be computed by integrating the MDB-splinefunctions of derivative spaces as in (18). We remark that a similar result was proven in [28] for a less general subclassof MD-splines. For instance, knowing the MDB-splines of space D S ≡ D S , = S , , which is defined on the sameinterval and breakpoint sequence as S , but has degrees d = (3 , ,
3) and continuities K = (1 , ffi cients α , i , i = , ,
6, to pass from S ≡ S , to S ≡ S , . With reference to the first line of the rhomboidscheme, observe how going from S , to S , it is necessary to calculate only one coe ffi cient, as one passes fromdegree 2 to 1. For the same reason it is necessary to calculate only one coe ffi cient α , for the RDE step from S , to S , . The RDE step from S , to S , , instead, does not involve non-trivial RDE coe ffi cients, as they are all equal to 0or 1. Remark 6.
The procedure can be modified in such a way to avoid any subtraction operation and therefore improveits numerical stability. In fact, the coe ffi cients in the first row of the rhomboid scheme (in the example α , and α , )are determined by (13) , whose numerator can be computed by the middle line of (19) (Algorithm 3, lines 13 and 14).However, we can as well further di ff erentiate these spaces, in such a way that the first two rows of the scheme become: (1 − −
1) (1 − −
1) (1 − − −
1) (1 − − − S , : (2 − S , : (2 − S , : (2 − S , : (2 − α , α , In this extended version of the scheme also the coe ffi cients α , and α , can be determined through (18) , avoidingthe aforementioned di ff erences. This variant of Algorithm 3 can be obtained by defining r asr ≔ max { d i } − and summarizing lines from 12 to 18 of the algorithm by lines 16 and 17 only. The RDE-based algorithm is numerically stable for the same considerations made in the RKI case and all thenumerical tests carried out have verified its excellent accuracy in the calculation of both the MDB-spline functionsand the representation matrix.
Remark 7 (Mixed RDE-RKI Algorithm) . It is also possible to design an algorithm that simultaneously performsRDE and RKI steps, like the one proposed in [24]. In this case we shall choose the initial space S , containing S , insuch a way that S can be reached through a sequence of successive steps of RDE and RKI type. We shall hence breakthe target space S into sections, each of which will be generated from the corresponding section of S via RDE usingAlgorithm 3 or through RKI joins of the corresponding sections in S via Algorithm 2. At this point it is necessary toproceed by first addressing all the sections requiring RDE, obtaining the representation matrices of the correspondingMDB-spline bases, and then joining by RKI the resulting MD-spline spaces, starting from the one with the highestcontinuity up to the one with the least continuity. . Conclusions We have presented an algorithm for the e ffi cient evaluation of multi-degree B-splines, which, unlike previousapproaches, is numerically stable. This has been emphasized via theoretical analysis of the involved operations, as wellas by numerical experiments and comparisons with previous methods. From the point of view of numerical stability,the proposed method is at present the most e ff ective tool for evaluating multi-degree splines. Furthermore, similarideas could be employed in the more general context of piecewise Chebyshevian splines of variable dimensions, whichhave been the subject of recent studies [19, 23]. Acknowledgements
The authors gratefully acknowledge support from INdAM-GNCS Gruppo Nazionale per il Calcolo Scientifico.
ReferencesReferences [1] de Boor, C.. A Practical Guide to Splines. New York: Springer-Verlag; 1978. doi:10.2307 / / CBO9780511618994.[3] Cox, M.. The numerical evaluation of B-splines. J Inst Maths Applics 1972;10:134–149. doi:10.1093 / imamat / / / j.cma.2016.11.009.[6] N¨urnberger, G., Schumaker, L.L., Sommer, M., Strauss, H.. Generalized Chebyshevian splines. SIAM J Math Anal 1984;15(4):790–804.doi:10.1137 / / S0167-8396(03)00096-7.[8] Liu, L., Casquero, H., Gomez, H., Zhang, Y.J.. Hybrid-degree weighted t-splines and their application in isogeometric analysis. Computers& Fluids 2016;141:42 – 53. doi:https: // doi.org / / j.compfluid.2016.03.020. Advances in Fluid-Structure Interaction.[9] Thomas, D.C., Engvall, L., Schmidt, S.K., Tewa, K., Scott, M.A.. U-splines: Splines over unstructured meshes; 2018. Coreform report.[10] Shen, W., Wang, G.. A basis of multi-degree splines. Comput Aided Geom Design 2010;27(1):23–35. doi:10.1016 / j.cagd.2009.08.005.[11] Shen, W., Wang, G.. Changeable degree spline basis functions. J Comput Appl Math 2010;234(8):2516–2529.doi:10.1016 / j.cam.2010.03.015.[12] Shen, W., Wang, G., Yin, P.. Explicit representations of changeable degree spline basis functions. J Comput Appl Math 2013;238(1):39–50.doi:10.1016 / j.cam.2012.08.017.[13] Beccari, C.V., Casciola, G.. A Cox-de Boor-type recurrence relation for C multi-degree splines. Comput Aided Geom Design2019;75:101784–101784. doi:https: // doi.org / / j.cagd.2019.101784.[14] Li, X., Huang, Z.J., Liu, Z.. A geometric approach for multi-degree spline. Journal of Computer Science and Technology 2012;27(4):841–850. doi:10.1007 / s11390-012-1268-2.[15] Buchwald, B., M¨uhlbach, G.. Construction of B-splines for generalized spline spaces generated from local ECT-systems. J Comput ApplMath 2003;159(2):249–267. doi:10.1016 / S0377-0427(03)00533-8.[16] Beccari, C., Casciola, G., Morigi, S.. On multi-degree splines. Comput Aided Geom Design 2017;58:8–23. doi:10.1016 / j.cagd.2017.10.003.[17] Antonelli, M., Beccari, C.V., Casciola, G.. A general framework for the construction of piecewise-polynomial local interpolants of minimumdegree. Adv Comput Math 2014;40(4):945–976. doi:10.1007 / s10444-013-9335-y.[18] Beccari, C.V., Casciola, G., Romani, L.. Construction and characterization of non-uniform local interpolating polynomial splines. J ComputAppl Math 2013;240:5–19. doi:10.1016 / j.cam.2012.06.025.[19] Beccari, C.V., Casciola, G., Romani, L.. Computation and modeling in piecewise Chebyshevian spline spaces; 2017. ArXiv:1611.02068.[20] Mazure, M.L.. How to build all Chebyshevian spline spaces good for geometric design? Numer Math 2011;119(3):517–556.doi:10.1007 / s00211-011-0390-3.[21] Toshniwal, D., Speleers, H., Hiemstra, R.R., Manni, C., Hughes, T.J.. Multi-degree B-splines: Algorithmic computation and properties.Comput Aided Geom Design 2020;76:101792–101792. doi:https: // doi.org / / j.cagd.2019.101792.[22] Speleers, H.. Algorithm 999: Computation of multi-degree B-splines. ACM Transactions on Mathematical Software 2019;45(4):1–15.doi:10.1145 / ffi an extension of multi-degree B-splines: Algorithmiccomputation and properties. SIAM Journal on Numerical Analysis 2020;2(58):1138–1163. doi:https: // doi.org / / // doi.org / / j.cam.2020.113007.[25] Butterfield, K.R.. The computation of all derivatives of a B-spline basis. J Inst Maths Applics 1976;17:15–25. doi:10.1093 / imamat / /
27] Carnicer, J.M., Mainar, E., Pe˜na, J.M.. Greville abscissae for totally positive bases. Comput Aided Geom Design 2016;48:60–74.doi:10.1016 / j.cagd.2016.09.001.[28] Shen, W., Yin, P., Tan, C.. Degree elevation of changeable degree spline. Journal of Computational and Applied Mathematics 2016;300:56– 67. doi:http: // dx.doi.org / / j.cam.2015.11.030.j.cam.2015.11.030.