C² Cubic Algebraic Hyperbolic Spline Interpolating Scheme by Means of Integral Values

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 $C^2$ everywhere and yields optimal order. We provide some numerical tests to illustrate the good performance of the novel approximation scheme.

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,6,7,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] 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 $\left\{1,x,\dots ,x^{n-2},u_1\left(x\right),u_2\left(x\right)\right\}$, where $u_1$ and $u_2$ are smooth functions. The two functions $u_1$ and $u_2$ 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: $\left(u_1\left(x\right),u_2\left(x\right)\right)=\left(sin\left(x\right),cos\left(x\right)\right)$ and $\left(u_1\left(x\right),u_2\left(x\right)\right)=\left(exp\left(x\right),exp(-x)\right)$, 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,21,22,23,24,25,26]. The algebraic hyperbolic spaces spanned by the functions 1, x, …, $x^{n-2}$, $cosh\left(x\right)$, $sinh\left(x\right)$, for $x\in \mathbb{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 $\left\{1,x,sinh,cosh\right\}$, and let $X_n={\left\{x_i:=a+ih\right\}}_{i=0}^n$ be a uniform partition of a bounded interval $\mathfrak{I}=\left[a,b\right]$, with $h=\frac{b-a}{n}$. Given values $f_i^j$, $i=0,\dots ,n$, $j=0,1$, there exists a unique cubic AH Hermite interpolant $s_i\in span\left\{1,x,sinh,cosh\right\}$ such that $s_i^{\left(j\right)}\left(x_{i+k}\right)=f_{i+k}^j$, $i=0,\dots ,n-1$, $j,k=0,1$. The spline s defined from the local interpolants $s_i$ is a $C^1$ continuous AH spline that interpolates the data $f_i^j$.

In this work, we suppose that the data $f_i^j$, $i=0,\dots ,n$, $j=0,1$, are not given, and we assume that the integral values over subintervals $[x_i,x_{i+1}]$, $i=0,\dots ,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 $f_i^j$ will be computed by means of $C^2$ smoothness conditions at the knots $x_i$ and the integral values over subintervals $[x_i,x_{i+1}]$. This is done by solving a three-diagonal linear system. Some final conditions are required. In particular, we assume that the three values $f_0^0$, $f_n^0$ and $f_0^1$ or $f_n^1$ 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 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 $C^2$ 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 $C^1$ 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 $C^2$ 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 $X_n={\left\{x_i:=a+ih\right\}}_{i=0}^n$ be a uniform partition of a bounded interval $\mathfrak{I}=\left[a,b\right]$, with $h=\frac{b-a}{n}$.

2.1. Cubic Algebraic Hyperbolic (AH) Splines of Class $C^1$

The construction of a $C^1$ cubic AH spline interpolant on the partition $X_n$ should be locally expressed in each sub-interval ${\ell }_i:=\left[x_i,x_{i+1}\right]$ in terms of the function and the first derivative values of the approximated function at knots $x_i$ and $x_{i+1}$.

The space of cubic AH splines on $X_n$ with global $C^1$ continuity is denoted as

$$S_3^1\left(X_n\right):=\left\{s\in C^1\left(\mathfrak{I}\right);s_{\mid {\ell }_i}\in {\Gamma }_3\mathrm{for}\mathrm{all}i=0,\dots ,n\right\},$$

where ${\Gamma }_3:=span\left\{1,x,sinh,cosh\right\}$ 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\left(x\right)\in S_3^1\left(X_n\right)$ such that

$$s^{\left(j\right)}\left(x_i\right)=f_i^j,i=0,\dots ,n,j=0,1,$$

(1)

for any given set of $f_i^j$-values. Therefore, each spline $s\left(x\right)\in S_3^1\left(X_n\right)$ restricted to the sub-interval ${\ell }_i$ can be represented as follows:

$$s_{\mid {\ell }_i}\left(x\right):=\sum _{k=0}^1\sum _{j=0}^1{f}_{i+k}^j{\varphi }_{i+k}^j\left(x\right),$$

(2)

in which ${\varphi }_{i+k}^j$, $k=0,1$, $j=0,1$ are classical Hermite basis functions of $S_3^1\left(X_n\right)$ restricted to the sub-interval ${\ell }_i$. The basis functions ${\varphi }_{i+k}^j$, $k=0,1$, $j=0,1$ are the unique solution of the interpolation problem given by (1) in $s\left(x\right)\in S_3^1\left(X_n\right)$, where $\left(f_i^0,f_i^1,f_{i+1}^0,f_{i+1}^1\right)=\left(1,0,0,0\right)$, $\left(f_i^0,f_i^1,f_{i+1}^0,f_{i+1}^1\right)=\left(0,1,0,0\right)$, $\left(f_i^0,f_i^1,f_{i+1}^0,f_{i+1}^1\right)=\left(0,0,1,0\right)$ and $\left(f_i^0,f_i^1,f_{i+1}^0,f_{i+1}^1\right)=\left(0,0,0,1\right)$, respectively. More precisely,

$$\begin{array}{cc}\hfill & {\varphi }_i^0\left(x_i\right)=1,{\left({\varphi }_i^0\right)}^{\prime }\left(x_i\right)=0,{\varphi }_i^0\left(x_{i+1}\right)=0,{\left({\varphi }_i^0\right)}^{\prime }\left(x_{i+1}\right)=0,\hfill \\ \hfill & {\varphi }_i^1\left(x_i\right)=0,{\left({\varphi }_i^1\right)}^{\prime }\left(x_i\right)=1,{\varphi }_i^1\left(x_{i+1}\right)=0,{\left({\varphi }_i^1\right)}^{\prime }\left(x_{i+1}\right)=0,\hfill \\ \hfill & {\varphi }_{i+1}^0\left(x_i\right)=0,{\left({\varphi }_{i+1}^0\right)}^{\prime }\left(x_i\right)=0,{\varphi }_{i+1}^0\left(x_{i+1}\right)=1,{\left({\varphi }_i^0\right)}^{\prime }\left(x_{i+1}\right)=0,\hfill \\ \hfill & {\varphi }_{i+1}^1\left(x_i\right)=0,{\left({\varphi }_{i+1}^1\right)}^{\prime }\left(x_i\right)=0,{\varphi }_{i+1}^1\left(x_{i+1}\right)=0,{\left({\varphi }_{i+1}^1\right)}^{\prime }\left(x_{i+1}\right)=1.\hfill \end{array}$$

They can be expressed as follows:

$$\begin{array}{cc}\hfill {\varphi }_i^0\left(x\right)& :=\frac{1}{\theta }(1,x,sinh\left(x\right),cosh\left(x\right))\left(\begin{array}{c}-cosh\left(h\right)+sinh\left(h\right)\left(h+x_i\right)+1\\ -sinh\left(h\right)\\ sinh\left(h+x_i\right)-sinh\left(x_i\right)\\ cosh\left(x_i\right)-cosh\left(h+x_i\right)\end{array}\right),\hfill \\ \hfill {\varphi }_i^1\left(x\right)& :=\frac{1}{\theta }(1,x,sinh\left(x\right),cosh\left(x\right))\left(\begin{array}{c}-cosh\left(h\right)-sinh\left(h\right)x_i+1\\ sinh\left(h\right)\\ sinh\left(x_i\right)-sinh\left(h+x_i\right)\\ cosh\left(h+x_i\right)-cosh\left(x_i\right)\end{array}\right),\hfill \\ \hfill {\varphi }_{i+1}^0\left(x\right)& :=\frac{1}{\theta }(1,x,sinh\left(x\right),cosh\left(x\right))\left(\begin{array}{c}hcosh\left(h\right)-sinh\left(h\right)+(cosh\left(h\right)-1)x_i\\ 1-cosh\left(h\right)\\ cosh\left(x_i\right)-cosh\left(h+x_i\right)+hsinh\left(h+x_i\right)\\ -hcosh\left(h+x_i\right)-sinh\left(x_i\right)+sinh\left(h+x_i\right)\end{array}\right),\hfill \\ \hfill {\varphi }_{i+1}^1\left(x\right)& :=\frac{1}{\theta }(1,x,sinh\left(x\right),cosh\left(x\right))\left(\begin{array}{c}-h+sinh\left(h\right)+(cosh\left(h\right)-1)x_i\\ 1-cosh\left(h\right)\\ -cosh\left(x_i\right)+cosh\left(h+x_i\right)-hsinh\left(x_i\right)\\ hcosh\left(x_i\right)+sinh\left(x_i\right)-sinh\left(h+x_i\right)\end{array}\right),\hfill \end{array}$$

where $\theta =hsinh\left(h\right)-2cosh\left(h\right)+2$.

Figure 1 shows the typical plots of the basis functions ${\varphi }_{i+k}^j$, $k=0,1$, $j=0,1$, on the sub-interval $[x_i,x_{i+1}]=[0,1]$.

Consider that the values of $f_i^j$ are related to an explicit function f, i.e., $f_i^j=f^{\left(j\right)}\left(x_i\right)$, $j=0,1$. Then, it holds that

$${\parallel s-f\parallel }_{\infty ,[x_i,x_{i+1}]}=\mathcal{O}\left(h^4\right).$$

(3)

In order to prove this statement, consider an arbitrary but fixed value $\tilde{t}$ different from $x_i$ and $x_{i+1}$, and define

$$R\left(t\right)=s\left(t\right)-f\left(t\right)-\rho \prod _{k=0}^1{\left(t-x_{i+k}\right)}^2,$$

in which the constant $\rho$ is chosen such that $R\left(\tilde{t}\right)=0$, that is,

$$\rho =\frac{s\left(\tilde{t}\right)-f\left(\tilde{t}\right)}{{\prod }_{k=0}^1{\left(\tilde{t}-x_{i+k}\right)}^2}.$$

The function R has at least three roots in $[x_i,x_{i+1}]$, which are $x_i$, $x_{i+1}$ and $\tilde{t}$.

According to Rolle’s theorem, $R^{\prime }$ has at least two roots in $[x_i,x_{i+1}]$ that are different from $x_i$, $x_{i+1}$, $\tilde{t}$ and also $R^{\prime }\left(x_i\right)=R^{\prime }\left(x_{i+1}\right)=0$, which means that $R^{\prime }$ has at least four roots in $[x_i,x_{i+1}]$. Analogously and progressively, it is shown that $R^{\left(2\right)}$ has three roots in the interval $[x_i,x_{i+1}]$, $R^{\left(3\right)}$ has two and $R^{\left(4\right)}$ has only one root, say $\xi$. It then states

$$R^{\left(4\right)}\left(\xi \right)=s^{\left(4\right)}\left(\xi \right)-f^{\left(4\right)}\left(\xi \right)-24\rho =0.$$

It holds that

$$s\left(\tilde{t}\right)-f\left(\tilde{t}\right)=\frac{1}{24}\left(s^{\left(4\right)}\left(\xi \right)-f^{\left(4\right)}\left(\xi \right)\right)\prod _{k=0}^1{\left(\tilde{t}-x_{i+k}\right)}^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 $C^2$ smoothness at the set of break points. Next, we express the derivative values in the mean integral values provided.

The $C^2$ smoothness of $s^{″}\left(x\right)$ at $x_i$, $i=1,\dots ,n-1$ yields the following consistency relations:

$$\alpha f_{i-1}^1+\beta f_i^1+\alpha f_{i+1}^1=f_{i+1}^0-f_{i-1}^0,i=1,\dots ,n-1.$$

(4)

where $\alpha =\frac{sinh\left(h\right)-h}{cosh\left(h\right)-1}$ and $\beta ={\mathrm{csch}}^2\left(\frac{h}{2}\right)(hcosh\left(h\right)-sinh\left(h\right))$.

The error related to the approximation scheme (4) can be derived from the following result.

Lemma 1.

Let $f\in {\mathcal{C}}^3\left(\left[a,b\right]\right)$; then, the local truncation errors ${\mathfrak{t}}_i$, $i=1,2,\dots ,n$, associated with the scheme (4) are given by the expressions

$${\mathfrak{t}}_i=-\frac{1}{6}{h}^2{f}^{\left(3\right)}\left(x_i\right)(h(cosh\left(h\right)+2)-3sinh\left(h\right))+\mathcal{O}\left(h^2\right).$$

Proof. 

The function f is supposed to be of class ${\mathcal{C}}^3\left(\left[a,b\right]\right)$, that is,

$$\begin{array}{cc}\hfill f\left(x_i+h\right)& =f\left(x_i\right)+f^{\prime }\left(x_i\right)h+f^{″}\left(x_i\right)\frac{h^2}{2}+f^{\left(3\right)}\left(x_i\right)\frac{h^3}{6}+\mathcal{O}\left(h^3\right),\hfill \\ \hfill f\left(x_i-h\right)& =f\left(x_i\right)-f^{\prime }\left(x_i\right)h+f^{″}\left(x_i\right)\frac{h^2}{2}-f^{\left(3\right)}\left(x_i\right)\frac{h^3}{6}+\mathcal{O}\left(h^3\right),\hfill \\ \hfill f^{\prime }\left(x_i+h\right)& =f^{\prime }\left(x_i\right)+f^{″}\left(x_i\right)h+f^{\left(3\right)}\left(x_i\right)\frac{h^2}{2}+\mathcal{O}\left(h^2\right),\hfill \\ \hfill f^{\prime }\left(x_i-h\right)& =f^{\prime }\left(x_i\right)-f^{″}\left(x_i\right)h+f^{\left(3\right)}\left(x_i\right)\frac{h^2}{2}+\mathcal{O}\left(h^2\right).\hfill \end{array}$$

Using Equation (4), it results that

$${\mathfrak{t}}_i=\alpha f^{\prime }\left(x_{i-1}\right)+\beta f^{\prime }\left(x_i\right)+\alpha f^{\prime }\left(x_{i+1}\right)-\frac{1}{{\mathrm{csch}}^2\left(\frac{h}{2}\right)}\left(f\left(x_{i+1}\right)-f\left(x_{i-1}\right)\right).$$

By replacing $f^{\prime }\left(x_{i-1}\right)$, $f^{\prime }\left(x_{i+1}\right)$ and $f\left(x_{i+1}\right)-f\left(x_{i-1}\right)$ 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\in$ ${\mathcal{C}}^3\left(\left[a,b\right]\right)$, then

$$\left|f_i^1-f^{\prime }\left(x_i\right)\right|=\mathcal{O}\left(h^2\right).$$

At this point, we have simply provided a scheme that approximates the derivative values from the function values by imposing $C^2$ 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 $\left[x_{i-1},x_i\right]$, $i=1,\dots ,n$ are provided and are equal to

$$t_i={\int }_{x_{i-1}}^{x_i}f\left(x\right)dx,i=1,\dots ,n$$

(5)

In the sequel, we will provide a scheme that approximates derivative values from the integrals $t_i$, $i=1,\dots ,n$.

By integrating s over $[x_{i-1},x_i]$ and $[x_i,x_{i+1}]$, one can obtain

$$2t_i=h\left(f_{i-1}^0+f_i^0\right)+\left(2-hcoth\left(\frac{h}{2}\right)\right)\left(f_i^1-f_{i-1}^1\right),i=1,\dots ,n$$

(6)

$$2t_{i+1}=h\left(f_i^0+f_{i+1}^0\right)+\left(2-hcoth\left(\frac{h}{2}\right)\right)\left(f_{i+1}^1-f_i^1\right),i=0,\dots ,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 $f_i^0$, as well as to achieve new relations that link only the unknowns $f_i^1$ with the provided data $t_i$. It results that

$$\mu f_{i-1}^1+\lambda f_i^1+\mu f_{i+1}^1=2\left(t_{i+1}-t_i\right),i=1,\dots ,n-1,$$

(8)

with

$$\mu =\frac{1}{2}\left(4-h^2{\mathrm{csch}}^2\left(\frac{h}{2}\right)\right)\mathrm{and}\lambda =\left(\left(h^2-2\right)cosh\left(h\right)+2\right){\mathrm{csch}}^2\left(\frac{h}{2}\right)$$

This yields a system of $n-1$ linear equations, while there are $n+1$ unknowns $f_i^1$. 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 $f^{\prime }\left(a\right)=f_a^1$ and $f^{\prime }\left(b\right)=f_b^1$ are provided. Then, a $(n-1)\times (n-1)$ linear tridiagonal system results.

$$\left(\begin{array}{cccccc}\lambda & \mu & & & & \\ \mu & \lambda & \mu & & & \\ & \mu & \lambda & \mu & & \\ & & \ddots & \ddots & \ddots & \\ & & & \mu & \lambda & \mu \\ & & & & \mu & \lambda \end{array}\right)\left(\begin{array}{c}{m}_1\\ \\ ⋮\\ \\ m_{n-1}\end{array}\right)=\left(\begin{array}{c}{b}_1-\mu c\\ b_2\\ ⋮\\ b_{n-2}\\ b_{n-1}-\mu d\end{array}\right),$$

(9)

with $b_i=2\left(t_{i+1}-t_i\right)$.

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 $\lambda >2\mu$. Indeed, a simple calculation gives us the following equality:

$$\begin{array}{cc}\hfill \lambda -2\mu & =\frac{\left(\left(h^2-2\right)cosh\left(h\right)+2\right){\mathrm{csch}}^2\left(\frac{h}{2}\right)}{h\left(hcoth\left(\frac{h}{2}\right)-2\right)}-\frac{2\left(4-h^2{\mathrm{csch}}^2\left(\frac{h}{2}\right)\right)}{2h\left(hcoth\left(\frac{h}{2}\right)-2\right)}\hfill \\ & =\frac{4}{h}+2coth\left(\frac{h}{2}\right)\hfill \end{array}$$

Since $coth\left(\frac{h}{2}\right)>0$, then $\lambda -2\mu >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 $\frac{1}{5}$ for small enough $h>0$. In fact, let $A:={\left(a_{ij}\right)}_{1\le i,j\le n}$ be the matrix coefficient of system (9). The index on the diagonally dominant property of matrix A is given by

$$\begin{array}{c}\hfill D_h\left(A\right):=\underset{i=1,\dots ,n}{max}\frac{1}{\left|a_{ii}\right|}\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^n\left|a_{i,j}\right|=2\underset{h>0}{max}\frac{\mu \left(h\right)}{\lambda \left(h\right)}=\frac{\frac{2}{h}-\frac{h}{15}+\frac{13h^3}{6300}+\mathcal{O}\left(h^5\right)}{\frac{10}{h}+\frac{4h}{15}-\frac{11h^3}{3150}+\mathcal{O}\left(h^5\right)}\end{array}$$

Its limit when h is close to zero is equal to $\frac{1}{5}$.

Once the values of $f_i^1$, $i=1,\dots ,n-1$ are determined, we can then compute $f_i^1$, $i=1,\dots ,n-1$ by means of (6). However, we still need another end condition. Suppose that one of the values $f\left(a\right)$ and $f\left(b\right)$ 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\in {\mathcal{C}}^3\left(\left[a,b\right]\right)$; then, the local truncation errors ${\tilde{\mathfrak{t}}}_i$, $i=1,2,\dots ,n$, associated with the scheme (9), are given by the expressions

$${\tilde{\mathfrak{t}}}_i=-\frac{1}{6}{h}^2{f}^{\left(3\right)}\left(x_i\right)\left(h^2+3h^2{csch}^2\left(\frac{h}{2}\right)-12\right)+\mathcal{O}\left(h^2\right).$$

Theorem 5.

Let f $\in {\mathcal{C}}^4\left(\left[a,b\right]\right)$ and s be the interpolation operator defined by (2), (6) and (9). For a uniform step size h, we have

$${∥f\left(x\right)-s\left(x\right)∥}_{\infty ,\left[a,b\right]}=\mathcal{O}\left(h^2\right).$$

Proof. 

Consider the sub-interval $[x_i,x_{i+1}]$. Let p be the spline in $S_3^1$, which satisfies $p^{\left(j\right)}\left(x_{i+k}\right)=f^{\left(j\right)}\left(x_{i+k}\right)$, $j,k=0,1$. Then, it results that

$$\begin{array}{cc}\hfill \parallel s-f\parallel & =\parallel s-p+p-f\parallel \hfill \\ \hfill & \le \parallel s-p\parallel +\parallel p-f\parallel \hfill \\ \hfill & \le \mathcal{O}\left(h^2\right)+\mathcal{O}\left(h^4\right)\hfill \\ \hfill & \equiv \mathcal{O}\left(h^2\right).\hfill \end{array}$$

This due to the fact that $\parallel s-p\parallel \le max\left\{|f^{\prime }\left(x_i\right)-f_i^1|,|f^{\prime }\left(x_{i+1}\right)-f_{i+1}^1|\right\}$, 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\left(a\right)$, $f^{\prime }\left(a\right)$ and $f\left(b\right)$ by means of integral values.

Lemma 6.

For a given function $f\in C^2$,

$$\begin{array}{cc}\hfill f\left(a\right)& =\frac{11t_1-7t_2+2t_3}{6h},f^{\prime }\left(a\right)=-\frac{2t_1-3t_2+t_3}{h^2},f^{\prime }\left(b\right)=-\frac{2t_n-3t_{n-1}+t_{n-2}}{h^2}.\hfill \end{array}$$

Proof. 

The Taylor expansion of f around a is as follows:

$$f\left(x\right)=f\left(a\right)+f^{\prime }\left(a\right)\left(x-a\right)+\frac{1}{2}{f}^{″}\left(a\right){\left(x-a\right)}^2+\mathcal{O}{\left(x-a\right)}^2.$$

$$\begin{array}{cc}\hfill t_1& =hf\left(a\right)+\frac{h^2}{2}{f}^{\prime }\left(a\right)+\frac{h^3}{6}{f}^{″}\left(a\right)+\mathcal{O}\left(h^3\right),\hfill \\ \hfill t_2& =hf\left(a\right)+\frac{3h^2}{2}{f}^{\prime }\left(a\right)+\frac{7h^3}{6}{f}^{″}\left(a\right)+\mathcal{O}\left(h^3\right),\hfill \\ \hfill t_3& =hf\left(a\right)+\frac{5h^2}{2}{f}^{\prime }\left(a\right)+\frac{19h^3}{6}{f}^{″}\left(a\right)+\mathcal{O}\left(h^3\right).\hfill \end{array}$$

Therefore, the following system is obtained:

$$\left(\begin{array}{c}{t}_1\\ t_2\\ t_3\end{array}\right)=\left(\begin{array}{ccc}h& \frac{h^2}{2}& \frac{h^3}{6}\\ h& \frac{3h^2}{2}& \frac{7h^3}{6}\\ h& \frac{5h^2}{2}& \frac{19h^3}{6}\end{array}\right)\left(\begin{array}{c}f\left(a\right)\\ f^{\prime }\left(a\right)\\ f^{″}\left(a\right)\end{array}\right)$$

Through a straightforward computation, we can obtain the expressions of $f\left(a\right)$ and $f^{\prime }\left(a\right)$.

By the same approach, we can obtain the expression of $f^{\prime }\left(b\right)$, 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

$$\begin{array}{cc}\hfill f_1\left(x\right)& =\frac{3}{4}{e}^{-2{(9x-2)}^2}-\frac{1}{5}{e}^{-{(9x-7)}^2-{(9x-4)}^2}+\frac{1}{2}{e}^{-{(9x-7)}^2-\frac{1}{4}{(9x-3)}^2}+\frac{3}{4}{e}^{\frac{1}{10}(-9x-1)-\frac{1}{49}{(9x+1)}^2},\hfill \\ \hfill f_2\left(x\right)& =\frac{1}{2}x{cos}^4\left(4\left(x^2+x-1\right)\right),\hfill \\ \hfill f_3\left(x\right)& =-\frac{exp\left(-x^2\right)\left(log\left(x^5+6\right)+sin\left(3\pi x\right)\right)}{cos\left(2\pi x\right)+2},\hfill \end{array}$$

whose plots appear in Figure 2. The two first functions are the 1D versions of the Franke [37] and Nielson [38] functions.

Let us consider the interval $I=[0,$ $1]$. The tests are carried out for a sequence of uniform mesh ${\mathfrak{I}}_n$ associated with the break points $ih$, $i=0,\dots ,$ n, where $h=\frac{1}{n}$.

The interpolation error is estimated as

$${\mathcal{E}}_n\left(f\right)=\underset{0\le \ell \le 200}{max}|s\left(x_{\ell }\right)-f\left(x_{\ell }\right)|,$$

where $x_{\ell }$, $\ell =0,\dots ,$ 200, are equally spaced points in I. The estimated numerical convergence order (NCO) is given by the rate

$$\mathrm{NCO}:=\frac{log\left(\frac{{\mathcal{E}}_n}{{\mathcal{E}}_{2n}}\right)}{log\left(2\right)}.$$

In

Table 1. Estimated errors for functions $f_1$, $f_2$ and $f_3$, and NCOs with different values of n.

n	${\mathcal{E}}_{\mathit{n}}\left({\mathit{f}}_1\right)$	NCO	${\mathcal{E}}_{\mathit{n}}\left({\mathit{f}}_2\right)$	NCO	${\mathcal{E}}_{\mathit{n}}\left({\mathit{f}}_3\right)$	NCO
10	$1.7857\times {10}^{-1}$	$--$	$9.1243\times {10}^{-2}$	$--$	$7.4186\times {10}^{-3}$	$--$
20	$8.7411\times {10}^{-3}$	$4.3525$	$9.8171\times {10}^{-3}$	$3.2163$	$2.8348\times {10}^{-4}$	$4.7097$
40	$1.8198\times {10}^{-4}$	$5.5859$	$2.3654\times {10}^{-4}$	$5.3751$	$1.0365\times {10}^{-5}$	$4.7734$
80	$8.6397\times {10}^{-6}$	$4.3967$	$1.1330\times {10}^{-5}$	$4.3838$	$5.7600\times {10}^{-7}$	$4.1695$

, the estimated quasi-interpolation errors and NCOs for functions $f_1$, $f_2$ and $f_3$ are shown.

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

$$g_1\left(x\right)=cos\left(\pi x\right)\mathrm{and}{g}_2\left(x\right)=xsin\left(x\right).$$

In

Table 2. Estimated errors for function $g_1$, and NCOs with different values of n.

n	${\mathcal{E}}_{\mathit{n}}\left({\mathit{g}}_1\right)$	NCO	Method in [34]	NCO	Method in [39]	NCO
10	$3.00\times {10}^{-5}$	$--$	$6.00\times {10}^{-5}$	$--$	$6.25\times {10}^{-4}$	$--$
20	$1.86\times {10}^{-6}$	$4.01$	$3.28\times {10}^{-6}$	$4.19$	$4.05\times {10}^{-5}$	$3.94$
40	$1.16\times {10}^{-7}$	$4.00$	$2.43\times {10}^{-7}$	$3.75$	$2.56\times {10}^{-6}$	$3.98$

and

Table 3. Estimated errors for function $g_2$, and NCOs with different values of n.

n	${\mathcal{E}}_{\mathit{n}}\left({\mathit{g}}_2\right)$	NCO	Method in [40]	NCO
10	$1.66\times {10}^{-6}$	$--$	$2.80\times {10}^{-5}$	$--$
20	$1.04\times {10}^{-7}$	$4.00$	$1.75\times {10}^{-7}$	$3.99$
40	$6.51\times {10}^{-9}$	$4.00$	$1.10\times {10}^{-8}$	$3.99$

, we list the resulting errors for the approximation of the functions $g_1$ and $g_2$, respectively, by using the cubic spline operator provided here and those in References [34,39,40]. Table 2 and Table 3 show that the novel numerical scheme improves the results in References [34,39,40].

Next, we deal with the following four test functions also defined on $[0,1]$:

$$\begin{array}{cc}\hfill & k_1\left(x\right)=\frac{\sqrt{x+2}exp\left(2x^2\right)sin\left(4\pi x\right)}{{\left(x^2+3\right)}^{5/7}},k_2\left(x\right)=\frac{sinh\left(x^2\right)sin\left(2\pi \sqrt{cosh\left(2x\right)}\right)}{x^6+1},\hfill \\ \hfill & k_3\left(x\right)=\frac{exp\left(\frac{1}{x^2+1}\right)tanh\left(\frac{x}{10\pi }\right)}{16x^3+1},\mathrm{and}{k}_4\left(x\right)=cosh\left(x\right)exp(sinh\left(x\right)).\hfill \end{array}$$

Their typical plots are shown in Figure 3.

Our goal now is to compare the method introduced herein with the methods that use only algebraic splines, as well as those that combine the features of algebraic and hyperbolic functions. To this end, we consider the approaches presented by A. Boujraf et al. in [29], D. Barrera et al. in [21], J. Wu and X. Zhang in [40] and S. Eddargani et al. in [34], so we implement the approaches described in [29,34,40] to be able to execute any test function, because those involved in the cited references are simple examples and one of them (exponential function) is reproduced by our approach.

In

Table 4. Estimated errors for function $k_1$, and NCOs with different values of n.

n	${\mathcal{E}}_{\mathit{n}}\left({\mathit{k}}_1\right)$	NCO	Method in [29]	NCO
8	$3.6083\times {10}^{-2}$	$--$	$7.8090\times {10}^{-1}$	$--$
16	$2.5592\times {10}^{-3}$	$3.81755$	$6.3999\times {10}^{-2}$	$3.60902$
32	$1.6951\times {10}^{-4}$	$3.91625$	$4.6292\times {10}^{-3}$	$3.78922$
64	$1.0783\times {10}^{-5}$	$3.9745$	$1.5127\times {10}^{-3}$	$1.61358$
128	$6.8819\times {10}^{-7}$	$3.96986$	$5.2353\times {10}^{-4}$	$1.53083$

and

Table 5. Estimated errors for function $k_2$, and NCOs with different values of n.

n	${\mathcal{E}}_{\mathit{n}}\left({\mathit{k}}_2\right)$	NCO	Method in [29]	NCO
8	$5.0763\times {10}^{-3}$	$--$	$1.2442\times {10}^{-1}$	$--$
16	$3.6283\times {10}^{-4}$	$3.80642$	$1.0854\times {10}^{-2}$	$3.51895$
32	$1.8540\times {10}^{-5}$	$4.29057$	$7.5616\times {10}^{-4}$	$3.84338$
64	$9.9072\times {10}^{-7}$	$4.22602$	$4.7513\times {10}^{-5}$	$3.99231$
128	$7.4838\times {10}^{-8}$	$3.72664$	$3.0252\times {10}^{-6}$	$3.97319$

, we list the resulting errors for the approximation of the functions $k_1$ and $k_2$, respectively, using our approach and the one provided in [29].

Table 6. Estimated errors for function $k_3$, and NCOs with different values of n.

n	${\mathcal{E}}_{\mathit{n}}\left({\mathit{k}}_3\right)$	NCO	Method in [21]	NCO	Method in [40]	NCO	Method in [34]	NCO
8	$7.78\times {10}^{-5}$	$--$	$1.24\times {10}^{-3}$	$--$	$1.64\times {10}^{-2}$	$--$	$1.69\times {10}^{-2}$	$--$
16	$1.93\times {10}^{-6}$	$5.33$	$7.81\times {10}^{-5}$	$3.98$	$1.00\times {10}^{-3}$	$4.03$	$1.01\times {10}^{-3}$	$4.27$
32	$1.03\times {10}^{-7}$	$4.21$	$4.88\times {10}^{-6}$	$4.00$	$6.26\times {10}^{-5}$	$3.99$	$6.28\times {10}^{-5}$	$4.00$
64	$6.02\times {10}^{-9}$	$4.10$	$3.05\times {10}^{-7}$	$4.00$	$3.91\times {10}^{-6}$	$4.00$	$3.92\times {10}^{-6}$	$4.00$
128	$4.91\times {10}^{-10}$	$3.61$	$1.90\times {10}^{-8}$	$4.00$	$2.44\times {10}^{-7}$	$4.00$	$2.45\times {10}^{-7}$	$4.00$

shows the resulting errors for the approximation of the function $k_3$, using the approach described here and those in References [21,34,40]. In

Table 7. Estimated errors for function $k_4$, and NCOs with different values of n.

n	${\mathcal{E}}_{\mathit{n}}\left({\mathit{k}}_3\right)$	NCO	Method in [40]	NCO	Method in [34]	NCO
8	$9.41\times {10}^{-5}$	$--$	$3.25\times {10}^{-4}$	$--$	$3.83\times {10}^{-4}$	$--$
16	$7.70\times {10}^{-6}$	$3.61$	$2.69\times {10}^{-5}$	$3.59$	$3.16\times {10}^{-5}$	$3.59$
32	$5.19\times {10}^{-7}$	$3.88$	$1.95\times {10}^{-6}$	$3.78$	$2.29\times {10}^{-6}$	$3.79$
64	$3.06\times {10}^{-8}$	$4.08$	$1.32\times {10}^{-7}$	$3.88$	$1.54\times {10}^{-7}$	$3.90$

, we list the resulting errors for the approximation of the function $k_4$ by using the approach provided here and those in References [34,40].

It is clear that the proposed scheme improves the results in the previous papers by at least two orders of magnitude. Consequently, in certain contexts, it is highly recommended to use splines that benefit from the features of both algebraic and hyperbolic functions instead of using only the features of algebraic ones.

4. Conclusions

Approximation from integral values represents a critical topic because of its extensive application in many different areas. This paper considers a cubic Hermite spline interpolant that reproduces linear polynomials and hyperbolic functions. The proposed interpolant is $C^2$ everywhere and is defined from the value and the first derivative value at each knot of the partition. The function and derivative values are assumed to be unknowns, and then they are determined using the $C^2$ smoothness conditions and the mean integral values. The numerical results illustrate the good performance of the novel approximation scheme. The construction used herein requires the resolution of a system of linear equations, which can be computationally expensive, especially when dealing with a large number of data. Future work will address the issue of avoiding this limitation.