C2 Cubic Algebraic Hyperbolic Spline Interpolating Scheme by Means of Integral Values

  Citation: Eddargani, S.; Oraiche, M.; Lamnii, A.; Louzar, M. C2Cubic Algebraic Hyperbolic Spline Interpolating Scheme by Means of Integral Values. Mathematics 2022,10, 1490. https://doi.org/10.3390/ math10091490 Academic Editor: Ioannis K. Argyros Received: 21 March 2022 Accepted: 27 April 2022 Published: 29 April 2022 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). mathematics Article C2Cubic Algebraic Hyperbolic Spline Interpolating Scheme by Means of Integral Values Salah Eddargani 1,* , Mohammed Oraiche 2, Abdellah Lamnii 3and Mohamed Louzar 2 1Department of Applied Mathematics, University of Granada, 18071 Granada, Spain 2 MISI Laboratory, Faculty of Sciences and Techniques, Hassan First University of Settat, Settat 26000, Morocco; [email protected] (M.O.); [email protected] (M.L.) 3LaSAD Laboratory, Ecole Normale Supérieure, Abdelmalek Essaadi University, Tetouan 93030, Morocco; [email protected] *Correspondence: [email protected] Abstract: In this paper, a cubic Hermite spline interpolating scheme reproducing both linear polynomials and hyperbolic functions is considered. The interpolating scheme is mainly defined by means of integral values over the subintervals of a partition of the function to be approximated, rather than the function and its first derivative values. The scheme provided is C2 everywhere and yields optimal order. We provide some numerical tests to illustrate the good performance of the novel approximation scheme. Keywords: algebraic hyperbolic splines; integro cubic interpolation; Hermite representation MSC: 65D07 1. Introduction Nowadays, numerical methods are a common tool, just a click away from the user. Interpolation is a particular and very important numerical method, which is widely used to address the solution of theoretical problems and show their full potential to numerically solve problems that occur in many different branches of science, engineering and economics. The interpolation approximants should be easily evaluated, differentiated and integrated. Spline functions, i.e., smooth piecewise polynomial functions, fulfil all these requirements. Since the introduction of the systematic study of spline functions by I. J. Schoenberg in the 1940s [ 1 ], they have become an indispensable tool in approximation theory and numerical computation, including computer-aided geometric design (CAGD) [ 2 , 3 ], the numerical solution of PDEs, numerical quadratures, interpolation and quasi-interpolation, regularization, least squares, isogeometric analysis and image processing. Spline functions have been the subject of many research results that have been presented in well-known books [ 4 , 5 ]. Furthermore, thousands of papers related to spline functions and their applications have been published in the last five years, in fields ranging from computer science, engineering, physics and astronomy to entertainment and conceptual design assistance. Polynomial spline functions are the most commonly used class, especially due to the fact that they admit normalized bases on any bounded interval [a , b] (Bernstein bases and B-spline bases); for more details, see [ 5 – 8 ]. Indeed, Bernstein bases and B-splines possess several interesting properties, such as non-negativity, local support, partition of unity and totally positive bases [ 9 ]. Bernstein basis functions of degree n are the best among all bases of the polynomial space of degree less than or equal to n . This means that it is the basis relative to which the control polygon of any curve yields the best information on the curve itself. The non-polynomial B-splines have also been studied in the literature. For instance, the trigonometric B-splines were presented in [ 10 , 11 ]. The authors in [ 12 ] Mathematics 2022,10, 1490. https://doi.org/10.3390/math10091490 https://www.mdpi.com/journal/mathematics Mathematics 2022,10, 1490 2 of 13 have established a recurrence relation for the trigonometric B-splines of arbitrary order. A complete construction of exponential tension B-splines of arbitrary order was given in [ 13 ]. A more recent study on this kind of spline was given in [ 14 , 15 ]. The splines associated with the exponential B-spline space are often referred to as hyperbolic splines [16]. Unfortunately, neither Bernstein bases nor B-splines are suitable to perfectly describe conic sections, which are shapes of major interest in certain engineering applications. This resulted in the introduction of NURBS [ 17 ], which can be seen as a generalization of B-splines, inheriting from them important properties and with the additional benefit of making possible the exact representation of conic sections. On the other hand, the NURBS representation suffers from some drawbacks that are considered critical in CAD. In fact, the necessity of weights does not have an evident geometric meaning and their selection is often unclear. Furthermore, it behaves awkwardly with respect to differentiation and integration, which are indispensable operators in analysis. On this concern, it is sufficient to think about the complex structure of the derivative of a NURBS curve of a given order. An alternative is to use the so-called generalized B-splines; see [ 18 , 19 ] and references therein. The generalized B-splines belong to the extended space spanned by 1, x, . . . , xn−2,u1(x),u2(x) , where u1 and u2 are smooth functions. The two functions u1 and u2 can be selected to achieve the exact representation of salient profiles of interest and/or to obtain particular features. The most popular choices of these functions are: (u1(x),u2(x))=(sin(x), cos(x)) and (u1(x),u2(x))=(exp(x), exp(−x)) , which yield algebraic trigonometric and algebraic exponential splines, respectively. The algebraic exponential splines are often referred as algebraic hyperbolic splines. Algebraic trigonometric and hyperbolic splines allow an exact representation of conic sections, as well as of some transcendental curves, such as helix and cycloid curves. In fact, they are in a position to provide parametrizations of conic sections that are significantly more related to the arc length than NURBS. These classes of splines are also known as cycloidal spaces, and they have become the subjects of a considerable amount of research [ 2 , 20 – 26 ]. The algebraic hyperbolic spaces spanned by the functions 1, x , . . . , xn−2 , cosh(x) , sinh(x) , for x∈R , have been widely considered in the literature; see [ 21 ] and references quoted. They yield the tension splines, which are extremely useful for avoiding undesirable oscillations in the interpolation curves [27]. In this paper, we consider the algebraic hyperbolic (AH) cubic spline space, spanned by {1, x, sinh, cosh} , and let Xn={xi:=a+ih}n i=0 be a uniform partition of a bounded interval I=[a,b] , with h=b−a n . Given values fj i , i= 0, . . . , n , j= 0, 1, there exists a unique cubic AH Hermite interpolant si∈span{1, x, sinh, cosh} such that s(j) i(xi+k) = fj i+k , i= 0, . . . , n− 1, j , k= 0, 1. The spline s defined from the local interpolants si is a C1 continuous AH spline that interpolates the data fj i. In this work, we suppose that the data fj i , i= 0, . . . , n , j= 0, 1, are not given, and we assume that the integral values over subintervals [xi , xi+1] , i= 0, . . . , n− 1 are given. Then, our purpose is the construction of the Hermite interpolant s using only this information. This kind of approximation arises in various fields, such as mechanics, mathematical statistics, electricity, environmental science, climatology, oceanography and so on (for more details, see References [ 28 , 29 ] and references therein). More precisely, the values fj i will be computed by means of C2 smoothness conditions at the knots xi and the integral values over subintervals [xi , xi+1] . This is done by solving a three-diagonal linear system. Some final conditions are required. In particular, we assume that the three values f0 0 , f0 n and f1 0 or f1 n are given. In general, these three values are not always available. We suggest a modified scheme that does not involve any final conditions to avoid this limitation. Integro spline approximation was treated in various works in the literature. The author in [ 30 , 31 ] developed two types of integro spline approximants, cubic and quintic cases, respectively. The two schemes introduced in [ 30 , 31 ] require various end conditions and the solution of a three-diagonal system of linear equations. Solving a linear system Mathematics 2022,10, 1490 3 of 13 of equations sometimes is very expensive, so the authors in [ 32 ] developed cubic integro splines quasi-interpolant without solving any system of equations. An integro quartic spline scheme has been constructed in [ 33 ]. The authors in [ 21 , 34 ] provided some integro spline schemes for the case of non-polynomial splines. More recent work on the integro spline approximation is given in [35,36] In this work, a new class of integro spline approximant is introduced. The proposed operator is C2 smooth everywhere and exactly reproduces both linear polynomials and hyperbolic functions, which is useful to avoid undesirable oscillations in curves’ interpolation. Some end conditions are needed, and to avoid this inconvenience, we have proposed a modified scheme that does not require additional end conditions. This paper is organized as follows: in Section 2, we first present the C1 cubic AH interpolating scheme, which is defined by means of the value and the derivative value at each knot of the partition. We also study the error bound of the presented scheme. Then, the conditions for achieving the C2 smoothness are described. Thus, an approach to define the integro spline scheme by means of integral values is proposed. In Section 3, we provide some numerical tests. Finally, we present some conclusions. 2. Algebraic Hyperbolic Spline Interpolation Let Xn={xi:=a+ih}n i=0 be a uniform partition of a bounded interval I=[a,b] , with h=b−a n. 2.1. Cubic Algebraic Hyperbolic (AH) Splines of Class C1 The construction of a C1 cubic AH spline interpolant on the partition Xn should be locally expressed in each sub-interval `i:=[xi,xi+1] in terms of the function and the first derivative values of the approximated function at knots xiand xi+1. The space of cubic AH splines on Xnwith global C1continuity is denoted as S1 3(Xn):=ns∈C1(I);s|`i∈Γ3for all i=0, . . . , no, where Γ3:=span{1, x, sinh, cosh} stands for the linear space of cubic AH splines. Its dimension equals 2 (n+ 1 ) . The following Hermite interpolation problem can then be considered: there exists a unique spline s(x)∈S1 3(Xn)such that s(j)(xi) = fj i,i=0, . . . , n,j=0, 1, (1) for any given set of fj i− values. Therefore, each spline s(x)∈S1 3(Xn) restricted to the sub-interval `ican be represented as follows: s|`i(x):= 1 ∑ k=0 1 ∑ j=0 fj i+kφj i+k(x), (2) in which φj i+k , k= 0, 1, j= 0, 1 are classical Hermite basis functions of S1 3(Xn) restricted to the sub-interval `i . The basis functions φj i+k , k= 0, 1, j= 0, 1 are the unique solution of the interpolation problem given by (1) in s(x)∈S1 3(Xn) , where f0 i,f1 i,f0 i+1,f1 i+1=(1, 0, 0, 0) , f0 i,f1 i,f0 i+1,f1 i+1=(0, 1, 0, 0) , f0 i,f1 i,f0 i+1,f1 i+1=(0, 0, 1, 0) and f0 i,f1 i,f0 i+1,f1 i+1= (0, 0, 0, 1), respectively. More precisely, Mathematics 2022,10, 1490 4 of 13 φ0 i(xi) = 1, φ0 i0(xi) = 0, φ0 i(xi+1) = 0, φ0 i0(xi+1) = 0, φ1 i(xi) = 0, φ1 i0(xi) = 1, φ1 i(xi+1) = 0, φ1 i0(xi+1) = 0, φ0 i+1(xi) = 0, φ0 i+10(xi) = 0, φ0 i+1(xi+1) = 1, φ0 i0(xi+1) = 0, φ1 i+1(xi) = 0, φ1 i+10(xi) = 0, φ1 i+1(xi+1) = 0, φ1 i+10(xi+1) = 1. They can be expressed as follows: φ0 i(x):=1 θ(1, x, sinh(x), cosh(x))    −cosh(h) + sinh(h)(h+xi)+1 −sinh(h) sinh(h+xi)−sinh(xi) cosh(xi)−cosh(h+xi)     , φ1 i(x):=1 θ(1, x, sinh(x), cosh(x))    −cosh(h)−sinh(h)xi+1 sinh(h) sinh(xi)−sinh(h+xi) cosh(h+xi)−cosh(xi)     , φ0 i+1(x):=1 θ(1, x, sinh(x), cosh(x))    hcosh(h)−sinh(h) + (cosh(h)−1)xi 1−cosh(h) cosh(xi)−cosh(h+xi)+hsinh(h+xi) −hcosh(h+xi)−sinh(xi)+sinh(h+xi)     , φ1 i+1(x):=1 θ(1, x, sinh(x), cosh(x))    −h+sinh(h) + (cosh(h)−1)xi 1−cosh(h) −cosh(xi)+cosh(h+xi)−hsinh(xi) hcosh(xi)+sinh(xi)−sinh(h+xi)     , where θ=hsinh(h)−2 cosh(h) + 2. Figure 1shows the typical plots of the basis functions φj i+k , k= 0, 1, j= 0, 1, on the sub-interval [xi,xi+1]=[0, 1]. ϕi+1 1.0 -0.2 0.2 0.4 0.6 0.8 1.0 Figure 1. Plots of the basis functions on [xi,xi+1]=[0, 1]. Consider that the values of fj i are related to an explicit function f , i.e., fj i=f(j)(xi) , j=0, 1. Then, it holds that ks−fk∞,[xi,xi+1]=Oh4. (3) Mathematics 2022,10, 1490 5 of 13 In order to prove this statement, consider an arbitrary but fixed value ˜ t different from xiand xi+1, and define R(t) = s(t)−f(t)−ρ 1 ∏ k=0 (t−xi+k)2, in which the constant ρis chosen such that R(˜ t) = 0, that is, ρ=s(˜ t)−f(˜ t) ∏1 k=0(˜ t−xi+k)2. The function Rhas at least three roots in [xi,xi+1], which are xi,xi+1and ˜ t. According to Rolle’s theorem, R0 has at least two roots in [xi , xi+1] that are different from xi , xi+1 , ˜ t and also R0(xi) = R0(xi+1) = 0, which means that R0 has at least four roots in [xi , xi+1] . Analogously and progressively, it is shown that R(2) has three roots in the interval [xi,xi+1],R(3)has two and R(4)has only one root, say ξ. It then states R(4)(ξ)=s(4)(ξ)−f(4)(ξ)−24ρ=0. It holds that s(˜ t)−f(˜ t) = 1 24s(4)(ξ)−f(4)(ξ)1 ∏ k=0 (˜ t−xi+k)2. This proves Equation (3). The proposed interpolant s is defined from the values and first derivative values at the break points. Unfortunately, this dataset is not always at hand. This paper deals with the case where neither the values nor the derivative values are known. Instead, we assume that the integral values over the sub-intervals are available. The strategy pursued in this work is the following: we first highlight the relationship between function and first derivative values by imposing the C2 smoothness at the set of break points. Next, we express the derivative values in the mean integral values provided. The C2 smoothness of s00(x) at xi , i= 1, ..., n− 1 yields the following consistency relations: αf1 i−1+βf1 i+αf1 i+1=f0 i+1−f0 i−1,i=1, ..., n−1. (4) where α=sinh(h)−h cosh(h)−1and β=csch2h 2(hcosh(h)−sinh(h)). The error related to the approximation scheme (4) can be derived from the following result. Lemma 1. Let f∈ C3([a,b]) ; then, the local truncation errors ti , i= 1, 2, . . . , n , associated with the scheme (4) are given by the expressions ti=−1 6h2f(3)(xi)(h(cosh(h) + 2)−3 sinh(h)) + Oh2. Mathematics 2022,10, 1490 6 of 13 Proof. The function fis supposed to be of class C3([a,b]), that is, f(xi+h)=f(xi)+f0(xi)h+f00(xi)h2 2+f(3)(xi)h3 6+Oh3, f(xi−h)=f(xi)−f0(xi)h+f00(xi)h2 2−f(3)(xi)h3 6+Oh3, f0(xi+h)=f0(xi)+f00(xi)h+f(3)(xi)h2 2+Oh2, f0(xi−h)=f0(xi)−f00(xi)h+f(3)(xi)h2 2+Oh2. Using Equation (4), it results that ti=αf0(xi−1)+βf0(xi)+αf0(xi+1)−1 csch2h 2(f(xi+1)−f(xi−1)). By replacing f0(xi−1) , f0(xi+1) and f(xi+1)−f(xi−1) by their Taylor expansions, the intended result can be achieved, which completes the proof. The next result can be easily deduced from the previous lemma. Theorem 2. Let f ∈ C3([a,b]), then f1 i−f0(xi)=Oh2. At this point, we have simply provided a scheme that approximates the derivative values from the function values by imposing C2 smoothness at the break points. In what follows, we will deal with the case where the function values themselves are not available, whereas the integrals over the sub-intervals are known. 2.2. Cubic Algebraic Hyperbolic Spline Interpolant Based on Mean Integral Value In traditional spline interpolation problems, it is assumed that the function values at the knots are given. In this subsection, the function values are supposed to be unknowns and we assume that the integrals over the sub-intervals [xi−1,xi] , i= 1, . . . , n are provided and are equal to ti=Zxi xi−1 f(x)dx,i=1, . . . , n(5) In the sequel, we will provide a scheme that approximates derivative values from the integrals ti,i=1, . . . , n. By integrating sover [xi−1,xi]and [xi,xi+1], one can obtain 2ti=hf0 i−1+f0 i+2−hcothh 2f1 i−f1 i−1,i=1, . . . , n(6) 2ti+1=hf0 i+f0 i+1+2−hcothh 2f1 i+1−f1 i,i=0, . . . , n−1 (7) respectively. Subtracting (6) from (7) as a first step, and then applying (4) as a second step, allows us to eliminate the unknowns f0 i , as well as to achieve new relations that link only the unknowns f1 iwith the provided data ti. It results that µf1 i−1+λf1 i+µf1 i+1=2(ti+1−ti),i=1, . . . , n−1, (8) with µ=1 24−h2csch2h 2and λ=h2−2cosh(h) + 2csch2h 2 Mathematics 2022,10, 1490 7 of 13 This yields a system of n− 1 linear equations, while there are n+ 1 unknowns f1 i . Then, two additional end conditions are required to determine the unknowns. The end conditions are the first derivative values at the end points a and b . Assume that f0(a) = f1 a and f0(b) = f1 bare provided. Then, a (n−1)×(n−1)linear tridiagonal system results.          λ µ µ λ µ µ λ µ ......... µ λ µ µ λ                 m1 . . . mn−1        =        b1−µc b2 . . . bn−2 bn−1−µd        , (9) with bi=2(ti+1−ti). The following result shows that the linear system (9) has a unique solution. Theorem 3. For h >0, the matrix system in (9) is a strictly diagonally dominant matrix. Proof. Let h> 0. It is easy to show that λ> 2 µ . Indeed, a simple calculation gives us the following equality: λ−2µ=h2−2cosh(h) + 2csch2h 2 hhcothh 2−2−24−h2csch2h 2 2hhcothh 2−2 =4 h+2 cothh 2 Since cothh 2>0, then λ−2µ>0, which completes the proof. The LU factorization is adequate for solving (9) because the index on the diagonally dominant property of the matrix A is equal to 1 5 for small enough h> 0. In fact, let A:=aij1≤i,j≤n be the matrix coefficient of system (9). The index on the diagonally dominant property of matrix Ais given by Dh(A):=max i=1,...,n 1 |aii| n ∑ j=1 j6=iai,j=2 max h>0 µ(h) λ(h)= 2 h−h 15 +13h3 6300 +Oh5 10 h+4h 15 −11h3 3150 +O(h5) Its limit when his close to zero is equal to 1 5. Once the values of f1 i , i= 1, . . . , n− 1 are determined, we can then compute f1 i , i= 1, . . . , n− 1 by means of (6). However, we still need another end condition. Suppose that one of the values f(a) and f(b) is provided; then, we can start from it and use (6) to obtain the remaining unknowns in an iterative way. To analyze the interpolation error, we need the following lemma to establish an error bound for our operator. Lemma 4. Let f∈ C3([a,b]) ; then, the local truncation errors ˜ ti , i= 1, 2, . . . , n , associated with the scheme (9), are given by the expressions ˜ ti=−1 6h2f(3)(xi)h2+3h2csch2h 2−12+Oh2. Theorem 5. Let f∈ C4([a,b]) and s be the interpolation operator defined by (2), (6) and (9). For a uniform step size h, we have kf(x)−s(x)k∞,[a,b]=O(h2). Mathematics 2022,10, 1490 8 of 13 Proof. Consider the sub-interval [xi , xi+1] . Let p be the spline in S1 3 , which satisfies p(j)(xi+k) = f(j)(xi+k),j,k=0, 1. Then, it results that ks−fk=ks−p+p−fk ≤ ks−pk+kp−fk ≤Oh2+Oh4 ≡Oh2. This due to the fact that ks−pk ≤ max|f0(xi)−f1 i|,|f0(xi+1)−f1 i+1| , which concludes the proof. In the general context, the end conditions may not be provided. Thus, to avoid this limitation, we will provide explicit expressions for the end conditions f(a) , f0(a) and f(b) by means of integral values. Lemma 6. For a given function f ∈C2, f(a)=11t1−7t2+2t3 6h,f0(a)=−2t1−3t2+t3 h2,f0(b)=−2tn−3tn−1+tn−2 h2. Proof. The Taylor expansion of faround ais as follows: f(x)=f(a)+f0(a)(x−a)+1 2f00(a)(x−a)2+O(x−a)2. t1=h f (a)+h2 2f0(a)+h3 6f00(a)+Oh3, t2=h f (a)+3h2 2f0(a)+7h3 6f00(a)+Oh3, t3=h f (a)+5h2 2f0(a)+19h3 6f00(a)+Oh3. Therefore, the following system is obtained:   t1 t2 t3  =   hh2 2h3 6 h3h2 27h3 6 h5h2 219h3 6     f(a) f0(a) f00(a)   Through a straightforward computation, we can obtain the expressions of f(a) and f0(a). By the same approach, we can obtain the expression of f0(b) , which concludes the proof. 3. Numerical Results This section provides some numerical results to illustrate the performance of the above Hermite interpolation operator. To this end, we will use the test functions f1(x) = 3 4e−2(9x−2)2−1 5e−(9x−7)2−(9x−4)2+1 2e−(9x−7)2−1 4(9x−3)2+3 4e1 10 (−9x−1)−1 49 (9x+1)2, f2(x) = 1 2xcos44x2+x−1, f3(x) = −exp−x2logx5+6+sin(3πx) cos(2πx) + 2, Mathematics 2022,10, 1490 9 of 13 whose plots appear in Figure 2. The two first functions are the 1D versions of the Franke [ 37 ] and Nielson [38] functions. 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0 1.2 0.2 0.4 0.6 0.8 1.0 0.1 0.2 0.3 0.4 0.2 0.4 0.6 0.8 1.0 -1.0 -0.8 -0.6 -0.4 -0.2 Figure 2. Plots of test functions: f1(left), f2(center) and f3(right). Let us consider the interval I= [ 0, 1 ] . The tests are carried out for a sequence of uniform mesh Inassociated with the break points ih,i=0, . . . , n, where h=1 n. The interpolation error is estimated as En(f)=max 0≤`≤200 |s(x`)−f(x`)|, where x` , `= 0, . . . , 200, are equally spaced points in I . The estimated numerical convergence order (NCO) is given by the rate NCO :=logEn E2n log(2). In Table 1, the estimated quasi-interpolation errors and NCOs for functions f1 , f2 and f3are shown. Table 1. Estimated errors for functions f1,f2and f3, and NCOs with different values of n. nEn(f1)NCO En(f2)NCO En(f3)NCO 10 1.7857 ×10−1−− 9.1243 ×10−2−− 7.4186 ×10−3−− 20 8.7411 ×10−34.3525 9.8171 ×10−33.2163 2.8348 ×10−44.7097 40 1.8198 ×10−45.5859 2.3654 ×10−45.3751 1.0365 ×10−54.7734 80 8.6397 ×10−64.3967 1.1330 ×10−54.3838 5.7600 ×10−74.1695 Now, we will compare the numerical method proposed here with the results obtained with different methods in other papers in the literature, although the test functions in these papers are extremely simple, namely g1(x) = cos(πx)and g2(x) = xsin(x). In Tables 2and 3, we list the resulting errors for the approximation of the functions g1 and g2 , respectively, by using the cubic spline operator provided here and those in References [ 34 , 39 , 40 ]. Tables 2and 3show that the novel numerical scheme improves the results in References [34,39,40]. Table 2. Estimated errors for function g1, and NCOs with different values of n. nEn(g1)NCO Method in [34] NCO Method in [39] NCO 10 3.00 ×10−5−− 6.00 ×10−5−− 6.25 ×10−4−− 20 1.86 ×10−64.01 3.28 ×10−64.19 4.05 ×10−53.94 40 1.16 ×10−74.00 2.43 ×10−73.75 2.56 ×10−63.98