Monotone cubic spline interpolation for functions with a strong gradient
aa r X i v : . [ m a t h . NA ] F e b Monotone cubic spline interpolation for functions with a strong gradient
Francesc Ar`andiga a , Antonio Baeza a , Dionisio F. Y´a˜nez a a Departament de Matem`atiques. Universitat de Val`encia (EG) (Spain)
Abstract
Spline interpolation has been used in several applications due to its favorable properties regarding smoothnessand accuracy of the interpolant. However, when there exists a discontinuity or a steep gradient in the data, someartifacts can appear due to the Gibbs phenomenon. Also, preservation of data monotonicity is a requirement insome applications, and that property is not automatically verified by the interpolator. In this paper, we studysufficient conditions to obtain monotone cubic splines based on Hermite cubic interpolators and propose different waysto construct them using non-linear formulas. The order of approximation, in each case, is calculated and severalnumerical experiments are performed to contrast the theoretical results.
Keywords:
Monotonicity, Cubic Hermite Interpolants, Cubic Spline Interpolants, Non-linear computation ofderivatives.
1. Introduction and review: Hermite cubic interpolation
Approximation techniques are used in applications as design of curves, surfaces, robotics, creation of pieces inindustry and many others due to the fact that they present certain regularity properties. In particular, Hermiteinterpolatory polynomials have been developed to obtain interpolants of class C that have been applied, for example,to the numerical solution of differential equations (see [1, 9]). We consider the problem of piecewise cubic Hermiteinterpolation, that can be stated as follows: let x < x < . . . < x n be a partition of the interval [ x , x n ] and let f i = f ( x i ) be the values of a certain function at the knots. Given approximate values of the first derivative of f at the knots { x i } ni =1 , denoted by { ˙ f i } ni =1 , construct a piecewise cubic polynomial function P ( x ) , conformed by n − P i ( x ) defined on the ranges [ x i , x i +1 ] , that satisfy P ( x i ) = f i , P ( x i +1 ) = f i +1 ,P ′ ( x i ) = ˙ f i , P ′ ( x i +1 ) = ˙ f i +1 . (1)We will use the following notation: the undivided differences of a function f are denoted by ∆ f i = f i +1 − f i and m i = ∆ f i /h i denotes its divided differences, where h i = x i +1 − x i are the mesh spacings and ˆ h = max i =1 ,...,n − ( h i ).The i -th polynomial P i ( x ) , x ∈ [ x i , x i +1 ] (see [7] for details), has the form P i ( x ) = c i + c i ( x − x i ) + c i ( x − x i ) + c i ( x − x i ) ( x − x i +1 ) , (2)where: c i := f i , c i := f [ x i , x i ] = ˙ f i , c i := f [ x i , x i , x i +1 ] = ( m i − ˙ f i ) /h i ,c i := f [ x i , x i , x i +1 , x i +1 ] = ( ˙ f i +1 + ˙ f i − m i ) /h i . Hence, a procedure to compute { ˙ f i } defines an algorithm for constructing a cubic Hermite interpolant.In some problems, it is required that the interpolant employed preserves the monotonicity of the data. Thisproblem has been tackled in the literature (see, e.g., [12, 19, 6, 8, 4]) leading to several options for monotonic Hermiteinterpolation. In the remaining of this section we cite some known results dealing with conditions for a cubic Hermiteinterpolants to be monotonicity preserving and about its accuracy. ⋆ This research has been supported by Spanish MINECO projects MTM2017-83942-P.
Email addresses: [email protected] (Francesc Ar`andiga), [email protected] (Antonio Baeza), [email protected] (Dionisio F.Y´a˜nez)
Preprint submitted to Elsevier Received: date / Accepted: date heorem 1.1. (Necessary conditions for monotonicity.) Let P i be a monotone cubic Hermite interpolant of the data { ( x i , f i , ˙ f i ) , ( x i +1 , f i +1 , ˙ f i +1 ) } . Then: sign ( ˙ f i ) = sign ( ˙ f i +1 ) = sign ( m i ) . (3) Furthermore, if m i = 0 then P i is monotone (constant) if and only if ˙ f i = ˙ f i +1 = 0 . Theorem 1.2. (Sufficient conditions for monotonicity [11].) Let I i = [ x i , x i +1 ] and P i be a cubic Hermite interpolantof the data { ( x i , f i , ˙ f i ) , ( x i +1 , f i +1 , ˙ f i +1 ) } , and let α i := ˙ f i /m i , β i := ˙ f i +1 /m i . If ≤ α i , β i ≤ , (4) then the resulting cubic Hermite interpolant (2) is monotone on I i . The following theorem provides a global result (see [13]). This is a generalization of Theorem 1.2.
Theorem 1.3. ([13]) With the notation of Theorem 1.2, if for all i , ≤ i ≤ n − , | ˙ f i | ≤ | m i − | , | m i | ) (5) then (2) is monotone in each [ x i , x i +1 ] , ≤ i ≤ n − . Finally, we show a more general theorem proved in [11]:
Theorem 1.4. (Sufficient conditions for monotonicity) Let I i = [ x i , x i +1 ] and P i be a cubic Hermite interpolant ofthe data { ( x i , f i , ˙ f i ) , ( x i +1 , f i +1 , ˙ f i +1 ) } , and let α i := ˙ f i /m i , β i := ˙ f i +1 /m i . If one of the following conditions aresatisfied ≤ α i + β i ≤ ,α i + α i ( β i −
6) + ( β i − < , (6) ∀ i = 1 , . . . , n − then the resulting cubic Hermite interpolant (2) is monotone on I i . There exist many methods in the literature that deal with the problem of computing approximate derivative valuesin a way such that the resulting polynomials keep high (ideally, maximal) order of approximation and at the sametime produce monotonicity-preserving reconstructions. They all face the problem that in order to ensure high orderaccuracy, the monotonicity–preserving property is lost and conversely. In [11] it is proved that the use of non-lineartechniques is necessary to obtain third-order accurate interpolants ([3]) using the following lemma:
Lemma 1.5. ([7]) Let us assume that f ( x ) is smooth. If ˙ f i = f ′ ( x i ) + O ( h q i ) and ˙ f i +1 = f ′ ( x i +1 ) + O ( h q i +1 ) , thenthe piecewise cubic Hermite interpolant (2) satisfies, on the interval [ x i , x i +1 ] , P i ( x ) = f ( x ) + O ( h q ) where q = min(4 , q i + 1 , q i +1 + 1) . (7)In addition to monotonicity preservation, in some applications some regularity is demanded. In order to obtain C approximations we introduce cubic spline interpolation in Section 2. Also, we determine sufficient conditions to obtainmonotone spline cubic interpolants using the theorems presented above. In Section 3 we construct new monotoneinterpolants and study their properties regarding the order of approximation. Some numerical experiments are shownin Section 4 in order to confirm the properties of the proposed algorithms. Finally, some conclusions are presented inSection 5
2. Cubic spline interpolation
In this section, we construct a cubic spline P i ( x ) verifying the conditions P ( k ) i ( x i +1 ) = P ( k ) i +1 ( x i +1 ) , k = 0 , , , (8) P ′ ( x ) = ˙ f := f ′ ( x ) = f ′ , (9) P ′ n − ( x n ) = ˙ f n := f ′ ( x n ) = f ′ n , (10)with i = 1 , . . . , n − k = 0 , k = 2 will be used to obtain the appropriate approximations to thevalues of the first derivatives. 2 .1. Spline cubic interpolation from Hermite cubic form As indicated above we start from Eq. (2) and we impose that: P ′′ i ( x i +1 ) = P ′′ i +1 ( x i +1 ) , i = 1 , . . . , n − , (11)thus, we have that: P ′′ l ( x ) = 2 c l + 4 c l ( x − x l ) + 2 c l ( x − x l +1 ) , l = i, i + 1 . Then, using (3) we get P ′′ i ( x i +1 ) = 2 c i + 4 c i h i = 2 m i − ˙ f i h i ! + 4 ˙ f i +1 + ˙ f i − m i h i ! = 2 ˙ f i h i + 4 ˙ f i +1 h i − m i h i ,P ′′ i +1 ( x i +1 ) = 2 c i +13 + 2 c i +14 ( x i +1 − x i +2 ) = 2 m i +1 h i +1 − f i +1 h i +1 − ˙ f i +2 + ˙ f i +1 − m i +1 h i +1 ! = − f i +1 h i +1 − f i +2 h i +1 + 6 m i +1 h i +1 . By Eq. (11), we obtain2 ˙ f i h i + 4 ˙ f i +1 h i − m i h i = − f i +1 h i +1 − f i +2 h i +1 + 6 m i +1 h i +1 ⇒ ˙ f i h i + 2 (cid:18) h i + 1 h i +1 (cid:19) ˙ f i +1 + ˙ f i +2 h i +1 = 3 (cid:18) m i h i + m i +1 h i +1 (cid:19) ⇒ h i +1 h i + h i +1 ˙ f i + 2 ˙ f i +1 + h i h i + h i +1 ˙ f i +2 = 3 (cid:18) h i +1 h i + h i +1 m i + h i h i + h i +1 m i +1 (cid:19) . (12)If we take the boundary conditions given in Eq. (9)–(10), i. e. ˙ f = f ′ and ˙ f n = f ′ n , we have the following system f + µ ˙ f = b ,λ i ˙ f i + 2 ˙ f i +1 + µ i ˙ f i +2 = b i +1 , i = 2 , . . . , n − ,λ n − ˙ f n − + 2 ˙ f n − = b n − , (13)where λ i = h i +1 h i + h i +1 , µ i = h i h i + h i +1 , λ i + µ i = 1, i = 1 , . . . , n −
1, and b = 3 ( λ m + µ m ) − λ f ′ ,b i = 3 ( λ i − m i − + µ i − m i ) , i = 3 , . . . , n − ,b n − = 3 ( λ n − m n − + µ n − m n − ) − µ n − f ′ n . Thus, the system obtained is A ˙ F = B, (14)with A = µ λ µ λ µ . . . . . . . . . λ n − µ n − λ n − , ˙ F = ˙ f ˙ f ˙ f ...˙ f n − ˙ f n − ˙ f n − , B = λ m + µ m ) − λ f ′ λ m + µ m )3( λ m + µ m )...3( λ n − m n − + µ n − m n − )3( λ n − m n − + µ n − m n − )3( λ n − m n − + µ n − m n − ) − µ n − f ′ n . (15)The equality λ i + µ i = 1 implies that the matrix A is irreducibly diagonally dominant and hence non-singular. Inorder to calculate the order of the approximate derivative values computed by solving (14), we prove the followinglemma and theorem using the ideas presented in [18] (Eq. 2.4.2.14).3 emma 2.1. ([18]) Let be < m ∈ N , ≤ µ i , λ i ≤ with ≤ i ≤ m such that λ i + µ i = 1 i = 1 , . . . , m, and A ∈ R m × m defined as: A := µ λ µ λ µ . . . . . . . . . λ m − µ m − λ m . (16) Given w ∈ R m , if z ∈ R m solves Az T = w then || z || ∞ ≤ || w || ∞ . Proof.
Let i be such that | z i | = || z || ∞ , then || w || ∞ ≥ | w i | = | λ i +1 z i − + 2 z i + µ i +1 z i +1 |≥ | z i | − λ i +1 | z i − | − µ i +1 | z i +1 |≥ | z i | − ( λ i +1 + µ i +1 ) | z i | = | z i | = || z || ∞ (17) (cid:4) Lemma 2.2.
Let us assume that f ( x ) ∈ C ([ x , x n ]) and let L > be such that such that | f (4) ( x ) | ≤ L , for all x ∈ [ x , x n ] . If there exists K > such that ˆ h/h j ≤ K for all j = 1 , . . . , n − and R ( i ) := 3 λ i − m i − + 3 µ i − m i − λ i − f ′ i − − f ′ i − µ i − f ′ i +1 , ≤ i ≤ n − , (18) with λ i , m i , µ i , ≤ i ≤ n − , previously defined, then: | R ( i ) | ≤ (cid:18) K + K
16 + 1 (cid:19) L ˆ h = O (ˆ h ) , ≤ i ≤ n − . (19) Proof.
Let be 2 ≤ i ≤ n −
1, using Taylor’s expansions we have that there exist τ ij ∈ [ x , x n ] , j = 1 . . . , m i − = f ( x i − + h i − ) − f ( x i − ) h i − = f ′ i − + h i − f ′′ i − + h i − f (3) i − + h i − f (4) ( τ i ) m i = (cid:18) ( h i − + h i ) h i f ′ i − + ( h i − + h i ) h i f ′′ i − + ( h i − + h i ) h i f (3) i − + ( h i − + h i ) h i f (4) ( τ i ) (cid:19) − h i − h i f ′ i − − h i − h i f ′′ i − − h i − h i f (3) i − − h i − h i f (4) ( τ i )= f ′ i − + (2 h i − + h i )2 f ′′ i − + (cid:0) h i − + 3 h i − h i + h i − (cid:1) f (3)1 + ( h i − + h i ) h i f (4) ( τ i ) − h i − h i f (4) ( τ i ) f ′ i = f ′ i − + h i − f ′′ i − + h i − f (3) i − + h i − f (4) ( τ i ) f ′ i +1 = f ′ i − + ( h i − + h i ) f ′′ i − + ( h i − + h i ) f (3) i − + ( h i − + h i ) f (4) ( τ i ) (20)4hen by (20) we have: | R ( i ) | = | λ i − m i − + 3 µ i − m i − λ i − f ′ i − − f ′ i − µ i − f ′ i +1 | = | λ i − (3 m i − − f ′ i − ) + µ i − (3 m i − f ′ i +1 ) − f ′ i | = (cid:12)(cid:12)(cid:12)(cid:12) h i − + h i ) ( (cid:16) h i f ′ i − + (3 h i − h i ) f ′′ i − + ( h i − h i ) f (3) i − (cid:17) ++ (cid:16) h i − f ′ i − + (cid:0) h i − + h i − h i (cid:1) f ′′ i − + (cid:0) h i − + 2 h i − h i (cid:1) f (3) i − (cid:17) ) − (cid:18) f ′ i − + h i − f ′′ i − + h i − f (3) i − (cid:19) + h i − h i h i − + h i ) f (4) ( τ i ) + h i − ( h i − + h i ) h i f (4) ( τ i ) − h i − h i ( h i − + h i ) f (4) ( τ i ) − h i − ( h i − + h i ) f (4) ( τ i ) − h i − f (4) ( τ i ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:18) h i − h i h i − + h i ) + h i − ( h i − + h i ) h i + h i − h i ( h i − + h i ) + h i − ( h i − + h i ) h i − (cid:19) L ≤ K ˆ h
16 + K ˆ h + K ˆ h
16 + 2ˆ h h ! L = (cid:18) K + K
16 + 1 (cid:19) L ˆ h = O (ˆ h ) . (cid:4) Theorem 2.3. ([18]) Let us assume that ˆ h < , f ( x ) ∈ C ([ x , x n ]) and let L > be such that | f (4) ( x ) | ≤ L , for all x ∈ [ x , x n ] . F ′ = [ f ′ ( x ) , . . . , f ′ ( x n − )] T , ˙ F is the solution of system (14) and there exist K > , such that ˆ h/h j ≤ K for all j = 1 , . . . , n − then: || ˙ F − F ′ || ∞ ≤ O (ˆ h ) . Proof.
We define r = A ( ˙ F − F ′ ) = A ˙ F − AF ′ = B − AF ′ , by Lemma 2.2 we have: | r | = | b − f ′ − µ f ′ | = | λ m + 3 µ m − λ f ′ − f ′ − µ f ′ | = | R (2) | ≤ (cid:18) K + K
16 + 1 (cid:19) L ˆ h = O (ˆ h ) , | r i − | = | b i − λ i − f ′ i − − f ′ i − µ i − f ′ i +1 | = | λ i − m i − + 3 µ i − m i − λ i − f ′ i − − f ′ i − µ i − f ′ i +1 | = | R ( i ) |≤ (cid:18) K + K
16 + 1 (cid:19) L ˆ h = O (ˆ h ) , ≤ i ≤ n − , | r n − | = | b n − − λ n − f ′ n − − f ′ n − | = | λ n − m n − + µ n − m n − ) − µ n − f ′ n − λ n − f ′ n − − f ′ n − | = | R ( n − | ≤ (cid:18) K + K
16 + 1 (cid:19) L ˆ h = O (ˆ h ) . Finally, by Lemma 2.1, we obtain: || r || ∞ = || A ( ˙ F − F ′ ) || ∞ ⇒ || ( ˙ F − F ′ ) || ∞ ≤ || r || ∞ ≤ O (ˆ h ) . (cid:4) In the case of equally-spaced grids [16] the estimate in Theorem 2.3 can be improved to || ˙ F − F ′ || ∞ ≤ O (ˆ h ) (21) Corollary 2.4.
The cubic Hermite interpolant (2) obtained using the approximations to the derivatives solving thesystem (14) has order of accuracy O (ˆ h ) .Proof. This is direct by Lemma 1.5 and Theorem 2.3. (cid:4)
Now, we indicate the following conditions on the values ˙ f i , i = 2 , . . . , n −
2, in order to obtain a monotoneinterpolator in each interval [ x i , x i +1 ]. We can prove the following result using Theorem 1.3.5 heorem 2.5. If for all ≤ i ≤ n − , (cid:12)(cid:12)(cid:12) ˙ f i (cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) n − X j =1 A − ij b j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ | m i − | , | m i | ) (22) then the resulting cubic Hermite interpolant (2) is monotone. These conditions, in general, are not satisfied.If there exist any points where the computed approximations of the derivative do not satisfy the sufficient conditionswe will analyze two possibilities:1. Change some of the derivative approximations, obtained by solving (14), by other values that produce a monotoneinterpolant. Hence, we define for i = 2 , . . . , n − f Mi := ( P n − j =1 A − ij b j , if | P n − j =1 A − ij b j | ≤ | m i − | , | m i | );˜ f i , otherwise, (23)where ˜ f i is a value calculated using non-linear techniques that will be explained in Section 3. We only modifythe approximations of the derivative values at the points where monotonicity constraints were not satisfied. Inthis way the interpolant keeps the maximum order in every interval. With this method, the regularity is reducedto C in a neighborhood of each point where the approximation of the derivative is modified. We analyze thiscase in Section 2.2.2. Another possibility is to try to keep C regularity in the complete interval [ x , x n ] except at the single pointswhere the approximations of the derivatives have been changed. Our proposal is to change the values wheremonitonicity was lost as in the previous case, rewrite system (14) but eliminating the modified points, and solveit. Afterwards, we again study if the new values satisfy the monotonicity conditions and repeat the process. Wewill prove in Section 2.3 that the order is lost in a neighborhood of the conflicting points, but is conserved atthe rest. Assume that there exists a point x i , with 1 < i < n , where the approximation to the derivative does not satisfythe conditions of the Theorem 1.3. In that case, we change the value ˙ f i by another value ˜ f i . As a result the followingequalities are not necessarily satisfied: P ′′ i − ( x i − ) = P ′′ i − ( x i − ) ,P ′′ i − ( x i ) = P ′′ i ( x i ) ,P ′′ i ( x i +1 ) = P ′′ i +1 ( x i +1 ) . Thus, the regularity is C in all points excepted at x i + l with l = − , ,
1. Finally, by Lemma 1.5 the order is 4 exceptin the intervals I i + l with l = − , ,
1. We recapitulate these results in the following proposition.
Proposition 2.6.
Let us assume that f ( x ) ∈ C ([ x , x n ]) and let L > be such that | f (4) ( x ) | ≤ L , for all x ∈ [ x , x n ] .Let p i , T > and A , ˙ F , B be as defined in Eq. (15) , which satisfy that A ˙ F = B . If we define ˙ F i = (cid:26) ˜ f i , i = i ;˙ f i , i = i , such that | f ′ i − ˜ f i | = T · ˆ h p i and there exists K > such that ˆ h/h i ≤ K for all i = 1 , . . . , n then: | ˙ f i − f ′ i | = (cid:26) O (ˆ h p i ) , i = i ; O (ˆ h ) , i = i . (24) Also, the cubic spline interpolator defined in Eq. (2) using F i as an approximation of the values of the first derivatives,has C regularity except at the points x i + l , l = − , , . .3. Monotone spline with maximum regularity As a second option we replace the values for which the approximate first derivative does not satisfy the conditionsin Theorem 1.3, and recalculate the remaining values by rewriting system (14) in a way such that the equationscorresponding to the changed values are removed from the system. Thus, we suppose that the approximation to thefirst derivative does not satisfy the sufficient conditions at i with 1 < i < n . Consequently, we calculate ˙ f Mi = ˜ f i and define A i − ∈ R ( i − × ( i − and A i +0 ∈ R ( n − − i ) × ( n − − i ) by: A i − = µ λ µ λ µ . . . . . . . . . λ i − µ i − λ i − , A i +0 = µ i λ i +1 µ i +1 λ i +2 µ i +2 . . . . . . . . . λ n − µ n − λ n − , (25)and the vectors: ˙ F i − = ˙ f ˙ f ˙ f ...˙ f i − ˙ f i − ˙ f i − , B i − = b b b ... b i − b i − b i − ˙ F i +0 = ˙ f i +1 ˙ f i +2 ˙ f i +3 ...˙ f n − ˙ f n − ˙ f n − , B i +0 = b i +1 b i +2 b i +3 ... b n − b n − b n − , (26)being: b = 3( λ m + µ m ) − λ f ′ ,b n − = 3( λ n − m n − + µ n − m n − ) − µ n − f ′ n ,b i − = 3( λ i − m i − + µ i − m i − ) − µ i − ˜ f i ,b i +1 = 3( λ i m i + µ i m i +1 ) − λ i ˜ f i ,b i − = 3( λ i − m i − + µ i − m i ) , for 2 < i < n − , i = i + l, l = − , , . (27)With these variables we can rewrite the new system as A i ˙ F i = B i where: A i = " A i − A i +0 , ˙ F i = " ˙ F i − ˙ F i +0 , B i = " B i − B i +0 . (28)We will prove that there exists a set of intervals around of x i where the order of the approximation to the firstderivative is affected because of the modification of ˙ f i but is maintained at the points that are sufficiently separatedfrom the discontinuity.We adapt the results obtained in [14, 15, 16]. For this, we divide the system A i ˙ F i = B i in two subsystemsand analyze them separately. Each subsystem is similar to the system obtained to construct a spline with differentboundary conditions.The following result [14] provides a bound of the elements of the inverse. Lemma 2.7.
Let be < m ∈ N , ≤ µ i , λ i ≤ with ≤ i ≤ m such that λ i + µ i = 1 i = 1 , . . . , m, and A ∈ R m × R m defined by 16. If the elements of A − are denoted by A − ij then | A − ij | ≤ · −| i − j | , ≤ i, j ≤ m. emark 2.1. In the case of uniform grid, i.e., when h i = h i − , for all i = 2 , . . . , n , then the bound can be improved.In [14], it is proved that: | A − ij | ≤ · (2 + √ −| i − j | , ≤ i, j ≤ m. (29) Proposition 2.8.
Let us assume that ˆ h < , f ( x ) ∈ C ([ x , x i ]) and let L > be such that | f (4) ( x ) | ≤ L , forall x ∈ [ x , x i − ] . Let p i , T > and A i − , ˙ F i − , B i − be as defined in Eqs. (25) , (26) and (27) , which satisfy | f ′ i − ˜ f i | = T · ˆ h p i and A i − ˙ F i − = B i − . If i > − log (ˆ h ) and there exists K > , such that ˆ h/h i ≤ K for all i = 1 , . . . , i − then:and ˆ h < There exists an integer l < i − , such that: | ˙ f i − f ′ i | = (cid:26) O (ˆ h min(3 ,p i +1) ) , ≤ i ≤ l ; O (ˆ h p i ) , l < i ≤ i − . (30) Proof.
We define r i − = (( r i − ) , . . . , ( r i − ) i − ) T as: r i − = A i − ( ˙ F i − − F ′ i − ) = A i − ˙ F i − − AF ′ i − = B i − − AF ′ i − . Let i be such that 2 < i ≤ i −
2. By Lemma 2.2: | ( r i − ) i | = | b i − λ i − f ′ i − − f ′ i − µ i − f ′ i +1 | = | λ i − m i − + 3 µ i − m i − λ i − f ′ i − − f ′ i − µ i − f ′ i +1 | = | R ( i ) | ≤ (cid:18) K + K
16 + 1 (cid:19) L ˆ h = O (ˆ h ) . (31)The result is proved analogously for i = 2. In the case i = i − | ( r i − ) i − | = | b i − − λ i − f ′ i − − f ′ i − | = | b i − − λ i − f ′ i − − f ′ i − + ( − µ i − f ′ i + µ i − f ′ i ) | = (cid:12)(cid:12)(cid:12)(cid:0) λ i − m i − + µ i − m i − ) − λ i − f ′ i − − f ′ i − − µ i − f ′ i (cid:1) + µ i − ( f ′ i − ˜ f i ) (cid:12)(cid:12)(cid:12) ≤| R ( i − | + µ i − (cid:12)(cid:12)(cid:12) f ′ i − ˜ f i (cid:12)(cid:12)(cid:12) ≤ (cid:18) K + K
16 + 1 (cid:19) L ˆ h + µ i − · T · h p i = O (ˆ h ) + µ i − · T · h p i . (32)Now, from 2 i − ( i − ≤ ˆ h ↔ ( i − ( i − ≤ log(ˆ h ) ↔ i ≤ ( i −
2) + log(ˆ h )log(2) = ( i −
2) + log (ˆ h ) =: ˜ l Also, we impose that 1 < ˜ l = ( i −
2) + log (ˆ h ) → − log (ˆ h ) < i . Then, by Lemma 2.7 we have for 1 ≤ i ≤ ˜ l : | A − i i − | ≤
23 2 i − ( i − ≤ ˆ h and by Eqs. (31) and (32) if 1 ≤ i ≤ ˜ l | ˙ f i +1 − f ′ i +1 | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) i − X j =1 A − i,j r j +1 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ O (ˆ h ) + | A − i i − | · T · ˆ h p i ≤ O (ˆ h ) + ˆ h · O (ˆ h p i ) = O (ˆ h min (3 ,p i +1) ) . Thus, we define l = ˜ l + 1 and obtain: | ˙ f i − f ′ i | = (cid:26) O (ˆ h min(3 ,p i +1) ) , ≤ i ≤ l ; O (ˆ h p i ) , l < i ≤ i − . (33) (cid:4) roposition 2.9. Let us assume that ˆ h < , f ( x ) ∈ C ([ x i +1 , x n ]) and let L > be such that | f (4) ( x ) | ≤ L ,for all x ∈ [ x i +1 , x n ] . Let p i , T > and A i +0 , ˙ F i +0 , B i +0 be as defined in Eqs. (25) , (26) and (27) , which satisfy | f ′ i − ˜ f i | = T · ˆ h p i and A i +0 ˙ F i +0 = B i +0 . If n > i + 1 − log (ˆ h ) and there exist K > , such that ˆ h/h i ≤ K for all i = i + 1 , . . . , n then:There exists an integer i + 1 < l , such that: | ˙ f i − f ′ i | = (cid:26) O (ˆ h min(3 ,p i +1) ) , l ≤ i ≤ n − O (ˆ h p i ) , i + 1 ≤ i < l . (34)The proof is similar to Prop. 2.8, the key is again to use Lemma 2.7.The following corollary summarizes the order obtained in the approximation of the derivatives if the originalsystem (15) is changed by (28). The order ot the approximation is reduced in a neighborhood of a point x i wheremonotonicity constraints do not hold when the original system is used, and, in exchange, C regularity is maintained,except at the point x i itself. the corollary is a direct consequence of Propositions. 2.8 and 2.9. Corollary 2.10.
Let us assume that ˆ h < , f ( x ) ∈ C ([ x , x n ]) and let L > be such that | f (4) ( x ) | ≤ L , for all x ∈ [ x , x n ] . Let p i , T > and A i , ˙ F i , B i be as defined in Eqs. (25) , (26) and (27) which satisfy | f ′ i − ˜ f i | = T · ˆ h p i and A i ˙ F i = B i . If − log (ˆ h ) < i < ( n − (ˆ h ) and there exists K > , such that ˆ h/h i ≤ K for all i = 1 , . . . , n then:There exist integers l < i < l , such that: | ˙ f i − f ′ i | = O (ˆ h min(3 ,p i +1) ) , ≤ i ≤ l ; O (ˆ h p i ) , l < i < l ; O (ˆ h min(3 ,p i +1) ) , l ≤ i ≤ n − . (35) Moreover the cubic spline interpolator defined in Eq. (2) using F i as an approximation of the values of the firstderivatives, has C regularity except at point x i . In the case of a piecewise C function that has a smaller smoothness at an interval [ x i , x i +1 ], the following resultholds [16]: Proposition 2.11.
Let us assume that ˆ h < , f ( x ) ∈ C ([ x , x i ] ∪ [ x i +1 , x n ]) and let L > be such that | f (4) ( x ) | ≤ L , for all x ∈ [ x , x i ] ∪ [ x i +1 , x n ] . Let r > , p i , p i +1 ≥ and A i − , ˙ F i − , B i − , A ( i +1) + , ˙ F ( i +1) + , B ( i +1) + be as defined in Eqs. (25) , (26) and (27) which satisfy | f ′ l − ˜ f l | = T · ˆ h p l , l = i , i + 1 and A i ˙ F i = B i . If − r log (ˆ h ) < i < ( n −
1) + r log (ˆ h ) and there exists K > such that ˆ h/h i ≤ K for all i = 1 , . . . , n then there existintegers ˜ l < i , i + 1 < ˜ l , such that: | ˙ f i − f ′ i | = O (ˆ h min(3 ,r + p i ) ) , ≤ i ≤ ˜ l ; O (ˆ h min( p i ,p i ) , ˜ l < i < ˜ l ; O (ˆ h min(3 ,r + p i ) ) , ˜ l ≤ i ≤ n − , (36) being: ˜ l = ( i −
1) + r log (ˆ h ) , ˜ l = ( i + 2) − r log (ˆ h ) .
3. Non-linear computation of derivatives
In Sections 2.2 and 2.3 we have shown two ways to construct monotonicity-preserving cubic splines by replacingsome derivative values whenever necessary. In both cases, we have shown the importance of the order of accuracyof the approximate derivative values used as a replacement. In this section we review different ways to design thesevalues and we study the respective order obtained by the interpolant.Let us start by a formula designed by Fritsch and Butland [10].˙ f F Bi := m i − m i m i − +2 m i , if | m i | ≤ | m i − | ; m i − m i m i +2 m i − , if | m i | > | m i − | ;0 , if m i − m i ≤
0. (37)9n this case the value ˙ f F Bi is defined in a way that automatically satisfies (5). If f is smooth and m i − m i > f F Bi = f ′ ( x i ) + O (ˆ h ) . (38)By Lemma 1.5, the cubic Hermite interpolant (2) is at least second-order accurate (see [3]).A second possibility is Brodlie’s formula [10]:˙ f Bi := ( wl i + wr i ) m i − m i wl i m i + wr i m i − , m i m i − ≥ m i m i − <
0; (39)where wl i = h i − + 2 h i , wr i = 2 h i − + h i . Formula (39) is implemented in the PCHIP program of Matlab ([17]).Using the properties of the weighted harmonic mean the following results are proved in [3]: Lemma 3.1.
Let us assume that f is smooth and m i m i − > . Then, If h i − = h i then ˙ f Bi = f ′ ( x i ) + (cid:26) O (ˆ h ) if h i − = h i = ˆ h,O (ˆ h ) if h i − = h i . (40) Theorem 3.2.
Let us assume that m i m i − > . If the derivatives ˙ f i are computed using Brodlie’s formula (39),then the cubic Hermite interpolant P i ( x ) constructed from (2) is monotone. Moreover, if f is smooth then P i ( x ) isthird-order accurate when h i − = h i ∀ i and second-order accurate otherwise. Brodlie’s formula produces a third order interpolant only in the case of using equally-spaced grids. In [5] (see also[3]) Ar`andiga and Y´a˜nez introduce a new method to compute the approximated derivatives based on the weightedharmonic mean that achieve third order of accuracy for non-uniform grids and preserves monotonicity. The proposedformula is: ˙ f AYi := sign ( m i ) ( h i + h i − ) /p i | m i − | | m i | ( h i − | m i − | p i + h i | m i | p i ) /p i if m i m i − ≥ m i m i − < , (41)where p i = max(1 , log( w i )log(3) ) and w i = 2 max( h i − , h i ) / min( h i − , h i ). With this formula, the following propositionholds: Proposition 3.3. ([5]) Let us assume that m i m i − > and there exists K > such that ˆ h/h j < K , for all j =1 , . . . , n . If we obtain the derivatives using Ar`andiga-Y´a˜nez’s formula, (41), then the cubic Hermite interpolant (2) ismonotone. Moreover, if f is smooth then P i ( x ) is third-order accurate. Note that in the case of equally-spaced grids the formulas by Ar`andiga and Y´a˜nez and Brodlie coincide.
4. Numerical experiments
In this section we present some experiments to verify the theoretical results previously obtained. In particular,we will divide our experiments in two subsections: in 4.1 we study the order of approximation of the different re-constructions using smooth or piecewise smooth functions. We perform two experiments with both equally- and notequally-spaced grids. In the first case methods B and AY are the same.On the other hand, in 4.2 we check the monotonicity property in cases where the function is unknown and onlynodal values are given.In this section each method will be identified by an acronym, being: S: Cubic spline with boundary conditions S ′ ( a ) = f ′ ( a ), S ′ ( b ) = f ′ ( b ). If we do not know these boundary conditionswe impose S ′ ( a ) = m , S ′ ( b ) = m n − . O k : Method explained in Section 2.2. The approximations of the first derivative are computed using system (14)–(15) and those values which do not satisfy the conditions by Theorem 1.4 are replaced by new approximationsobtained through the methods explained in Section 3. R k : Method explained in Section 2.3. It is a cubic spline with boundary conditions S ′ ( a ) = f ′ ( a ), S ′ ( b ) = f ′ ( b ) butconstructed using system (14)–(15) instead of (28) because an approximate derivative value is modified.For methods O and R we introduce the subscript k to indicate the approximation to the derivative used, thus k = F B, B, or AY , Eqs. (37), (39) and (41) respectively.10 .1. Accuracy We divide this section in two parts: Firstly, we analyze the case of equally-spaced grids. We will check that theorder of accuracy of the approximation to the derivatives’ values is four at smooth parts. Secondly, we perform someexperiments using a non-uniform grid to discretize the functions. In both cases, we explore the order of approximationat the points depending on the distance to the discontinuity.
In this section we consider two experiments: in the first one the function is smooth and we replace the approximationof the derivative at a single point to check the effect of this new value in the smoothness and accuracy of the spline;in the second one we consider a piecewise smooth function with a jump discontinuity.Experiment 1. In order to check the order of approximation of the methods we consider the following smoothfunction: f ( x ) = x + sin( x ) , (42)and discretize it on [0 ,
2] using a uniform grid: x lj = j/ l , j = 0 , . . . , l +1 , being l a fixed positive integer. We establisha window W ⊆ { , . . . , l +1 } , that selects a subset of the points in the discretization. Errors and numerical orders ofthe various methods are computed in the selected points, in order to verify the properties stated Sections 2.2, 2.3 and3. The errors are computed in the window using: e Wl = max i ∈ W | f ′ ( x i ) − ˙ f i | and the order of accuracy of the approximation are estimated by computing o W = log e Wl e Wl − ! . For methods O and R we replace the derivative value ˙ f l corresponding to the point i = 2 l by new values computedby the methods in Section 3, so as to verify the accuracy and smoothness properties stated in Sections 2.2 and 2.3.We first consider the window W = { i : 0 ≤ i ≤ l +1 } With this setup the order of accuracy is determined bythe approximation of the first derivative made in the point i . As shown in Table 1, in the case of B and AY , weobtain second order in accordance with Eq. (40); for F B method, it is reduced by Eq. (38). Finally we remark that,according to (21) the order is four for the S algorithm as the grid is uniform. h S R
F B O F B R B = R AY O B = O AY . e − . . . . . . e − . . . . . . e − . . . . . . e − . . . . . Table 1: Experiment 1 with h = 2 − l (equally-spaced grid) and estimated orders log ( e W l /e W l − ), 5 ≤ l ≤ W = { i : 2 ≤ i ≤ l +1 − } . If the window considered for order estimation is reduced so as to exclude i , that is, W = W \ { i } , the order forthe O methods is increased up to four, in agreement with Prop. 2.6. In contrast, for the R methods the order doesnot increase as the order reduction affects points close to i , according to Cor. 2.10 (see Table 2). h S R
F B O F B R B = R AY O B = O AY . e − . . . . . . e − . . . . . . e − . . . . . . e − . . . . . Table 2: Experiment 1 with h = 2 − l (equally-spaced grid) and estimated orders log ( e W l /e W l − ), 5 ≤ l ≤ W = W \ { i } . l = ( i −
1) + log (ˆ h ) = (2 l −
1) + log (2 − l ) = 2 l − l − ,l = ( i + 1) − log (ˆ h ) = (2 l + 1) − log (2 − l ) = 2 l + l + 1 , and define W = { i : 1 ≤ i ≤ l or l ≤ i ≤ l +1 − } . As we can see in Table 3, the order increases in R methods fromtwo up to four in B = AY methods and from one to three for R F B . In this case, the size of the chosen window aroundthe discontinuity is sufficiently large as to increase the order of accuracy at the rest of points from the order obtainedfor all the points showed in Table 2. It would be possible to reduce this interval if we take the bound indicated in Eq.(29) for equally-spaced grids. In this way, the constructed spline has maximum order and regularity at the interval[ x , x l ] ∪ [ x l , x n ] for all methods, according to Lemma 1.5. h S R
F B O F B R B = R AY O B = O AY . e − . . . . . . e − . . . . . . e − . . . . . . e − . . . . . Table 3: Experiment 1 with h = 2 − l (equally-spaced grid) and estimated orders log ( e W l /e W l − ), 5 ≤ l ≤ W = { i : 1 ≤ i ≤ l or l ≤ i ≤ l +1 − } . Experiment 2. In this experiment we consider a piecewise smooth function with a jump discontinuity located in theinterval [ x i , x i +1 ]. The presence of the discontinuity produces two effects: first, the approximation of the derivativesobtained from (14) will suffer the Gibbs phenomenon and produce some spurious oscillations near the discontinuity andhence monotonicity will not be preserved; second, the methods discussed in Section 3 will not attain their maximumaccuracy.We consider the function: g ( x ) = (cid:26) x + sin( x ) , ≤ x ≤ , x + cos( x ) , < x ≤ , (43)and discretize it in [0 ,
2] analogously to previous subsection.This function has a jump discontinuity at x = 1 thatcorresponds to the node x l . In Table 4 we display the numerical orders obtained in the derivative approximation when we take ˜ W = W \{ l + l + 1 } with W being as defined in the previous experiment, and in the left column of Fig. 1 the reconstructionsproduced by the different methods are shown. The spline S produces oscillations due to the violation of monotonicityconstraints and produces a poor order of accuracy in the reconstruction.The same results are obtained for the O methods because the aproximation of the first derivative is no modified inthe nodes belonging to ˜ W . Finally there is an improvement in the order of accuracy if the R methods are applied. Inthe right plots of Fig. 1 the errors obtained in the derivative computation are shown. It can be seen that all methodsproduce big errors around the discontinuity, being the ones corresponding to the spline reconstruction one order ofmagnitude bigger than the rest. On the other hand, methods B and AY produce errors that are smaller by a factorof around 1 / F B . Also, we observe that, as expected, the O methods suffer an accuracyloss in a bigger neighbourhood around [ x i , x i +1 ] than the R methods. h S R
F B O F B R B = R AY O B = O AY . e − . . . . . . e − . . . . . . e − . . . . . . e − . . . . . Table 4: Experiment 2 with h = 2 − l and estimated orders log ( e ˜ W l /e ˜ W l − ), 5 ≤ l ≤
8, ˜ W = { i : 1 ≤ i ≤ l or l ≤ i ≤ l +1 − } , whenthe grid is uniform. S R F B O F B R B = R A Y O B = O A Y Figure 1: Experiment 2. l = 4 (a) Reconstructions using the different methods, (b) | g ′ ( x i ) − ˙ g ki | , i = 0 , . . . ,
13s expected, the maximum order is not obtained. According to Prop. 2.11, we consider a window W = { i : 1 ≤ i ≤ ˜ l or ˜ l ≤ i ≤ l +1 − } for˜ l = ( i −
1) + r log (ˆ h ) = (2 l −
1) + r log (2 − l ) = 2 l − r · l − , ˜ l = ( i + 2) − r log (ˆ h ) = (2 l + 2) − r log (2 − l ) = 2 l + r · l + 2 , and r = 2. The results corresponding to this window are shown in Table 5. The accuracy orders for O and S coincide,whilst for R the order is determined by the method used in the derivative computation. h S R
F B O F B R B = R AY O B = O AY . e − . . . . . . e − . . . . . . e − . . . . . . e − . . . . . Table 5: Experiment 2 with h = 2 − l (equally-spaced grid) and estimated orders log ( e W l /e W l − ), 5 ≤ l ≤ W = { i : 1 ≤ i ≤ ˜ l or ˜ l ≤ i ≤ l +1 − } . In order to check the theoretical results in non-uniform grids, we consider the same functions as in subsection 4.1.2but using a different discretization.Experiment 1 with a non-uniform grid In this experiment we discretize the function f ( x ) in (42) at the interval[0 ,
2] using a non-uniform grid constructed according to the following procedure: Let l be a fixed positive integer, wedefine the points as: (cid:26) x l i = 2 − l · i,x l i +1 = 2 − l ( i + ) , (44)with i = 0 , . . . , l +1 − x l l +1 = 2. It is clear that ˆ h = − l . As in the case of uniform grids, for R and O methodsthe value of the derivative computed by 14 at the point ˙ f l is replaced by a new value computed using methods F B , B and AY . In this case the node immediately at the right of the discontinuity, x l l +1 , does not coincide with any nodein the grid corresponding to l + 1, but with the node x l +22 l +2 +2 , belonging to the grid corresponding to l + 2. Indeed: x l l +1 = 2 − l (cid:18) l − + 14 (cid:19) = 2 − l (cid:18) l +1 + 14 (cid:19) = 2 − ( l +2) (cid:0) l +1 + 1 (cid:1) = x l +22(2 l +1 +1) = x l +22 l +2 +2 . Taking it into account, in order to estimate the order of accuracy of the approximation we use the following formula:˜ o W = log e Wl e Wl − ! . We take the following windows, similar to the ones defined in Experiment 1 for uniform grids, i.e.: W n = { i : 0 ≤ i ≤ l +2 } ,W n = W n \ { i } ,W n = { i : 0 ≤ i ≤ l n or l n ≤ i ≤ l +2 } , being l n = ( i −
1) + log (ˆ h ) = (2 l +1 −
1) + log (cid:18)
34 2 − l (cid:19) ,l n = ( i + 1) − log (ˆ h ) = (2 l +1 + 1) − log (cid:18)
34 2 − l (cid:19) . (45)It can be seen in Tables 6, 7 and 8 that the orders of approximation obtained for R B , O B and R AY , O AY are differentand are in accordance with Lemma 3.2 and Prop. 3.3. For the S method the expected order, O (ˆ h ) is obtained.14 h S R
F B O F B R B R AY O B O AY . e −
02 2 . . . . . . . . e −
02 2 . . . . . . . . e −
03 3 . . . . . . . . e −
03 2 . . . . . . . Table 6: Experiment 1 with ˆ h = − l (non-equally-spaced grid) and estimated orders log ( e W n l /e W n l − ), 3 ≤ l ≤ W n = { i : 0 ≤ i ≤ l +2 } . If the window is changed by W n so as to avoid the point where the derivative was replaced, the order increases for O methods, in agreement with Prop. 2.6, and these methods achieve the maximum order as they coincide with the S method in the points of the window, see Table 7. However, the order is not improved using R methods because thesize of the window is not sufficiently large.ˆ h S R
F B O F B R B R AY O B O AY . e −
02 2 . . . . . . . . e −
02 2 . . . . . . . . e −
03 3 . . . . . . . . e −
03 2 . . . . . . . Table 7: Experiment 1 with ˆ h = − l (non-equally-spaced grid) and estimated orders log ( e W n l /e W n l − ), 3 ≤ l ≤ W n = W n \ { i } . When the window considered for the estimation of the order is W n , i.e., some more points around i are excludedfrom the order estimation according to Props. 2.8 and 2.9, the order of accuracy obtained in the points in the windowis optimal using any method. Even when the F B method is used the order obtained is 3, instead of the expectedsecond order (Cor. 2.10)). The reason can be either the regularity of the function used in this example or the factthat the values l n and l n taking from Eq. (45) are too large.ˆ h S R
F B O F B R B R AY O B O AY . e −
02 2 . . . . . . . . e −
02 2 . . . . . . . . e −
03 3 . . . . . . . . e −
03 2 . . . . . . . Table 8: Experiment 1 with ˆ h = − l (non-equally-spaced grid) and estimated orders log ( e W n l /e W n l − ), 3 ≤ l ≤ W n = { i : 0 ≤ i ≤ l n or l n ≤ i ≤ l +2 } . In order to analyze a function with a strong gradient, in this subsection, we discretize g ( x ), Eq. (43), at the interval[0 ,
2] using the non-uniform grid defined by Eq. (44). We take the following window:˜ W n = W n \ { l n } . In Table 9 we can see that the order of the O methods is one, due to the presence of the discontinuity, as thehypothesis of Prop. 3.3 are not satisfied. For R methods we obtain an improvement in the order but not the optimal.In order to get it we take the following window: W n = { i : 0 ≤ i ≤ ˜ l n or ˜ l n ≤ i ≤ l +2 } , with ˜ l n = ( i −
1) + 2 log (ˆ h ) = (2 l +1 −
1) + 2 log (cid:18)
34 2 − l (cid:19) , ˜ l n = ( i + 1) − (ˆ h ) = (2 l +1 + 1) − (cid:18)
34 2 − l (cid:19) . h S R
F B O F B R B R AY O B O AY . e −
02 0 . . . . . . . . e −
02 1 . . . . . . . . e −
03 1 . . . . . . . . e −
03 1 . . . . . . . Table 9: Experiment 2 with ˆ h = − l (non-equally-spaced grid) and estimated orders log ( e ˜ W n l /e ˜ W n l − ), 5 ≤ l ≤
8, ˜ W n = W n \ { l n } . The results corresponding to this setup are shown in Table 10. Third order of accuracy is obtained for all methods,in agreement with the theoretical results.As a conclusion, the kind of reconstructions proposed in the paper allow for replacing the approximations of thederivatives in some points, in order to ensure monotonicity preservation while maintaining optimal order at pointsthat are located at a certain distance from them.ˆ h S R
F B O F B R B R AY O B O AY . e −
02 3 . . . . . . . . e −
02 3 . . . . . . . . e −
03 3 . . . . . . . . e −
04 3 . . . . . . . Table 10: Experiment 2 with ˆ h = − l (non-equally-spaced grid) and estimated orders log ( e W n l /e W n l − ), 4 ≤ l ≤ W n = { i : 0 ≤ i ≤ ˜ l n or ˜ l n ≤ i ≤ l +2 } . In this section we show an example, indicated as experiment 3 (Table 11), in which we compare the reconstructionsobtained with the different methods considered in this paper. These data have been used in [11]. In this example weonly know the values of the function at certain nodes. The discretization is not equally spaced and, therefore, methods B and AY are different. i x i .
99 8 .
09 8 .
19 8 . . y i . e − . e − . . . . . . Table 11: Data for the experiment 3
The algorithm modifies the values at the nodes i = 2 , , ,
8. In this case, we obtain monotone reconstruction inthe complete interval (Figure 2). The results are very similar in all cases, i.e., both methods produce similar monotonecurves because the values which are modified are the most relevant in this example.
5. Conclusions
In this work, we have introduced two new algorithms to obtain monotone cubic spline interpolants. We haveconsidered the case in which there exists a discontinuity or a high gradient in the data, and considered two options,based on modifying the approximation to the derivatives in the points where monotonicity constraints are violated.In the first one we rewrite the spline system fixing the modified derivatives and recompute the derivatives using themodified spline system at both sides. Using this algorithm the regularity is C in all points except at the ones wherethe derivative approximation was modified, and the order is reduced in a neighborhood of these, but it is conserved inthe rest of the interval. In the second algorithm we again modify the approximations of the derivatives where requiredbut the rest of the values are kept as initially computed. In this case we conserve the order in all points but theregularity is lost in a neighborhood of the points where the derivative was modified. Some numerical tests confirmthese results. 16
10 12 14 16 18 2000.20.40.60.811.2 8 10 12 14 16 18 2000.20.40.60.811.2 8 10 12 14 16 18 2000.20.40.60.811.2 8 10 12 14 16 18 2000.20.40.60.811.2
S R
F B R B R AY Figure 2: Reconstructions obtained for Experiment 3 using the different methods
References [1]
A. A. Abushama and B. Bialecki (2008): “Modified nodal cubic spline collocation for poisson’s equation.”, SIAM J.Numer. Anal., 46(1), 397–418.[2]
H. Akima (1970) “A new method of interpolation and smooth curve fitting based on local procedures”, J. Assoc. Comput.Mach., 17, 589–602.[3]
F. Ar`andiga (2013): “On the order of nonuniform monotone cubic Hermite interpolation”, SIAM J. Numer. Anal.,51(5), 2613–2633.[4]
F. Ar`andiga, A. Baeza and D. F. Y´a˜nez (2013): “A new class of non-linear monotone Hermite interpolants”, Adv.Comput. Math., 39, 289–309.[5]
F. Ar`andiga and D. F. Y´a˜nez (2019): “Third-order accurate monotone cubic Hermite interpolants ”, Appl. Math.Letters, 94, 73–79.[6]
A. M. Bica (2012): “Fitting data using optimal Hermite type cubic interpolating splines”, Appl. Math. Letters, 25,2047-2051.[7]
C. de Boor : “A practical guide to splines”, Springer-Verlag, 2001.[8]
R. J. Cripps and M. Z. Hussain (2012): “C monotone cubic Hermite interpolant”, Appl. Math. Letters, 25, 1161-1165.[9] M. S. Floater and M.-J. Lai (2016): “Polygonal spline spaces and the numerical solution of the Poisson equation”,SIAM J. Numer. Anal., 54(2), 797–824.[10]
F. N. Fritsch and J. Butland (1984): “A method for constructing local monotone piecewise cubic interpolants”,SIAM J. Sci. Stat. Comput., 5, 2, 300-304.[11]
F. N. Fritsch and R. E. Carlson (1980): “Monotone piecewise cubic interpolation”, SIAM J. Numer. Anal., 17, 2,238-246.[12]
J. H. Hyman (1983): “Accurate monotonicity preserving cubic interpolation”, SIAM J. Numer. Anal., 4, 4, 645-654.[13]
H. T. Huynh (1993): “Accurate monotone cubic interpolation”, SIAM J. Numer. Anal., 30, 1, 57-100.[14]
D. Kershaw, (1970): “Inequalities on the elements of the inverse of a certain tridiagonal matrix”, Math. Comp., 24,155-158.[15]
D. Kershaw, (1971): “A note on the convergence of interpolatory cubic splines”, SIAM J. Numer. Anal., 8, 67-74.[16]
D. Kershaw, (1972): “The orders of approximation of the first derivative of cubic splines at the knots”, Math. Comp.,26, 191-198.[17]
C. Moller : “Numerical Computing with MATLAB”, SIAM, Philadelphia, 2004.[18]
J. Stoer and R. Bulirsch : “Introduction to Numerical Analysis”, Springer-Verlag, 1980.[19]
G. Wolberg and I. Alfy (2002): “An energy-minimization framework for monotonic cubic spline interpolation”, J.Comput. and Applied Math., 143, 145-188.(2002): “An energy-minimization framework for monotonic cubic spline interpolation”, J.Comput. and Applied Math., 143, 145-188.