Tension spline of order $k$ is a function that, for a given partition
mbox{$x_0<x_1<ldots<x_n$}, on each interval mbox{$[x_i,x_{i+1}]$} satisfies differential equation mbox{$(D^k-rho_{i}^2 D^{k-2})u=0$}, where $rho_i$\u27s are prescribed nonnegative real numbers.
Most articles deal with tension splines of order four, applied to the 
problem of convex (or monotone) interpolation or to the two-point boundary value problem for ordinary differential equations. Higher order tension splines are described in several papers, but no application is given. Possible reason for this is a lack of an appropriate algorithm for their evaluation.
Here we present an explicit algorithm for evaluation of tension splines of arbitrary order. We especially consider stable and accurate computation of hyperbolic-like functions used in our algorithm

I. Beroš

M. Marušić

Hrčak - Portal of scientific journals of Croatia

Mathematical Communications 4(1999), 73-81 73Evaluation of tension splines∗Ivo Berosˇ† and Miljenko Marusˇic´‡Abstract. Tension spline of order k is a function that, for a givenpartition x0 < x1 < . . . < xn, on each interval [xi, xi+1] satisfies differ-ential equation (Dk − ρ2iDk−2)u = 0, where ρi’s are prescribed nonneg-ative real numbers.Most articles deal with tension splines of order four, applied to theproblem of convex (or monotone) interpolation or to the two-point bound-ary value problem for ordinary differential equations. Higher order ten-sion splines are described in several papers, but no application is given.Possible reason for this is a lack of an appropriate algorithm for theirevaluation. Here we present an explicit algorithm for evaluation of ten-sion splines of arbitrary order. We especially consider stable and accu-rate computation of hyperbolic-like functions used in our algorithm.Key words: tension spline, B-spline, evaluationAMS subject classifications: 65D07, 65-041. PreliminariesTension spline of order k is a function that, for given partition x0 < x1 < . . . < xn,on each interval [xi, xi+1] satisfies differential equation(Dk − ρ2iDk−2)u = 0,where ρi’s are prescribed nonnegative real numbers. In other words, tension splinesof order k are functions whose restrictions on non-empty interval (xi, xi+1) lie inspan{1, x, . . . , xk−1}, for ρi = 0, or in span{eρix, e−ρix, 1, x, . . . , xk−3}, otherwise.Almost all papers dealing with tension splines consider tension spline of orderfour. Applications of such splines are mostly shape preserving approximation [9, 10]and a numerical solution of singularly perturbed two-point boundary value problemfor ODE [1, 5]. Higher order tension splines are described in several papers [2, 4],but no application is given. Their application may be foreseen in a problem of∗This paper was presented at the Special Section: Appl. Math. Comput. at the 7th InternationalConference on OR, Rovinj, 1998.†Department of Mathematics, University of Zagreb, Bijenicˇka 30, HR-10 000 Zagreb, Croatia,e-mail: ivo.beros@math.hr‡Department of Mathematics, University of Zagreb, Bijenicˇka 30, HR-10 000 Zagreb, Croatia,e-mail: miljenko.marusic@math.hr74 I. Berosˇ and M.Marusˇic´approximation by functions with positive higher derivative, as well as in singularlyperturbed boundary value problems of higher degree. A possible reason for theabsence of their application is a lack of an appropriate algorithm for their evalua-tion. Here we present an explicit algorithm for the evaluation of tension splines ofarbitrary order.2. B-splinesAs for the case of polynomial splines, we use a B-spline representation for tensionsplines (cf. [2, 3, 4, 6, 7]).Suppose m,n are positive integers and m ≥ 2. Let t1 ≤ t2 ≤ . . . ≤ tn+m be anon-decreasing sequence of knots, and for each non-empty interval (ti, ti+1) let ρi bea given non-negative number. Tension B-splines of order k are defined recursively[2, 8]Bj,k :=∫ xtjBj,k−1(y)dyσj,k−1−∫ xtj+1Bj+1,k−1(y)dyσj+1,k−1, j = m− k + 1, . . . n, (1)whereσj,k−1 :=∫ tj+k−1tjBj,k−1(y)dy .The recursion starts with functions Bi,2 defined byBi,2(x) :=⎧⎨⎩sinh(ρi(x− ti))/sinh(ρi∆i) ; ti ≤ x < ti+1 ,sinh(ρi+1(ti+2 − x))/sinh(ρi+1∆i+1) ; ti+1 ≤ x < ti+2 ,0 ; elsewhere,∆i = ti+1 − ti .In [2] it is proved that B-splines Bi,k have the following properties:1. suppBi,k = [ti, ti+k] (local support property)2. Bi,k(x) > 0 for ti < x < ti+k (positivity)3.∑i Bi,k ≡ 1, k ≥ 3 (partition of unity)4. Bi,k has k − 1− l continuous derivatives at a knot of multiplicity l.Moreover, B-splines Bj,k, j = m−k+1, . . . n, are linearly independent and theyform a basis for the space Tk = span {Bm−k+1,k, . . . , Bn,k}. Tk is a space of tensionsplines of order k, with smoothness at knots ti defined by multiplicity of the knots.3. B-spline representationNow we consider a local representation of the B-splines on the subinterval (ti, ti+1).We suppose that ti = ti+1 and that ρi is given.Evaluation of tension splines 75Restriction of each B-spline on subinterval (ti, ti+1) can be represented in theformBj,k(x) = βj,k−2,iϕk(x− ti) + βj,k−1,iϕk(ti+1 − x) +k−3∑l=0βj,kl,i(x− ti∆i)l, (2)where ϕk(x) := ϕk(x; ρi,∆i) is defined byϕk(x; ρ,∆) =⎧⎨⎩xk−1(k−1)!∆ , ρ = 0 ,Fk(ρx)−Pk(ρx)ρk−2 sinh(ρ∆) , ρ > 0,(3)andFk(t) ={sinh(t), k = 2l,cosh(t), k = 2l + 1, and Pk(t) =⎧⎨⎩∑l−2j=0t2j+1(2j+1)! , k = 2l,∑l−1j=0t2j(2j)! , k = 2l + 1.(4)We may note that ∫ u0ϕk(y)dy = ϕk+1(u) , u < ∆i . (5)Now, we are going to find a representation of∫ xtjBj,k(y)dy/σj,k on subinterval(ti, ti+1). Because of the local support property,∫ xtjBj,k(y)dy/σj,k = 0 for ti+1 < tj ,and it is equal to 1 for ti > tj+k. Otherwise, for x ∈ (ti, ti+1),∫ xtjBj,k(y)dyσj,k=(∫ xti(βj,k−2,iϕk(y − ti) + βj,k−1,iϕk(ti+1 − y))dy+∫ xtik−3∑l=0βj,kl,i(y − ti∆i)ldy +∫ titjBj,k(y)dy)/σj,k . (6)Taking into account that∫ xtiϕk(y − ti)dy = ϕk+1(x− ti),∫ xtiϕk(ti+1 − y)dy = ϕk+1(∆i)− ϕk+1(ti+1 − x) ,∫ xti(y − ti∆i)ldy =∆il + 1(x− ti∆i)l+1and definingσij,k =∫ titjBj,k(y)dy ,from (6) we obtain∫ xtjBj,k(y)dyσj,k=1σj,k(βj,k−2,iϕk+1(x− ti)− βj,k−1,iϕk+1(ti+1 − x)76 I. Berosˇ and M.Marusˇic´+ σij,k + βj,k−1,iϕk+1(∆i) + ∆ik−2∑l=1βj,kl−1,il(x− ti∆i)l ). (7)Using (1) and (7), we establish a relation among coefficients βj,k+1 of B-splineBj,k+1 and coefficients βj,k and βj+1,k of Bj,k and Bj+1,k, respectively:βj,k+1−2,i =βj,k−2,iσj,k− βj+1,k−2,iσj+1,kβj,k+1−1,i =βj+1,k−1,iσj+1,k− βj,k−1,iσj,kβj,k+10,i =σij,k + ϕk+1(∆i)σj,k− σij+1,k + ϕk+1(∆i)σj+1,k(8)βj,k+1l,i = ∆iβj,kl−1,ilσj,k−∆iβj+1,kl−1,ilσj+1,k, l = 1, . . . , k − 2 ,for k ≥ 3 and j = i− k + 1, . . . , i. The recursion starts withβi−1,2−2,i = 0, βi−1,2−1,i = 1, βi,2−2,i = 1 and βi,2−1,i = 0 . (9)Coefficients σij,k and σj,k are obtained fromσij,k =i−1∑l=j∫ tl+1tlBj,k(y)dy=i−1∑l=j((βj,k−2,l + βj,k−1,l)ϕk+1(∆l) + ∆lk−3∑s=0βj,ks,ls + 1),noting that σj,k = σj+kj,k .4. Computing functions ϕkFunctions ϕk, defined by (3) depend on three parameters, and for each parameterwe only know that it is nonnegative. It is more convenient to use functionsϕ˜k(p, t) =⎧⎨⎩tk−1(k−1)! , p = 0 ,Fk(pt)−Pk(pt)pk−2 sinh(p) , p > 0 .(10)By substitutingp = ρ∆ and t =x∆we obtain the relation:ϕk(x; ρ,∆) = ∆k−2ϕ˜k(ρ∆,x∆) . (11)Evaluation of tension splines 77Function ϕ˜k depends on two parameters. Parameter t is from interval [0, 1] and pis nonnegative.Computing hyperbolic-like functions ϕ˜k by using explicit formula (10) is verysensitive. For p or pt small (10) is numerically unstable, while for p or pt large (e.g.over 700 for double precision arithmetic), machine overflow occurs. So we have toconsider evaluation of ϕ˜k in the different situations. All calculation presented inthe following text is done in FORTRAN 77 for HP 9000.4.1. Case when p is smallWhen p is close to 0, there are two possible sources of error in evaluation of ϕ˜k by(10). The first is a loos of significant digits in subtraction Fk(pt)− Pk(pt) and thesecond is a possible underflow in evaluation of pk−2 sinh p. Figure 1 shows relativeerrors in evaluation of ϕ˜6 in this domain (small p). For example, we can see that(10) is inappropriate for p ≤ 0.5.Figure 1. Relative error for evaluating ϕ˜6 by (10), on [0, 0.5]× [0, 1], in doubleprecision arithmetic.To solve this problem, we expand the numerator and the denominator in Taylorseries. For k even we obtain:ϕ˜k(p, t) =Fk(pt)− Pk(pt)pk−2 sinh(p)=∑∞j=l−1(pt)2j+1(2j+1)!pk−1 sinh(p)p== tk−11(k−1)! +∑∞j=l(pt)2(j−l+2)(2j+1)!1 +∑∞j=1p2j(2j+1)!= tk−11(k−1)! +∑∞j=1(pt)2j(k+2j−1)!1 +∑∞j=1p2j(2j+1)!.78 I. Berosˇ and M.Marusˇic´The same formula holds for k odd, so ϕ˜k is of the formϕ˜k(p, t) = tk−11(k−1)! +∑∞j=1(pt)2j(k+2j−1))!1 +∑∞j=1p2j(2j+1)!, (12)for all k. The above representation does not have singularity for p = 0.To approximate (12), we use finite instead of infinite sums. In that case, (12) isapproximated byGmk (p, t) = tk−11(k−1)! +∑mj=1 (pt)2j/(k + 2j − 1)!1 +∑mj=1 p2j/(2j + 1)!. (13)We define numbers pmk :pmk := max{p ≥ 0 :∣∣∣∣ ϕ˜k(p, t)−Gmk (p, t)ϕ˜k(p, t)∣∣∣∣ <  , t ∈ [0, 1]}, (14)defining the region where a relative error in evaluation of ϕ˜k is less than . Thenumbers pmk are found experimentally calculating ϕ˜k(p, t) according to formula (10)in quadri precision arithmetic, while Gmk is calculated in double precision arithmetic.We have used tolerance  = 10−15. Table 1 shows values of pmk for different k andm. It is noteworthy that pmk almost does not depend on k.mk 4 5 6 7 8 92 0.08 0.21 0.41 0.68 1.02 1.413 0.08 0.21 0.41 0.68 1.02 1.414 0.08 0.21 0.41 0.68 1.02 1.425 0.08 0.21 0.41 0.68 1.02 1.426 0.08 0.21 0.41 0.68 1.02 1.427 0.10 0.21 0.41 0.69 1.02 1.428 0.10 0.21 0.40 0.68 1.01 1.40Table 1. Numbers pmk defined by (14), for  = 10−154.2. Case when pt is smallFigures 1 and 2 show a large relative error in evaluation of ϕ˜k not only for small pbut also for p relatively large and t small. When product pt is small, but p relativelylarge (e.g. p ≥ pmk ), we do not have a problem with underflow in the denominator(as in 4.1), but it is still possible to loose significant digits in evaluation of thenumerator. Figure 2 shows relative errors in evaluation of ϕ˜6 in this situation.Evaluation of tension splines 79Figure 2. Relative error for evaluating ϕ˜6 by (10), on [1, 10]× [0, 0.1] in doubleprecision arithmetic.To avoid a loss of significant digits, we expand the numerator in Taylor series:ϕ˜k(p, t) = 2e−pptk−1∑∞j=0(pt)2j(k+2j−1)!(1− e−p)(1 + e−p) . (15)Substitutionpk−2 sinh p = pk−2(1− e−p)(1 + e−p)2e−p,is used to avoid a potential overflow in the denominator for large p.As in (12), we use finite instead of infinite sum to approximate (15). In thatcase, (15) is approximated byG˜mk = 2e−pptk−1∑mj=0 (pt)2j/(k + 2j − 1)!(1− e−p)(1 + e−p) . (16)Region where a relative error for G˜mk is less than  is described by pt ≤ umk , whereumk = max{pt > 0 :∣∣∣∣∣ ϕ˜k(p, t)− G˜mk (p, t)ϕ˜k(p, t)∣∣∣∣∣ <  , t ∈ (0, 1], p > 0}. (17)As the numbers pmk (14), the numbers umk are found by a numerical experimentfor  = 10−15.mk 4 5 6 7 8 92 0.07 0.20 0.40 0.67 1.00 1.403 0.09 0.24 0.47 0.77 1.14 1.584 0.10 0.27 0.52 0.87 1.27 1.745 0.12 0.31 0.59 0.97 1.40 1.926 0.13 0.34 0.65 1.06 1.53 2.087 0.14 0.37 0.71 1.14 1.65 2.238 0.16 0.42 0.79 1.25 1.81 2.41Table 2. Numbers umk , defined by (17), for  = 10−154.3. Case when p is largeWhen p is large, we have the problem of machine overflow. For example, in calcu-lating ϕ˜6 by (10) overflow occurs for p ≥ 710 in double precision arithmetic, andfor p ≥ 89 in single precision arithmetic.80 I. Berosˇ and M.Marusˇic´For large p, we rewrite (10) asϕ˜k(p, t) =ept−(−1)ke−pt2 − Pk(pt)pk−2 ep−e−p2= e−p(1−t)1− 2e−ptPk(pt)− (−1)ke−2ptpk−2(1− e−p)(1 + e−p) . (18)Because oflimx→∞xke−x = 0 , for all k ,there exist numbers Vk such that2Pk(y)e−y < ε for y > Vk,where ε is machine precision. Since e−p < 2Pk(pt)e−pt holds for pt > Vk, ϕ˜k(p, t)can be approximated by an asymptotic formula:Hk(p, t) =e−p(1−t)pk−2. (19)By using (19) we avoid the problem of machine overflow.Numbers Vk are found experimentally. We found the smallest integer y satisfyingfl(1 + 2Pk(y)e−y) = 1. (20)Here fl(C) stands for “floating point result of C”. Values for Vk in double precisionarithmetic (ε ≈ 10−16) are presented in Table 3.k 2 3 4 5 6 7 8Vk 19 38 42 45 48 50 54Table 3. Numbers Vk, defined by (20), for double precision arithmetic4.4. Regular caseOutside the regions defined by 4.1–4.3 we use formula (10).Evaluation of tension splines 81Figure 3. Regions for evaluating of ϕ˜k: A – formula (13), B – formula (16), C –formula (19), D – formula (10).To summarize, evaluation of ϕ˜k(p, t) is divided in four different cases:A) When p ≤ pmk , for arbitrary m, we approximate ϕ˜k(p, t) by Gmk (p, t) (13).B) When p > pmk and pt ≤ ulk, for arbitrary m and l, we approximate ϕ˜k(p, t) byG˜lk(p, t) (16).C) When pt > Vk, we approximate ϕ˜k(p, t) by Hmk (p, t) (19).D) Otherwise, we evaluate ϕ˜k(p, t) by (10).AcknowledgmentsSupported by grant 037011 by Croatian Ministry of Science and Technology.References[1] J. E.Flaherty, W.Mathon, Collocation with polynomial and tension splinesfor singularly-perturbed boundary value problems, SIAM J. Sci. Stat. Comput.1(1980), 260–289.[2] P.E.Koch, T. Lyche, Construction of exponential tension B-spline of ar-bitrary order, in: Curves and Surfaces, (P. J. Laurent, A. LeMehaute andL. L. Schumaker, Eds.), Academic Press, N.Y., 1991, 255–258.[3] B. I.Kvasov, Local bases for generalized cubic splines, Russian J. Numer. Anal.Math. Modelling 10(1995), 49–80.[4] M.Marusˇic´, Stable calculation by splines in tension, Grazer Math. Ber.328(1996), 65–76.[5] M.Marusˇic´, M.Rogina, A collocation method for singularly perturbed two-point boundary value problems with splines in tension, Adv. Comput. Math.6(1996), 65–76.[6] M.Marusˇic´, M.Rogina, B-splines in tension, In: Proceedings of VII Confer-ence on Applied Mathematics, (R. Scitovski, Ed.), Osijek, 1989, 129–134.[7] B. J.McCartin, Theory of exponential splines J. Approx. Theory 66(1991),1–23.[8] M.Rogina, Basis of splines associated with some singular differential operators,BIT 32(1992), 496–505.[9] D.G. Schweikert, An interpolating curve using a spline in tension, J. Math.Physics 45(1966), 312–317.[10] H. Spa¨th, Exponential spline interpolation, Computing 4(1969), 225–233.

Evaluation of tension splines

Abstract

Similar works

Full text

Available Versions

Hrčak - Portal of scientific journals of Croatia