Adapting Cubic Hermite Splines to the Presence of Singularities Through Correction Terms

Abstract

Hermite interpolation is classically used to reconstruct smooth data when the function and its first order derivatives are available at certain nodes. If first order derivatives are not available, it is easy to set a system of equations imposing some regularity conditions at the data nodes in order to obtain them. This process leads to the construction of a Hermite spline. The problem of the described Hermite splines is that the accuracy is lost if the data contains singularities. The consequence is the appearance of oscillations, if there is a jump discontinuity in the function, that globally affects the accuracy of the spline, or the smearing of singularities, if the discontinuities are in the derivatives of the function. This paper is devoted to the construction and analysis of a new technique that allows for the computation of accurate first order derivatives of a function close to singularities using a Hermite spline. The idea is to correct the system of equations of the spline in order to attain the desired accuracy even close to the singularities. Once we have computed the first order derivatives with enough accuracy, a correction term is added to the Hermite spline in the intervals that contain a singularity. The aim is to reconstruct piecewise smooth functions with $O(h^4)$ accuracy even close to the singularities. The process of adaption will require some knowledge about the position of the singularity and the jumps of the function and some of its derivatives at the singularity. The whole process can be used as a post-processing, where a correction term is added to the classical cubic Hermite spline. Proofs for the accuracy and regularity of the corrected spline and its derivatives are given. We also analyse the mechanism that eliminates the Gibbs phenomenon close to jump discontinuities in the function. The numerical experiments presented confirm the theoretical results obtained.

1 Introduction

In the classic literature, we can find that a spline can be constructed defining polynomial pieces that join together under certain regularity or differentiability requirements at these joints. Such regularity properties are the reason why splines are used very often in computer aided geometric design (CAGD) for industrial and professional design applications. Some examples are the representation of geometrical objects in one or several dimensions [1,2,3,4,5], the solution of PDEs or ODEs [6,7,8], image processing [9,10,11] and many others.

Sometimes the data used to construct the spline presents singularities. In the cases where the discontinuity is in the function, some kinds of splines are known to introduce oscillations close to the discontinuities [12, 13]. These oscillations classically appear when truncating the Fourier series of a discontinuous function and have been widely studied [14,15,16,17]. Global expansions, such as Fourier series, are affected by the presence of local discontinuities. The usual result is non-uniform convergence and oscillations of the partial sums. D. Gottlieb and C.-W. Shu showed that it is possible to recover high order information from these slowly and non-uniformly convergent global approximations [18,19,20,21,22]. Some other interesting and related works are [16, 23, 24]. By splitting a function with singularities in two parts, a smooth one and another one containing the information of the singularity, K. S. Eckhoff [25, 26] presented a way of modifying the Fourier’s method to obtain numerical techniques that allow us to compute derivatives and integrals with a high order of accuracy.

In this work, our approach is different from all previous techniques, although the philosophy might be similar: to reconstruct piecewise smooth functions from discretized data, taking into account the presence of singularities. In our case, we approach the problem using cubic Hermite splines plus correction terms in order to attain adaption close to the singularities.

Recently, we have introduced some modified cubic splines that aim to solve the problems caused by the oscillations that appear close to the discontinuities [27, 28]. In those previous papers we try to exploit the properties of a nonlinear mean in the construction of the spline, through which we automatically adapt the interpolation to the presence of discontinuities in the function, thus eliminating the Gibbs oscillations. One drawback is that this technique cannot provide full accuracy close to the singularity, and it smears the discontinuities. In this article, it is also our aim to design a nonlinear cubic Hermite spline interpolation adapted to the presence of singularities and constructed using the point-values discretization, i.e. using the data values f(x) at the positions $x=x_j, j=0,\ldots , m$. In this case, we try to obtain fully accurate first order derivatives and interpolation even close to the singularities (by full accuracy we mean the highest theoretically possible accuracy that we can obtain through our spline), but trying to preserve the jumps in the function and its derivatives in the reconstruction. The resulting method presents a high order of accuracy in the whole domain and provides piecewise $C^2$ regularity at both sides of the singularities. In order to reach our objective, we are inspired by the construction of the Immersed Interface Method (IIM) [29,30,31], in the sense that we design correction terms that allow us to fulfil our objectives. This method has been mostly used for the solution of elliptic equations with singularities in the context of fluid–structure interaction and it is based on an accurate tracking of the singularities and the modification of the finite difference scheme close to them.

Let us start by describing in brief detail the discretization of data that we shall use. The point-values sampling process used to obtain the pairs of data $(x_j, f(x_j)), j=0, \ldots , m$, from a function f that may contain jump discontinuities implies the loss of the information regarding the exact position of these jump discontinuities, as it only preserves local information at the positions $x_j$. Even so, it is possible to detect the interval that contains the discontinuity [32]. The described kind of discretization would only allow for the detection and location of the position of discontinuities in the first order derivative (kinks in the function) [32, 33]. Yet, if we assume that the data come from a discretization through local averages of the function f(x), i.e. the cell-averages setting [32], then it is indeed possible to locate jump discontinuities in f(x). In this work we use data discretized in the point-values. Further research could cover the results using cell-averages.

Let us continue by describing the classical Hermite interpolation and by exposing some classical results. Classical cubic Hermite interpolation is based on the construction of a cubic polynomial using the values $f_j=f(x_j), f_{j+1}=f(x_{j+1})$ of the function f, and its first order derivatives, $f_j'=f'(x_j)$ at the positions $x_j, x_{j+1}$. These four pieces of data allow to obtain the four unknowns in a cubic polynomial. Then, continuity conditions on the function and its first and second order derivatives are imposed at the $x_j$ in order to set a system of equations that provides a global solution for the considered interval. That way, piecewise cubic Hermite interpolation is known to be fourth order accurate for smooth functions. Assuming that the data is smooth, if the spline uses approximations of the first order derivatives $\tilde{f}'_j\approx f'_j$, (from now on, if not stated otherwise, we use the tilde $\tilde{}$ to represent approximations of different values: the function, the derivatives, the location of the singularity, etc.), then the maximum accuracy of the piecewise cubic Hermite interpolant is known to be $O(h^4)$ if the approximation of the first order derivatives is $O(h^3)$ accurate. If f presents singularities between the nodes (from here on we will use the word node as a synonym for the points of union of the different polynomial pieces of the spline), even good approximations of the first order derivatives at the nodes will lead to inaccurate approximations of the function inside the interval that contains the singularity. Discontinuities in the data lead either to classical Hermite splines showing Gibbs oscillations close to discontinuities in the function, or the smearing of singularities if the discontinuity is in the derivatives. This is logical considering that Hermite interpolation can be expressed in terms of the Hermite basis, represented in Fig. 1, which is composed of polynomials that are smooth, two of them being non-monotone. The smearing of discontinuities in the first order derivative can be clearly explained through the smoothness of the base. The occurrence of oscillations close to jumps in the function will be justified later on, showing that the non-monotone elements of the base are multiplied by non-zero coefficients close to discontinuities.

Thus, our aim is twofold: to design a technique that allows for the obtention of the first order derivatives of the function at the nodes of the spline with full accuracy even close to singularities, and to use the computed first order derivatives to obtain a fully accurate spline. These two objectives will be reached through the computation of correction terms.

This manuscript is organized as follows. First of all, in Sect. 2 we introduce the classical way of constructing cubic Hermite splines. Section 3 explains how to obtain adapted first order derivatives using a Hermite spline. There we present the first main result of this article: a theorem about the accuracy of the adapted first order derivatives. Section 4 presents a study about the elimination of the Gibbs phenomenon in the classical spline when using the corrected first order derivatives. Section 5 introduces an adapted Hermite spline and analyses theoretically the accuracy of the interpolation near singularities, which is the second main result of this work. Section 6 exposes how the correction terms can be used as a post-processing of the classical cubic Hermite spline. Section 7 presents some numerical experiments which show how the new algorithm performs using univariate functions. In particular, experiments about the accuracy and regularity of the function and the two first order derivatives are presented, jointly with some tests that show the elimination of the Gibbs phenomenon close to jump discontinuities in the function. Finally, Sect. 8 presents the conclusions.

2 Some Preliminaries About Classical Cubic Hermite Splines

In this section we briefly introduce the classical way of constructing the cubic Hermite spline. Then, we introduce some new techniques that allow for the adaption of the classical interpolant to the presence of singularities. The objective is to design a technique that enables the obtention of sharp reconstructions of functions with discontinuities in the function or the derivatives, avoiding the smearing of the singularity or the appearance of Gibbs oscillations.

There is an extensive list of references that treat the field of the construction and the study of the properties of cubic Hermite splines. For example, [34] or [35] contain a in-depth revision of the classical theoretical results in the field. The idea behind the construction of the spline is to compute a piecewise polynomial function of degree three between the data nodes. The resulting function must be $C^2$, i.e. the function and the first two derivatives must be continuous. The basis used to construct the polynomial between the nodes varies in the literature. Some authors start from $m+1$ pairs of values $(x_j, y_j), j=0, \ldots , m$ and write the expression of the polynomial at a particular interval $[x_{j}, x_{j+1}]$ as

$$\begin{aligned} g_j(x)={a}_j(x-x_j)^3+{b}_j(x-x_j)^2+{c}_j(x-x_j)+{d}_j. \end{aligned}$$

(1)

Sometimes it is more convenient to use the alternative expression

$$\begin{aligned} g_j(x)={a}_j(x-x_j)^2(x-x_{j+1})+{b}_j(x-x_j)^2+{c}_j(x-x_j)+{d}_j, \end{aligned}$$

(2)

which easily provides a bound for the error of the Hermite interpolation. In any of both cases, regularity conditions must be imposed at the nodes in order to obtain a piecewise defined function that is $C^2$. Basically, we need to impose the continuity of the function and the two first order derivatives at the nodes. If $y_j, y_{j+1}$ denote the values of the function at the nodes $x_j, x_{j+1}$, and $D_j, D_{j+1}$ the values of the first order derivatives at the same nodes, imposing the continuity of the function and the first and second order derivatives at the nodes, it is straightforward to obtain the following expression for the coefficients in (1) (see for example [27] for a complete deduction of the formulas),

$$\begin{aligned} \begin{aligned} a_j&=\frac{h_{j+1}D_{j+1}+D_j h_{j+1}+2y_j-2y_{j+1}}{h_{j+1}^3},\\ b_j&=-\frac{h_{j+1}D_{j+1}+2D_j h_{j+1}+3y_j-3y_{j+1}}{h_{j+1}^2},\\ c_j&=D_j,\\ d_j&=y_j, \end{aligned} \end{aligned}$$

(3)

where $h_{j+1}=x_{j+1}-x_j$.

As mentioned before, the expression of the Hermite spline (1) can be expressed in terms of the Hermite basis, represented in Fig. 1 for cubic splines. By replacing the coefficients of the spline in (1) by the values found in (3) we obtain,

$$\begin{aligned} \begin{aligned} g_j(x)=&\,\frac{h_{j+1}D_{j+1}+D_j h_{j+1}+2y_j-2y_{j+1}}{h_{j+1}^3}(x-x_j)^3\\&-\frac{h_{j+1}D_{j+1}+2D_j h_{j+1}+3y_j-3y_{j+1}}{h_{j+1}^2}(x-x_j)^2+D_j(x-x_j)+y_j. \end{aligned} \end{aligned}$$

(4)

Moreover, the change of variables $s=\frac{x-x_j}{h_{j+1}}$ returns the expression of the spline in terms of the cubic Hermite basis,

$$\begin{aligned} { \begin{aligned} g_j(s)=&\, \left( h_{j+1}D_{j+1}+D_j h_{j+1}+2y_j-2y_{j+1}\right) s^3\\&-\left( h_{j+1}D_{j+1}+2D_j h_{j+1}+3y_j-3y_{j+1}\right) s^2+\,D_j h_{j+1} s+y_j\\ =&\, \left( 1-s^2 (3-2 s)\right) y_j+s^2 (3-2 s) y_{j+1}+\left( s^3-2s^2+s\right) h_{j+1}D_j\\&+\left( s^3-s^2\right) h_{j+1}D_{j+1}\\ =&\,{b_1(s) y_{j+1}+b_2(s) y_j}+b_3(s)h_{j+1}D_j+b_4(s)h_{j+1}D_{j+1}, \end{aligned} } \end{aligned}$$

(5)

for $s\in [0, 1]$.

Fig. 1

Representation of the Hermite basis for $s\in [0, 1]$

The values of the first order derivatives at the nodes $D_j$ can be known a priori, as with the Hermite interpolation, or they can be obtained by imposing the continuity of the second order derivatives at the nodes, as with the Hermite spline. Thus, the Hermite spline relies on the solution of a linear system of equations for the $D_j$, which needs two boundary conditions. Common options for the boundary conditions that can be found in the literature are: natural boundary conditions, not-a-knot condition, complete cubic spline, etc. The complete cubic spline consists of imposing slope boundary conditions,

$$\begin{aligned} g'_0(x_0)=D_0,\quad g'_{m-1}(x_m)=D_m, \end{aligned}$$

and we have chosen it for this article. If the first order derivatives are not available at the boundaries, they can be replaced by $O(h^3)$ approximations.

Imposing the continuity of the second order derivatives of the spline leads to the following equation for each interval between the nodes

$$\begin{aligned} \frac{D_{j-1}}{h_j}+2\left( \frac{1}{h_j}+\frac{1}{h_{j+1}}\right) D_{j}+\frac{D_{j+1}}{h_{j+1}}=3\left( \frac{y_{j+1}-y_{j}}{h_{j+1}^2}+\frac{y_{j}-y_{j-1}}{h_{j}^2}\right) . \end{aligned}$$

Thus we easily obtain a linear system for the first order derivatives $D_j$ at the nodes,

$$\begin{aligned} \begin{aligned}&\left[ \begin{array}{lclclclclcl} 2\left( \frac{1}{h_1}+\frac{1}{h_{2}}\right) &{}\frac{1}{h_2}&{}0&{}\cdots &{}0&{}0\\ \frac{1}{h_2}&{}2\left( \frac{1}{h_2}+\frac{1}{h_{3}}\right) &{}\frac{1}{h_3}&{}\cdots &{}0&{}0\\ \cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots \\ 0&{}0&{}0&{}\cdots &{}\frac{1}{h_{m-1}} &{}2\left( \frac{1}{h_{m-1}}+\frac{1}{h_{m}}\right) \end{array}\right] \left[ \begin{array}{lc|} D_1\\ D_2\\ \cdots \\ D_{m-1} \end{array}\right] \\&\quad =3\left[ \begin{array}{lc|} \frac{\delta _2}{h_2}+\frac{\delta _1}{h_1}-\frac{D_0}{3 h_1}\\ \frac{\delta _3}{h_3}+\frac{\delta _2}{h_2}\\ \cdots \\ \frac{\delta _m}{h_m}+\frac{\delta _{m-1}}{h_{m-1}}-\frac{D_m}{3 h_m} \end{array}\right] , \end{aligned} \end{aligned}$$

(6)

with $\delta _j=\frac{y_j-y_{j-1}}{h_{j}}$. We will consider a uniform grid spacing, but the results for non-uniform grid spacing can be obtained in an analogous way. With this consideration, the system in (6) transforms into

$$\begin{aligned} \begin{aligned} \left[ \begin{array}{lllllll} 4&{}1&{}0&{}0&{}\cdots &{}0&{}0\\ 1&{}4&{}1&{}0&{}\cdots &{}0&{}0\\ \cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots \\ 0&{}0&{}0&{}0&{}\cdots &{}1&{}4 \end{array}\right] \left[ \begin{array}{l} D_1\\ D_2\\ \cdots \\ D_{m-1}\\ \end{array}\right] =6\left[ \begin{array}{l} \frac{\delta _2+ \delta _1}{2}-\frac{D_0}{6}\\ \frac{\delta _3+\delta _2}{2}\\ \cdots \\ \frac{\delta _{m}+ \delta _{m-1}}{2}-\frac{D_m}{6}\\ \end{array}\right] . \end{aligned} \end{aligned}$$

(7)

The existence of a singularity at any place of the domain leads to the computation of inaccurate first order derivatives through the system in (7). In case of jump discontinuities in the function, the consequence is the appearance of global Gibbs oscillations in the reconstruction of the spline that will affect the entire domain. A discussion about the size of the oscillations close to the discontinuity can be found in [27]. If a discontinuity is found in the derivatives, the singularity is smeared. In the next sections we will propose correction terms for preserving the accuracy of the spline close to the singularities.

We can also consider now the accuracy of the Hermite interpolation and its second order derivatives. These are classical results that will be used later, and which proofs can be found in pages 58 and 59 of [36].

Theorem 1

(Error of classical cubic Hermite interpolation)

Given a sufficiently smooth function f, the error for the cubic Hermite interpolation is given by

$$\begin{aligned} |f(x)-g_j(x)|\le \displaystyle \frac{\max |f^{(4)}(\xi )|}{384}h^4,\quad x,\xi \in [x_j,x_{j+1}], \end{aligned}$$

(8)

for all $h>0$.

Corollary 1

(Error of the second order derivative of classical cubic Hermite interpolation)

Given a sufficiently smooth function f, the error for the second order derivative of the cubic Hermite interpolation is given by

$$\begin{aligned} |f''(x)-g''_j(x)|\le \frac{\max |f^{(4)}(\xi )|}{2}h^2,\quad x, \xi \in [x_j,x_{j+1}], \end{aligned}$$

(9)

for all $h>0$.

Remark 1

If the first order derivatives $D_j, D_{j+1}$ are approximated, an accuracy of $O(h^3)$ is enough to keep the $O(h^4)$ accuracy of the cubic Hermite interpolation at smooth zones. To confirm this observation, instead of considering the canonical base as in (1), let us write the Hermite polynomial in the interval $[x_j,x_{j+1}]$ using the expression in (2). As we know the values $f(x_j)=y_j, f(x_{j+1})=y_{j+1}$ and the first order derivatives of the function $f'(x_j)=D_j, f'(x_{j+1})=D_{j+1}$, we can just set a system of equations to obtain the Hermite polynomial. It takes the following expression:

$$\begin{aligned} g_j(x)= & {} \frac{(D_j+D_{j+1})h+2(y_{j+1}-y_j)}{h^3}(x-x_j)^2(x-x_{j+1})\nonumber \\{} & {} +\frac{y_{j+1}-y_j-D_jh}{h^2}(x-x_j)^2+D_j(x-x_j)+y_{j}. \end{aligned}$$

(10)

It is clear from the expression in (10) that $O(h^3)$ is enough in $D_j, D_{j+1}$ to preserve the $O(h^4)$ accuracy of the cubic Hermite interpolation.

3 Accurate Approximation of First Order Derivatives Using Splines and Data in the Point-Values

Fig. 2

In this figure we can see an example of discontinuity in the function and its derivatives placed in the interval $[x_j, x_{j+1}]$ at a position $x^*=x_j+\alpha$

Let us assume that there is a singularity at $x^*=x_j+\alpha$ in the interval $[x_j, x_{j+1}]$, just as shown in Fig. 2. For the moment, let us assume that we know the position of the singularity and the jumps in the function and its derivatives. As shown in Fig. 2, we label all the information about the function to the left of the singularity as −. The information to the right is labeled as $+$. Under this configuration, having the values of a sufficiently piecewise smooth function f at the nodes $x_0, x_1\ldots , x_m$ and the jumps in the function and some of its derivatives at $x^*$, $[f], [f'], [f''], [f''']$, it is possible to correct the approximations of the first and second order derivatives of the function f computed through the cubic spline (4) at the same nodes $x_0, x_1\ldots , x_m$ to obtain $O(h^3)$ accuracy. Of course, if there is a singularity, the system in (7) must be adapted to its presence in order to fulfil the accuracy and regularity requirements. In order to do this, we can modify the right hand side of the system in (7). For the general case presented in Fig. 2, we can adapt the system in (7) by modifying the equations corresponding to $D_j$ and $D_{j+1}$. At smooth zones of the data, it is clear that the system in (7) must hold, as the regularity conditions imposed at the nodes are met. If a singularity is placed in the interval $[x_{j}, x_{j+1}]$, the system in (7) is not valid, as the function cannot be approximated accurately through a polynomial. As each equation of the system in (7) is designed to obtain a particular first order derivative $D_j$ at the point $x_j$, if we encounter a singularity, we can just add a correction term at the right hand side of the system that takes into account the presence of the singularity and that allows for the obtention of the desired accuracy of the first order derivative. Let us summarize all the previous considerations in the following Lemma.

Lemma 1

If the singularity is placed in the interval $[x_{j}, x_{j+1}]$, then the local truncation error $C_j$ of the equation for the first order derivatives $D_j$ is equal to

$$\begin{aligned} \begin{aligned} C_j=&\,\left( 3\,{\frac{[f]}{h}}-\frac{1}{24}\,{\frac{ \left( -48\,h+72\,\alpha \right) [f']}{h}}-\frac{1}{24}\,{\frac{ \left( -36\,{\alpha }^{2}+48\, \alpha \,h-12\,{h}^{2} \right) [f'']}{h}}\right. \\&\left. -\,\frac{1}{24}\,{\frac{ \left( 12 \,{\alpha }^{3}+12\,{h}^{2}\alpha -24\,{\alpha }^{2}h \right) [f'''] }{h}}\right) -\frac{1}{24}\,{\frac{ \left( -6\,{h}^{2}{\alpha }^{2}-3\,{\alpha }^{4}+{h}^{4}+8\,{\alpha }^{3}h \right) [f^{(4)}]}{h}}\\&+O(h^4),\\ \end{aligned} \end{aligned}$$

(11)

and for the first order derivative $D_{j+1}$ it is equal to

$$\begin{aligned} \begin{aligned} C_{j+1}=&\,\left( 3\,{\frac{[f]}{h}}+\frac{1}{24}\,{\frac{ \left( -72\,\alpha +24\,h \right) [f']}{h}}+\frac{1}{24}\,{\frac{ \left( -24\,\alpha \,h+36\,{ \alpha }^{2} \right) [f'']}{h}}\right. \\&\left. +\,\frac{1}{24}\,{\frac{ \left( 12\,{\alpha } ^{2}h-12\,{\alpha }^{3} \right) [f''']}{h}}\right) +\frac{1}{24}\,{\frac{ \left( 3\,{\alpha }^{4}-4\,{\alpha }^{3}h \right) [f^{(4)}]}{h}}+O(h^4). \end{aligned} \end{aligned}$$

(12)

Proof

The equations for the first order derivatives $D_{j}^-$ and $D_{j+1}^+$ in (7) are

$$\begin{aligned} \begin{aligned} D_{j-1}^-+4D_j^-+D_{j+1}^+&=6\frac{f^+_{j+1}-f^-_{j-1}}{2h},\\ D_{j}^-+4D_{j+1}^++D_{j+2}^+&=6\frac{f^+_{j+2}-f^-_{j}}{2h}. \end{aligned} \end{aligned}$$

(13)

Taking into account the presence of the singularity, using Taylor expansions and considering the resulting local truncation error, we can write the expression in (13) using quantities only from one side of the discontinuity. Let us denote by $C^L_{j}$ the correction term for the left hand side of (13) and by $C^R_{j}$ the correction term for the right hand side of (13),

$$\begin{aligned} \begin{aligned} D_{j-1}^-+4D_j^-+D_{j+1}^+&=D_{j-1}^-+4D_j^-+D_{j+1}^-+C^L_j+O(h^4)=\frac{f^+_{j+1}-f^-_{j-1}}{2h}\\&=\frac{f^-_{j+1}-f^-_{j-1}}{2h}+C^R_j+O(h^4),\\ D_{j}^-+4D_{j+1}^++D_{j+2}^+&=D_{j}^++4D_{j+1}^++D_{j+2}^++C^L_{j+1}+O(h^4)=\frac{f^+_{j+2}-f^-_{j}}{2h}\\&=\frac{f^+_{j+2}-f^+_{j}}{2h}+C^R_{j+1}+O(h^4),\\ \end{aligned} \end{aligned}$$

(14)

so we can simply consider that, due to the presence of the singularity,

$$\begin{aligned} \begin{aligned} D_{j-1}^-+4D_j^-+D_{j+1}^-&=6\frac{f^-_{j+1}-f^-_{j-1}}{2h}+C_j+O(h^4),\\ D_{j}^++4D_{j+1}^++D_{j+2}^+&=6\frac{f^+_{j+2}-f^+_{j}}{2h}+C_{j+1}+O(h^4),\\ \end{aligned} \end{aligned}$$

(15)

where $C_j=C^R_j-C^L_j$ and $C_{j+1}=C^R_{j+1}-C^L_{j+1}$ amount to the local truncation error. These terms will depend on the position of the singularity and the jumps in the function and its derivatives. Thus, we will need the jump relations

$$\begin{aligned} \begin{aligned} \left[ f\right]&=f^{+}(x^{*}) - f^{-}(x^{*}),\\ \left[ f' \right]&=f_x^{+}(x^{*}) - f_x^{-}(x^{*}),\\ \left[ f'' \right]&=f_{xx}^{+}(x^{*}) - f_{xx}^{-}(x^{*}),\\ \left[ f''' \right]&=f_{xxx}^{+}(x^{*}) - f_{xxx}^{-}(x^{*}). \end{aligned} \end{aligned}$$

(16)

Looking at Fig. 2 and (13), we can see that, in the equation for $D_j^-$, $f_{j+1}^+$ and $D_{j+1}^+$ belong to the $+$ side of the domain, while $D_j^-$ belongs to the − side. Thus, we can use the Taylor expansions of $f_{j+1}^+$ and $D_{j+1}^+$ around $x^*$ and then use the jump relations $[f], [f'], [f''], [f''']$ in (16) to express the $+$ values in terms of the − values. Assuming that we know the jump conditions in (16), or that we have a good approximation, we are ready to obtain expressions for $f^+_{j+1}$ in terms of $f^-_{j+1}$ and for $f^-_{j}$ in terms of $f^+_{j}$. Using Taylor expansions, we can write

$$\begin{aligned} \begin{aligned} f^-(x_j)&=f^-_j=f^-(x^*)-f^-_x(x^*) \alpha +\frac{1}{2}f^-_{xx}(x^*) \alpha ^2-\frac{1}{3!}f^-_{xxx}(x^*) \alpha ^3+O(h^4),\\ f^+(x_j)&=f^+_j=f^+(x^*)-f^+_x(x^*) \alpha +\frac{1}{2}f^+_{xx}(x^*) \alpha ^2-\frac{1}{3!}f^+_{xxx}(x^*) \alpha ^3+O(h^4),\\ f^-(x_{j+1})&=f^-_{j+1}=f^-(x^*)+f^-_x(x^*) (h-\alpha ) +\frac{1}{2}f^-_{xx}(x^*) (h-\alpha )^2\\&\quad +\frac{1}{3!}f^-_{xxx}(x^*) (h-\alpha )^3+O(h^4),\\ f^+(x_{j+1})&=f^+_{j+1}=f^+(x^*)+f^+_x(x^*) (h-\alpha ) +\frac{1}{2}f^+_{xx}(x^*) (h-\alpha )^2\\&\quad +\frac{1}{3!}f^+_{xxx}(x^*) (h-\alpha )^3+O(h^4), \end{aligned} \end{aligned}$$

(17)

and subtracting we obtain

$$\begin{aligned} \begin{aligned} f_j^+&=f_j^-+[f]-[f']\alpha +\frac{1}{2}[f''] \alpha ^2-\frac{1}{3!}[f'''] \alpha ^3+O(h^4),\\ f_{j+1}^+&=f_{j+1}^-+[f]+[f'](h-\alpha )+\frac{1}{2}[f''] (h-\alpha )^2+\frac{1}{3!}[f'''] (h-\alpha )^3+O(h^4). \end{aligned} \end{aligned}$$

(18)

By replacing these expressions in the right hand side of (13), we obtain the local truncation error for this part of the equation:

$$\begin{aligned} \frac{f^+_{j+1}-f^-_{j-1}}{2h}&=\frac{f^-_{j+1}-f^-_{j-1}}{2h}+\frac{1}{2}\,{\frac{[f]}{h}}-\frac{1}{2}\,{\frac{ \left( \alpha -h \right) [f']}{h}}+\frac{1}{4}\,{\frac{ \left( {h}^{2}-2\,h\alpha +{\alpha }^{2} \right) [f'']}{h}}\nonumber \\&\quad -\frac{1}{12}\,{\frac{ \left( -{h}^{3}+3\,{h}^{2} \alpha -3\,h{\alpha }^{2}+{\alpha }^{3} \right) [f''']}{h}}\nonumber \\&\quad +\frac{1}{48}\,{ \frac{ \left( {h}^{4}-4\,{h}^{3}\alpha +6\,{h}^{2}{\alpha }^{2}-4\,h{ \alpha }^{3}+{\alpha }^{4} \right) [f^{(4)}]}{h}}+O(h^4)\nonumber \\&=\frac{f^-_{j+1}-f^-_{j-1}}{2h}+C^R_j+O(h^4),\nonumber \\ \frac{f^+_{j+2}-f^-_{j}}{2h}=&\,\frac{f^+_{j+2}-f^+_{j}}{2h}+\frac{1}{2}\,{\frac{[f]}{h}}-\frac{1}{2}\,{\frac{[f']\,\alpha }{h}}+\frac{1}{4}\,{ \frac{[f'']\,{\alpha }^{2}}{h}}-\frac{1}{12}\,{\frac{[f''']\,{ \alpha }^{3}}{h}}\nonumber \\&\quad +\frac{1}{48}\,{\frac{[f^{(4)}]\,{\alpha }^{4}}{h}}+O(h^4)\nonumber \\&=\frac{f^+_{j+2}-f^+_{j}}{2h}+C^R_{j+1}+O(h^4). \end{aligned}$$

(19)

The same thing can be done for the left hand side of (13). We can write that

$$\begin{aligned} f_x^-(x_j)&=f^-_x(x^*)-f^-_{xx}(x^*) \alpha +\frac{1}{2}f^-_{xxx}(x^*) \alpha ^2-\frac{1}{3!}f^-_{xxxx}(x^*) \alpha ^3+O(h^4), \nonumber \\ f_x^+(x_j)&=f_x^+(x^*)-f^+_{xx}(x^*) \alpha +\frac{1}{2}f^+_{xxx}(x^*) \alpha ^2-\frac{1}{3!}f^+_{xxxx}(x^*) \alpha ^3+O(h^4),\nonumber \\ f_x^-(x_{j+1})&=f_x^-(x^*)+f^-_{xx}(x^*) (h-\alpha ) +\frac{1}{2}f^-_{xxx}(x^*) (h-\alpha )^2\nonumber \\&\quad +\frac{1}{3!}f^-_{xxxx}(x^*) (h-\alpha )^3+O(h^4),\nonumber \\ f_x^+(x_{j+1})&=f_x^+(x^*)+f^+_{xx}(x^*) (h-\alpha ) ++\frac{1}{2}f^+_{xxx}(x^*) (h-\alpha )^2\nonumber \\&\quad +\frac{1}{3!}f^+_{xxxx}(x^*) (h-\alpha )^3+O(h^4). \end{aligned}$$

(20)

By subtracting and denoting $D_j^-=f_x^-(x_j), D_j^+=f_x^+(x_j), D_{j+1}^-=f_x^-(x_{j+1})$ and $D_{j+1}^+=f_x^+(x_{j+1})$ we obtain

$$\begin{aligned} \begin{aligned} D_j^+&=D_j^-+[f']-[f'']\alpha +\frac{1}{2}[f'''] \alpha ^2-\frac{1}{3!}[f^{(4)}] \alpha ^3+O(h^4),\\ D_{j+1}^+&=D_{j+1}^-+[f']+[f''](h-\alpha )+\frac{1}{2}[f'''] (h-\alpha )^2+\frac{1}{3!}[f^{(4)}] (h-\alpha )^3+O(h^4). \end{aligned} \end{aligned}$$

(21)

By replacing now (21) in the left hand side of (13), we get expressions for $C^L_j$ and $C^L_{j+1}$:

$$\begin{aligned} \begin{aligned} D_{j-1}^-+4D_j^-+D_{j+1}^+&=D_{j-1}^-+4D_j^-+D_{j+1}^--[f']+[f'']\alpha -\frac{1}{2}[f'''] \alpha ^2\\&\quad +\frac{1}{3!}[f^{(4)}] \alpha ^3+O(h^4)\\&=D_{j-1}^-+4D_j^-+D_{j+1}^-+C^L_j+O(h^4),\\ D_{j}^-+4D_{j+1}^++D_{j+2}^+&=D_{j}^++4D_{j+1}^++D_{j+2}^++[f']+[f''](h-\alpha )+\frac{1}{2}[f'''] (h-\alpha )^2\\&\quad +\frac{1}{3!}[f^{(4)}] (h-\alpha )^3+O(h^4)\\&=D_{j}^++4D_{j+1}^++D_{j+2}^++C^L_{j+1}+O(h^4).\\ \end{aligned} \end{aligned}$$

(22)

By replacing the expressions (19) in the right hand side of (13), and (22) in the left hand side, we recover the expression in (15):

$$\begin{aligned} \begin{aligned} D_{j-1}^-+4D_j^-+D_{j+1}^-&=6\frac{f^-_{j+1}-f^-_{j-1}}{2h}+C_j+O(h^4),\\ D_{j}^++4D_{j+1}^++D_{j+2}^+&=6\frac{f^+_{j+2}-f^+_{j}}{2h}+C_{j+1}+O(h^4),\\ \end{aligned} \end{aligned}$$

(23)

with the local truncation errors written as

$$\begin{aligned} \begin{aligned} C_j&=C^R_j-C^L_j=3\,{\frac{[f]}{h}}-\frac{1}{24}\,{\frac{ \left( -48\,h+72\,\alpha \right) [f']}{h}}-\frac{1}{24}\,{\frac{ \left( -36\,{\alpha }^{2}+48\, \alpha \,h-12\,{h}^{2} \right) [f'']}{h}}\\&\quad -\frac{1}{24}\,{\frac{ \left( 12 \,{\alpha }^{3}+12\,{h}^{2}\alpha -24\,{\alpha }^{2}h \right) [f'''] }{h}}-\frac{1}{24}\,{\frac{ \left( -6\,{h}^{2}{\alpha }^{2}-3\,{\alpha }^{4}+{h}^{4}+8\,{\alpha }^{3}h \right) [f^{(4)}]}{h}}\\&\quad +O(h^4),\\ C_{j+1}&=C^R_{j+1}-C^L_{j+1}=3\,{\frac{[f]}{h}}+\frac{1}{24}\,{\frac{ \left( -72\,\alpha +24\,h \right) [f']}{h}}+\frac{1}{24}\,{\frac{ \left( -24\,\alpha \,h+36\,{ \alpha }^{2} \right) [f'']}{h}}\\&\quad +\frac{1}{24}\,{\frac{ \left( 12\,{\alpha } ^{2}h-12\,{\alpha }^{3} \right) [f''']}{h}}+\frac{1}{24}\,{\frac{ \left( 3\,{\alpha }^{4}-4\,{\alpha }^{3}h \right) [f^{(4)}]}{h}}+O(h^4). \end{aligned} \end{aligned}$$

(24)

$\square$

Note that, if $x^*$ is unknown, we only need to locate the interval that contains the singularity to obtain an accurate computation of the first order derivatives at the nodes through the spline. Then, we can give a rough approximation $\tilde{x}^*$ of $x^*$ inside that interval, and use it to obtain accurate approximations of the jump in the function and its derivatives at the chosen $\tilde{x}^*$ (for example, using one-sided interpolation, as explained in [31]). For the cubic Hermite spline, we need $O(h^4)$ accuracy for [f], $O(h^3)$ accuracy for $[f']$, and so on. This observation is straightforward if we look at the expression of the correction terms in (24). Of course, the use of an inaccurate $\tilde{x}^*$ instead of $x^*$ (but still inside the interval that contains the singularity), will produce a large error in the approximation of the function in the interval that contains the singularity (typically an error of order O(1) for jumps in the function and O(h) for jumps in the first order derivative), but not in the approximation of the first order derivatives at the nodes. Moreover, it would also provide a reconstruction of the function without oscillations. Later on (in Remark 4) we will discuss how the approximation of the location of the singularity affects the interpolation of the function through the spline. The location of the interval that contains a jump discontinuity in the function, or the approximated location of a discontinuity in the first order derivative can be done using Harten’s ENO-SR (essentially non oscillatory-subcell resolution). A nice discussion about the process, with an improved algorithm for the location, is given in [32].

Let us consider now the Hermite spline in the interval [a, b], constructed using $m+1$ points belonging to a piecewise continuous function that contains a singularity at $x^*\in [x_j,x_{j+1}]$ and that is at least four times continuously differentiable except at $x^*$. Let us follow the same notation as before and denote the information to the left of the singularity with the − symbol and to the right with the $+$ symbol. The system of equations for the first order derivatives at the nodes can be expressed as

$$\begin{aligned} { \begin{aligned}&\left[ \begin{array}{lllllllllllll} 4&{}1&{}0&{}\cdots &{}0&{}0&{}0&{}0&{}\cdots &{}0&{}0\\ 1&{}4&{}1&{}\cdots &{}0&{}0&{}0&{}0&{}\cdots &{}0&{}0\\ \cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots \\ 0&{}0&{}0&{}\cdots &{}1&{}4&{}1&{}0&{}\cdots &{}0&{}0\\ 0&{}0&{}0&{}\cdots &{}0&{}1&{}4&{}1&{}\cdots &{}0&{}0\\ \cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots \\ 0&{}0&{}0&{}\cdots &{}0&{}0&{}0&{}0&{}\cdots &{}1&{}4 \end{array}\right] \left[ \begin{array}{l} D_1^-\\ D_2^-\\ \cdots \\ D_j^-\\ D_{j+1}^+\\ \cdots \\ D_{m-1}^+ \end{array}\right] =6\left[ \begin{array}{l} \frac{\delta _1^-+ \delta _2^-}{2}-\frac{D_0^-}{6}\\ \frac{\delta _2^-+\delta _3^-}{2}\\ \cdots \\ \frac{f^+_{j+1}-f^-_{j-1}}{2h}\\ \frac{f^+_{j+2}-f^-_{j}}{2h}\\ \cdots \\ \frac{\delta _{m-1}^++ \delta _{m}^+}{2}-\frac{D_m^+}{6} \end{array}\right] . \end{aligned} } \end{aligned}$$

(25)

Now we can prove the following theorem.

Theorem 2

(Accuracy of adapted first order derivatives close to singularities)

Let us consider a piecewise continuous function f that contains a singularity at $x^*$ and that is at least four times continuously differentiable on ${\mathbb {R}}\backslash \{x^*\}$ with uniformly bounded derivatives. The first order derivatives obtained through the non-corrected system of the spline for this function satisfy

$$\begin{aligned} ||f'-D||_{L^{\infty }}\le \max \{|C_j|,|C_{j+1}|\}, \end{aligned}$$

(26)

where $C_j, C_{j+1}$ are given in Lemma 1.

The addition of the correction terms $-C_j, -C_{j+1}$ up to $O(h^3)$, i.e. the subtraction of the local truncation errors in (24) up to $O(h^3)$, to the right hand side of the equations of the spline for $D_j$ and $D_{j+1}$ in the system (25) allows for the computation of first order derivatives that satisfy

$$\begin{aligned} ||f'-D||_{L^{\infty }}\le Ch^3 \sup _{c\in {\mathbb {R}}\backslash \{x^*\}}|f^{(4)}(c)|+Kh^3[f^{(4)}], \end{aligned}$$

(27)

for all $h>0$, with $K, C>0$ independent of f.

Proof

For simplicity we will take a uniform partition of the considered interval. The proof presented in what follows can be extended to a non-uniform partition. Let us represent the system in (25) by $AD=d$. Let F and r be the vectors,

$$\begin{aligned} F=\left[ \begin{array}{l} f'(x_1)\\ f'(x_2)\\ \cdots \\ f'(x_{m-1})\\ \end{array}\right] , \quad r=d-A F= A(D-F), \end{aligned}$$

(28)

where $f'(x)$ are the first order derivatives of f at the nodes.

For $i=1, \ldots , j-1,j+2,\ldots ,m-1$ we can use Taylor’s expansion to express $y_{i-1}=f(x_{i-1})$ and $y_{i+1}=f(x_{i+1})$ in terms of $f(x_i)$ and the derivatives of f at $x_i$,

$$\begin{aligned} r_i=6\left( \frac{y_{i+1}-y_{i-1}}{2h}\right) -D_{i-1}-4D_i-D_{i+1}. \end{aligned}$$

(29)

Now we can use Taylor’s expansion to express $y_{i-1}=f(x_{i-1})$, $y_{i+1}=f(x_{i+1})$, $D_{i-1}$ and $D_{i+1}$ in terms of $f(x_i)$ and the derivatives of f at $x_i$,

$$\begin{aligned} r_i= & {} \frac{3}{h}\Bigg (\left( f(x_i)+h f'(x_i)+\frac{h^2}{2}f''(x_i)+\frac{h^3}{6}f'''(x_i)+\frac{h^4}{24}f^{(4)}(\tau _1)\right) \nonumber \\{} & {} -\left( f(x_i)-hf'(x_i)+\frac{h^2}{2}f''(x_i)-\frac{h^3}{6}f'''(x_i)+\frac{h^4}{24}f^{(4)}(\tau _2)\right) \Bigg )\nonumber \\{} & {} -\left( f'(x_i)+h f''(x_i)+\frac{h^2}{2}f'''(x_i)+\frac{h^3}{6}f^{(4)}(\tau _3)\right) -4f'(x_i)\nonumber \\{} & {} -\left( f'(x_i)-hf''(x_i)+\frac{h^2}{2}f'''(x_i)-\frac{h^3}{6}f^{(4)}(\tau _4)\right) \nonumber \\= & {} \frac{1}{24}h^3\left( 3f^{(4)}(\tau _1)-3f^{(4)}(\tau _2)-4f^{(4)}(\tau _3)+4f^{(4)}(\tau _4)\right) , \end{aligned}$$

(30)

with $\tau _n\in [x_{i-1}, x_{i+1}], n=1,\ldots ,4$. Therefore, for $i=1, \ldots , j-1,j+2,\ldots ,m-1$ we obtain

$$\begin{aligned} |r_i|\le \frac{7}{12}\sup _{c\in [x_{i-1},x_{i+1}]}|f^{(4)}(c)|h^3. \end{aligned}$$

(31)

Note that if the boundary conditions $D_0^-, D_m^+$ are exact, the previous bound holds.

If we assume that the singularity is placed in the interval $[x_j, x_{j+1}]$, just as shown in Fig. 2, the process to obtain $r_j$ is similar, but we need to follow the same steps taken to obtain the local truncation error in (24): use Taylor expansions around $x^*$ to express the values $y^+_{j+1}$ and $f_x^+(x_{j+1})$ in terms of the − side (or viceversa) using the jump relations as in (18) and (21). After some basic algebra, we obtain that the $r_j$ that comes from (25) satisfies

$$\begin{aligned} \begin{aligned} r_j&=6\left( \frac{y^+_{j+1}-y^-_{j-1}}{2h}\right) -f_x^-(x_{j-1})-4f_x^-(x_j)-f_x^+(x_{j+1})\\&=6\left( \frac{y^-_{j+1}-y^-_{j-1}}{2h}\right) +C_j-f_x^-(x_{j-1})-4f_x^-(x_j)-f_x^-(x_{j+1})+O(h^4).\\ \end{aligned} \end{aligned}$$

By replicating the process for $r_{j+1}$, we easily obtain the bound in (26), corresponding to the non-corrected spline. Thus, if we subtract to $r_j$ the local truncation error $C_j$ in (24) up to $O(h^3)$, we obtain

$$\begin{aligned} { \begin{aligned} r_j-C_j&=6\left( \frac{y^-_{j+1}-y^-_{j-1}}{2h}\right) -f_x^-(x_{j-1})-4f_x^-(x_j)-f_x^-(x_{j+1})\\&\quad -\frac{1}{24}\,{\frac{ \left( -6\,{h}^{2}{\alpha }^{2}-3\,{\alpha }^{4}+{h }^{4}+8\,{\alpha }^{3}h \right) [f^{(4)}]}{h}}+O(h^4)\\&=\frac{3}{h}\left( y^-_{j+1}-y^-_{j-1}\right) -f_x^-(x_{j-1})-4f_x^-(x_j)-f_x^-(x_{j+1})\\&\quad -\frac{1}{24}\,{\frac{ \left( -6\,{h}^{2}{\alpha }^{2}-3\,{\alpha }^{4}+{h }^{4}+8\,{\alpha }^{3}h \right) [f^{(4)}]}{h}}+O(h^4)\\&=\frac{1}{24}h^3\left( 3f^{(4)}(\tau _1)-3f^{(4)}(\tau _2)-4f^{(4)}(\tau _3)+4f^{(4)}(\tau _4)\right) \\&\quad -\frac{1}{24}\,{\frac{ \left( -6\,{h}^{2}{\alpha }^{2}-3\,{\alpha }^{4}+{h }^{4}+8\,{\alpha }^{3}h \right) [f^{(4)}]}{h}}+O(h^4),\\ \end{aligned} } \end{aligned}$$

(32)

with $\tau _n\in [x_{i-1}, x_{i+1}], n=1,\ldots ,4$. For the equation of $D_{j+1}$ in (25), the result is equivalent. Taking into account the expression of A in (25), the norm of its inverse is $||A^{-1}||_{L^{\infty }}\le 1$ (see, for example, Theorem 2.1 of [37]), and since $r=A(D-F)$, then $||D-F||_{L^{\infty }}\le ||r||_{L^{\infty }}$. Therefore we obtain the bound in (27) by just applying the triangular inequality to the norm of the vector r.

Thus, the subtraction of the local truncation error to the right hand side of the system in (25),

$$\begin{aligned} { \begin{aligned}&\left[ \begin{array}{lllllllllllll} 4&{}1&{}0&{}\cdots &{}0&{}0&{}0&{}0&{}\cdots &{}0&{}0\\ 1&{}4&{}1&{}\cdots &{}0&{}0&{}0&{}0&{}\cdots &{}0&{}0\\ \cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots \\ 0&{}0&{}0&{}\cdots &{}1&{}4&{}1&{}0&{}\cdots &{}0&{}0\\ 0&{}0&{}0&{}\cdots &{}0&{}1&{}4&{}1&{}\cdots &{}0&{}0\\ \cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots &{}\cdots \\ 0&{}0&{}0&{}\cdots &{}0&{}0&{}0&{}0&{}\cdots &{}1&{}4\\ \end{array}\right] \left[ \begin{array}{l} D_1^-\\ D_2^-\\ \cdots \\ D_j^-\\ D_{j+1}^+\\ \cdots \\ D_{m-1}^+ \end{array}\right] \\&=6\left[ \begin{array}{l} \frac{\delta _1^-+ \delta _2^-}{2}-\frac{D_0^-}{6}\\ \frac{\delta _2^-+\delta _3^-}{2}\\ \cdots \\ \frac{f^+_{j+1}-f^-_{j-1}}{2h}-C_j\\ \frac{f^+_{j+2}-f^-_{j}}{2h}-C_{j+1}\\ \cdots \\ \frac{\delta _{m-1}^++ \delta _{m}^+}{2}-\frac{D_m^+}{6}\\ \end{array}\right] =6\left[ \begin{array}{l} \frac{\delta _1^-+ \delta _2^-}{2}-\frac{D_0^-}{6}\\ \frac{\delta _2^-+\delta _3^-}{2}\\ \cdots \\ \frac{f^+_{j+1}-f^-_{j-1}}{2h}\\ \frac{f^+_{j+2}-f^-_{j}}{2h}\\ \cdots \\ \frac{\delta _{m-1}^++ \delta _{m}^+}{2}-\frac{D_m^+}{6}\\ \end{array}\right] +6\left[ \begin{array}{l} 0\\ 0\\ 0\\ \cdots \\ -C_j\\ -C_{j+1}\\ \cdots \\ 0\\ 0\\ \end{array}\right] , \end{aligned} } \end{aligned}$$

(33)

allows for the obtention of first order derivatives with the bound for the error given in (27).

It is also important to note that, in order to solve the system in (33), we need to know its right hand side. The right hand side of the system in (33) is composed of finite difference approximations of the first order derivative. Thus, if the positions of the singularities are not known, the processing of these differences can be used to approximate their location. See for example the algorithm described in [32]. $\square$

Remark 2

If the data presents more than one singularity, the treatment is exactly the same. If an algorithm for locating the singularity is used, the only requirement is that the singularities are placed far enough from each other, i.e. at such a distance that the location algorithm is capable of distinguishing between two contiguous singularities. In the case of the algorithm presented in [32], singularities must be placed one stencil away from each other (four points for cubic polynomials).

4 Analysis of the Absence of the Gibbs Phenomenon Close to Jump Discontinuities in the Cubic Hermite Spline Using Corrected First Order Derivatives

In this section we will analyse the absence of the Gibbs phenomenon when using the corrected first order derivatives.

In Theorem 2 we have proved that the accuracy attained in the first order derivatives computed using the corrected system (33) is $O(h^3)$ in the infinity norm. With this result, it is posible to prove that, using these derivatives, the Hermite spline does not introduce the Gibbs phenomenon close to jump discontinuities.

We should start by introducing the definition of the Gibbs phenomenon given by D. Gottlieb and C.-W. Shu in [15].

Definition 1

The Gibbs phenomenon

Given a punctually discontinuous function f and its sampling, defined by $f(x_i)=f(ih)$, the Gibbs phenomenon deals with the convergence of the function g based on $f(x_i)$ towards f when h goes to zero. It can be characterized by two features [15]:

• 1. Away from the discontinuity, the convergence is rather slow and for any point x

  $$\begin{aligned} |f(x)-g(x)|=O(h). \end{aligned}$$

• 2. There is an overshoot, close to the discontinuity, that does not diminish with the reduction of h. Thus,

  $$\begin{aligned} \max |f(x)-g(x)| \text { does not tend to zero with } h. \end{aligned}$$

With the results about the accuracy of the first order derivatives obtained in Theorem 2 and the previous definition of the Gibbs phenomenon, we can state the following theorem.

Theorem 3

Absence of the Gibbs phenomenon

The Hermite spline obtained using the adapted first order derivatives computed through (33) does not introduce Gibbs oscillations close to jump discontinuities in the function.

Proof

We can analyse two cases. The first one is when the considered interval does not contain a singularity. In this case, Theorem 2 assures $O(h^3)$ accuracy for the first order derivatives. And, as mentioned in Remark 1, this is enough for the Hermite spline to attain $O(h^4)$ accuracy at smooth zones, so no Gibbs oscillations are possible.

The second case is when the considered interval contains a jump discontinuity in the function. Let us assume that the singularity is in the interval $[x_j, x_{j+1}]$. In (27), Theorem 2 assures that the first order derivatives attained through the corrected system in (33) also have $O(h^3)$ order of accuracy. Thus, it can be considered that the vector of first order derivatives D that results from solving (33) will be

$$\begin{aligned} ||{D}||_{L^\infty }=O(1). \end{aligned}$$

(34)

Now, considering (3) we can write

$$\begin{aligned} \begin{aligned} a_j&=\frac{hD_{j+1}+D_j h+2y_j-2y_{j+1}}{h^3}=\frac{O(h)+O(h) +2y_j-2y_{j+1}}{h^3},\\ b_j&=-\frac{hD_{j+1}+2D_j h+3y_j-3y_{j+1}}{h^2}=-\frac{O(h)+2O(h)+3y_j-3y_{j+1}}{h^2},\\ c_j&=D_j=O(1),\\ d_j&=y_j. \end{aligned} \end{aligned}$$

(35)

If we consider that there is a jump discontinuity in the interval $[x_j, x_{j+1}]$, then $y_{j+1}-y_j=O(1)$ and

$$\begin{aligned} \begin{aligned} d_j&=y_j,\\ c_j&=D_j=O(1),\\ b_j&=-\frac{hD_{j+1}+2D_j h+3y_j-3y_{j+1}}{h^2}=-\frac{O(h)+2O(h) -3O(1)}{h^2}=\frac{O(1)}{h^2}=O\left( \frac{1}{h^2}\right) ,\\ a_j&=\frac{hD_{j+1}+D_j h+2y_j-2y_{j+1}}{h^3}=\frac{O(h)+O(h) -2O(1)}{h^3}=\frac{O(1)}{h^3}=O\left( \frac{1}{h^3}\right) . \end{aligned} \end{aligned}$$

(36)

And if $x\in [x_j, x_{j+1}]$, the equation of the spline (1) transforms into

$$\begin{aligned} g_j(x)= & {} {a}_j(x-x_j)^3+{b}_j(x-x_j)^2+{c}_j(x-x_j)+{d}_j=O\left( \frac{1}{h^3}\right) O\left( h^3\right) \nonumber \\{} & {} +O\left( \frac{1}{h^2}\right) O\left( h^2\right) +O\left( h\right) +y_j=y_j+O(1). \end{aligned}$$

(37)

This means that the error attained by the Hermite spline with $O(h^3)$ accurate first order derivatives provides $O(h^4)$ accuracy when interpolating the data, except at the interval that contains the discontinuity, where the approximation is O(1) accurate.

Now we only need to check that the approximation provided by the spline is in the interval $[y_j, y_{j+1}]$ when h goes to zero. In order to verify this assumption, we can write the equation of the spline in (1) as

$$\begin{aligned} { \begin{aligned} g_j(x)=&{a}_j(x-x_j)^3+{b}_j(x-x_j)^2+{c}_j(x-x_j)+{d}_j\\ =&\frac{hD_{j+1}+D_j h+2y_j-2y_{j+1}}{h^3}(x-x_j)^3\\&-\frac{hD_{j+1}+2D_j h+3y_j-3y_{j+1}}{h^2}(x-x_j)^2\\&+D_j(x-x_j)+y_j. \end{aligned} } \end{aligned}$$

(38)

Introducing the change of variables $s=\frac{x-x_j}{h}$, we get

$$\begin{aligned} { \begin{aligned} g_j(x)&=\left( hD_{j+1}+D_j h+2y_j-2y_{j+1}\right) s^3\\&\quad -\left( hD_{j+1}+2D_j h+3y_j-3y_{j+1}\right) s^2+D_j h s+y_j,\\&=\left( 1-s^2 (3-2 s)\right) y_j+s^2 (3-2 s) y_{j+1}+\left( s^3-2s^2+s\right) hD_j+\left( s^3-s^2\right) hD_{j+1}\\&=b_1(s) y_{j+1}+b_2(s) y_j+b_3(s) hD_j+b_4(s) hD_{j+1}, \end{aligned} } \end{aligned}$$

(39)

for $s\in [0, 1]$. The expression can also be reformulated as follows:

$$\begin{aligned} \begin{aligned} g_j(x)&=y_j+\left( s^2 (3-2 s)\right) (y_{j+1}-y_j)+\left( s^3-2s^2+s\right) hD_j+\left( s^3-s^2\right) hD_{j+1}, \end{aligned} \end{aligned}$$

(40)

where the first two terms of (40) are just a dilation and a translation of the Hermite base function $b_1(s)=s^2 (3-2 s)$. This function presents a minimum at $s=0$ and a maximum at $s=1$, so no Gibbs oscillations can be introduced by this element of the Hermite base, as can be seen in Fig. 1. In this Figure it can be observed that the base functions $b_3(s)$ and $b_4(s)$ oscillate. Thus, Gibbs oscillations can appear if these bases appear multiplied by large coefficients. As we have already mentioned, in (40) $D_j$ and $D_{j+1}$ are O(1). With these considerations in mind, it is not difficult to conclude that the last two terms in (40) go to zero when $h\rightarrow 0$. So no Gibbs oscillation is possible.

It should be noted that an analogous analysis can be done for the classical spline. In this case, the vector of first order derivatives D that results from solving (33) without the correction terms will be

$$\begin{aligned} ||{D}||_{L^\infty }=O\left( \frac{1}{h}\right) , \end{aligned}$$

(41)

due to the bound in (26) and the expressions for $C_j, C_{j+1}$ given in Lemma 1. Following exactly the same arguments as before and rewriting the expressions from (35)–(40) taking into account (41) in the process, we can easily conclude that the classical spline introduces the Gibbs phenomenon close to jump discontinuities in the function. $\square$

5 Full Accuracy of the Adapted Hermite Interpolation

Once the first order derivatives have been calculated with $O(h^3)$ accuracy, if we use them directly in the expression in (4) without any other enhancement, the interpolation will lose its accuracy at the interval $[x_j, x_{j+1}]$ that contains the singularity, as it has been analysed in Theorem 3. This is logical, as the interpolation in the interval $[x_j, x_{j+1}]$ can be seen as a particularization of a linear combination of the polynomials of the Hermite basis, which are continuous. This situation can be solved by defining a piecewise continuous Hermite polynomial at the interval $[x_j, x_{j+1}]$.

Remark 3

The exact location of jump discontinuities and the definition of the function at $x^*$ is lost during the point-values discretization. Thus, to obtain a reconstruction at $x^*$, not only the location of the discontinuity must be provided, but also which of the two definitions corresponds to this point. If the discontinuities are in the derivatives, this problem does not arise, as the function is continuous.

Let us consider the following lemma.

Lemma 2

Let us assume that the first order derivatives of the function at the nodes have been obtained with $O(h^3)$ accuracy. Let us assume that the data presents a singularity at $x^*$ with $x_{j}\le x^*\le x_{j+1}$ and, without loss of generality, that we know that $f(x^*)$ belongs to the − side of the function. Then, the local truncation error of the non-corrected spline in the interval $[x_{j}, x^*]$ is

$$\begin{aligned} { \begin{aligned} C^-(x)&= \left( 3\,{\frac{ \left( x-x_j \right) ^{2}}{{h}^{2}}}-2\,{ \frac{ \left( x-x_j \right) ^{3}}{{h}^{3}}} \right) [f]\\&\quad + \left( -{\frac{ \left( -2\,h+3\,\alpha \right) \left( x-x_j \right) ^{2}}{{h}^{2}}}+{\frac{ \left( -h+2\,\alpha \right) \left( x-x_j \right) ^{3}}{{h}^{3}}} \right) [f']\\&\quad + \left( {\frac{ \left( - \left( h-\alpha \right) ^{2}+h \left( h-\alpha \right) \right) \left( x-x_j \right) ^{3}}{{h}^{3}}}-{\frac{ \left( -\frac{3}{2}\, \left( h-\alpha \right) ^{2}+h \left( h-\alpha \right) \right) \left( x-x_j \right) ^{2}}{{h}^{2}}} \right) [f'']\\&\quad + \left( {\frac{ \left( \frac{1}{2}\,h \left( h-\alpha \right) ^{2}-\frac{1}{3}\, \left( h- \alpha \right) ^{3} \right) \left( x-x_j \right) ^{3}}{{h}^{3}}} -{\frac{ \left( \frac{1}{2}\,h \left( h-\alpha \right) ^{2}-\frac{1}{2}\, \left( h- \alpha \right) ^{3} \right) \left( x-x_j \right) ^{2}}{{h}^{2}}} \right) \\&\qquad [f''']+O(h^4).\\ \end{aligned}} \end{aligned}$$

(42)

In the interval $(x^*, x_{j+1}]$ the local truncation error is

$$\begin{aligned} { \begin{aligned} C^+(x)&= C^-(x)-\left( [f]+(x-x^*)[f']+\frac{1}{2}(x-x^*)^2[f'']+\frac{1}{3!}(x-x^*)^3[f''']\right) \\&\quad +O(h^4). \end{aligned} } \end{aligned}$$

(43)

Proof

As mentioned in the lemma, we have two possible cases:

• The first one is when the interpolation is being obtained at a position $x_j\le x\le x^*$. In this case, the first order derivative $D^+_{j+1}$ and the data point $y^+_{j+1}$ belong to the $+$ side of the singularity while the interpolation point x is at the − side. We can use the same technique applied before to obtain the correction terms for the first order derivatives in order to assure the accuracy of the interpolation. Let us consider again Fig. 2. In the interval $[x_j, x_{j+1}]$, the Hermite interpolation in (4) will have the following expression:

  $$\begin{aligned} \begin{aligned} g_j(x)=&\,\frac{hD^+_{j+1}+D^-_j h+2y^-_j-2y^+_{j+1}}{h^3}(x-x_j)^3\\&-\frac{hD^+_{j+1}+2D^-_j h+3y^-_j-3y^+_{j+1}}{h^2}(x-x_j)^2 +D^-_j(x-x_j)+y^-_j. \end{aligned} \end{aligned}$$

  (44)

  Now we can use (18) and (21) to obtain an expression of the corrected spline in the interval $[x_j, x^*]$:

  $$\begin{aligned} g_j(x)&=\frac{hD^-_{j+1}+D^-_j h+2y^-_j-2y^-_{j+1}}{h^3}(x-x_j)^3\nonumber \\&\quad -\frac{hD^-_{j+1}+2D^-_j h+3y^-_j-3y^-_{j+1}}{h^2}(x-x_j)^2 +D^-_j(x-x_j)+y^-_j\nonumber \\&\quad + \left( 3\,{\frac{ \left( x-x_j \right) ^{2}}{{h}^{2}}}-2\,{ \frac{ \left( x-x_j \right) ^{3}}{{h}^{3}}} \right) [f]\nonumber \\&\quad + \left( -{\frac{ \left( -2\,h+3\,\alpha \right) \left( x-x_j \right) ^{2}}{{h}^{2}}}+{\frac{ \left( -h+2\,\alpha \right) \left( x-x_j \right) ^{3}}{{h}^{3}}} \right) [f']\nonumber \\&\quad + \left( {\frac{ \left( - \left( h-\alpha \right) ^{2}+h \left( h-\alpha \right) \right) \left( x-x_j \right) ^{3}}{{h}^{3}}}-{\frac{ \left( -\frac{3}{2}\, \left( h-\alpha \right) ^{2}+h \left( h-\alpha \right) \right) \left( x-x_j \right) ^{2}}{{h}^{2}}} \right) [f'']\nonumber \\&\quad + \left( {\frac{ \left( \frac{1}{2}\,h \left( h-\alpha \right) ^{2}-\frac{1}{3}\, \left( h- \alpha \right) ^{3} \right) \left( x-x_j \right) ^{3}}{{h}^{3}}} -{\frac{ \left( \frac{1}{2}\,h \left( h-\alpha \right) ^{2}-\frac{1}{2}\, \left( h- \alpha \right) ^{3} \right) \left( x-x_j \right) ^{2}}{{h}^{2}}} \right) \nonumber \\&\quad \quad [f''']+O(h^4)\nonumber \\&= g^-_j(x)+C^-(x)+O(h^4). \end{aligned}$$

  (45)

  In fact, the $O(h^4)$ term in the previous expression can be written as

  $$\begin{aligned} T_4^-(x)= & {} \left( {\frac{ \left( \frac{1}{6}\,h \left( h-\alpha \right) ^{3}-\frac{1}{12}\, \left( h-\alpha \right) ^{4} \right) \left( x-{ xj} \right) ^{3}}{{h}^{3}}}-{\frac{ \left( \frac{1}{6}\,h \left( h-\alpha \right) ^{3}-\frac{1}{8}\, \left( h-\alpha \right) ^{4} \right) \left( x-{ xj} \right) ^{2}}{{h}^{2}}} \right) \nonumber \\{} & {} \quad [f^{(4)}]+\cdots \end{aligned}$$

  (46)

• The second one is when the interpolation is being obtained at a position $x^*< x\le x_{j+1}$. In this case, the first order derivative $D^-_{j}$ and the data $y^-_{j}$ belong to the − side of the singularity while the interpolation point x is at the $+$ side. Proceeding as before, we have

  $$\begin{aligned} { \begin{aligned} g_j(x)=&\,\frac{hD^+_{j+1}+D^+_j h+2y^+_j-2y^+_{j+1}}{h^3}(x-x_j)^3\\&-\frac{hD^+_{j+1}+2D^+_j h+3y^+_j-3y^+_{j+1}}{h^2}(x-x_j)^2 +D^+_j(x-x_j)+y^+_j\\&+\left( -2\,{\frac{ \left( x-x_j \right) ^{3}}{{h}^{3}}}-1+3\,{ \frac{ \left( x-x_j \right) ^{2}}{{h}^{2}}} \right) [f]\\&+ \left( {\frac{ \left( -h+2\,\alpha \right) \left( x-x_j \right) ^{3}}{{h}^{3}}}-{\frac{ \left( -2\,h+3\,\alpha \right) \left( x-x_j \right) ^{2}}{{h}^{2}}}+\alpha -x+x_j \right) [f']\\&+ \left( {\frac{ \left( h\alpha -{\alpha }^{2} \right) \left( x-x_j \right) ^{3}}{{h}^{3}}}-{\frac{ \left( 2\,h\alpha -\frac{3}{2}\,{ \alpha }^{2} \right) \left( x-x_j \right) ^{2}}{{h}^{2}}}+\alpha \, \left( x-x_j \right) -\frac{1}{2}\,{\alpha }^{2} \right) [f'']\\&+ \left( {\frac{ \left( \frac{1}{3}\,{\alpha }^{3}-\frac{1}{2}\,h{\alpha }^{2} \right) \left( x-x_j \right) ^{3}}{{h}^{3}}}+\frac{1}{6}\,{\alpha }^{3}-\frac{1}{2}\,{ \alpha }^{2} \left( x-x_j \right) -{\frac{ \left( \frac{1}{2}\,{\alpha }^{ 3}-h{\alpha }^{2} \right) \left( x-x_j \right) ^{2}}{{h}^{2}}} \right) \\&\quad [f''']+O(h^4)\\ =&\, g^+_j(x)+C^+(x)+O(h^4)\\ =&\, g^+_j(x)+C^-(x)-\left( [f]+(x-x^*)[f']+\frac{1}{2}(x-x^*)^2[f'']+\frac{1}{3!}(x-x^*)^3[f''']\right) \\&\quad +O(h^4). \end{aligned} } \end{aligned}$$

  (47)

  In the previous expression, the $O(h^4)$ term can be written as

  $$\begin{aligned} T_4^+(x)= & {} \left( {\frac{ \left( -\frac{1}{12}\,{\alpha }^{4}+\frac{1}{6}\,h{\alpha }^{3} \right) \left( x-{ xj} \right) ^{3}}{{h}^{3}}}-{\frac{ \left( -\frac{1}{8}\,{\alpha }^{4}+\frac{1}{3}\,h{\alpha }^{3} \right) \left( x-{ xj} \right) ^{2}}{{h}^{2}}}-\frac{1}{24}\,{\alpha }^{4}+\frac{1}{6}\,{\alpha }^{3} \left( x -{ xj} \right) \right) \nonumber \\{} & {} \quad [f^{(4)}]+\cdots \end{aligned}$$

  (48)

  From (47), we can see that we obtain again a relation between the correction terms for the spline similar to that obtained for the data in (18),

  $$\begin{aligned} C^+(x)= & {} C^-(x)-\left( [f]+(x-x^*)[f']+\frac{1}{2}(x-x^*)^2[f'']+\frac{1}{3!}(x-x^*)^3[f''']\right) \\{} & {} +O(h^4). \end{aligned}$$

$\square$

Now we have all the necessary tools to prove the following theorem.

Theorem 4

(Accuracy of the corrected Hermite interpolation close to singularities)

Let us consider a piecewise continuous function f that contains a singularity at $x^*\in [x_j,x_{j+1}]$ and that is at least four times continuously differentiable on ${\mathbb {R}}\backslash \{x^*\}$. Let us assume, without loss of generality, that we know that $f(x^*)$ belongs to the − side of the function. Let us also consider that the first order derivatives at the nodes are available with at least $O(h^3)$ accuracy. The approximation obtained through the non-corrected Hermite spline for this function satisfies

$$\begin{aligned} |f(x)-g_j(x)|\le \max \{|C^+(x)|,|C^-(x)|\}, x\in [x_{j},x_{j+1}], \end{aligned}$$

(49)

where $C^+(x), C^-(x)$ are given in Lemma 2.

The addition of the correction terms $-C^-(x)$ when $x\in [x_j,x^*]$ and $-C^+(x)$ when $x\in (x^*,x_{j+1}]$, i.e. the subtraction of the local truncation error given in Lemma 2 to the approximation obtained by the non-corrected Hermite spline results in the accuracy

$$\begin{aligned} |f(x)-g_j(x)|\le \left\{ \begin{array}{ll} \displaystyle \max \left\{ \frac{13}{24}h^4 |[f^{(4)}]|, \frac{1}{384} h^4\sup _{\xi \in [x_{j},x^*]}|f^{(4)}(\xi )|\right\} , &{} x\in [x_{j},x^*],\\ \displaystyle \max \left\{ \frac{2}{3}h^4 |[f^{(4)}]|, \frac{1}{384} h^4\sup _{\xi \in (x^*,x_{j+1}]}|f^{(4)}(\xi )|\right\} , &{} x\in (x^*,x_{j+1}],\\ \displaystyle \frac{1}{384} h^4\sup _{\xi \in [x_j, x_{j+1}]}|f^{(4)}(\xi )|, &{} \text {in any other case,}\\ \end{array} \right. \end{aligned}$$

(50)

for all $h>0$, independent of f.

Proof

For smooth functions, the error has been analysed in Theorem 1 and is given in (8).

Now we can consider the case when there is a singularity in the interval $[x_j,x_{j+1}]$. Looking at the results obtained in expressions (42) and (43) of Lemma 2, we can see the local truncation error obtained when approximating the data using the uncorrected spline. We can just change the sign of the local truncation error, i.e isolate $g^-_j(x)$ or $g^+_j(x)$, to write that

$$\begin{aligned} g^-_j(x)= & {} g_j(x)-C^-(x)+O(h^4),\quad x_j<x\le x^*, \end{aligned}$$

(51)

$$\begin{aligned} g^+_j(x)= & {} g_j(x)-C^+(x)+O(h^4),\quad x^*<x<x_{j+1}, \end{aligned}$$

(52)

and to conclude that the addition of the correction terms $-C^-(x)$ in the interval $[x_j,x^*]$ or $-C^+(x)$ in the interval $(x^*,x_{j+1}]$ allows us to keep the accuracy of the spline close to the singularities. In fact, looking at the $O(h^4)$ terms in the local truncation error, we can see that the term in (46) can be roughly bounded by

$$\begin{aligned} |T_4^-(x)|\le \frac{13}{24} |[f^{(4)}]|h^4, \end{aligned}$$

(53)

and the one in (48) can be bounded in the same way by

$$\begin{aligned} |T_4^+(x)|\le \frac{2}{3} |[f^{(4)}]|h^4. \end{aligned}$$

(54)

$\square$

Remark 4

Let us assume now that we do not know the exact location of the singularity $x^*$ and we approximate it by $\tilde{x}^*$. For simplicity, we can assume that $\tilde{x}^*$ is located to the right of $x^*$ at a distance $\beta$ (which is the error of location). Then, using Taylor’s expansion on the jump in the function or the derivatives, we get

$$\begin{aligned} \begin{aligned} \left[ f(\tilde{x}^*)\right]&=[f]+\beta [f']+\frac{\beta ^2}{2}[f'']+\cdots ,\\ \left[ f'(\tilde{x}^*)\right]&=[f']+\beta [f'']+\frac{\beta ^2}{2}[f''']+\cdots \end{aligned} \end{aligned}$$

(55)

Therefore, from replacing these expressions in the local truncation error in (42) or (43), we derive that an error in the location of the discontinuity leads to an error of the size of the jump of the function in the interval between $\tilde{x}^*$ and ${x}^*$. In the case of jumps in the first order derivative, the problem can be solved computing an accurate enough $\tilde{x}^*$. Thus, in order to keep order of accuracy $O(h^4)$ for the cubic spline, we just need the error of location $\beta$ to be of order $O(h^{4})$. It is noteworthy that if there is a false detection of the singularity at a smooth zone, the jumps should be zero (if known a priori) or close enough to zero (if approximated) [31].

Remark 5

Although we have introduced the corrected cubic Hermite spline for cases where the location of the singularities and the jump in the function and its derivatives are known, this information can be approximated numerically in certain instances.

When working with data derived from the discretization of piecewise smooth functions, sometimes it is possible to detect the singularities and compute their location. Some kinds of singularities are jump discontinuities in the function and kinks, meaning jumps in the function or the first order derivative respectively. The location of these singularities can be found depending on the discretization of the data used. Kinks can be located using the point-value discretization, i.e. a sampling of the function. There is no way of locating the position of jump discontinuities using this kind of discretization: the exact position is lost during the discretization process. For data discretized using point-values, the location of kinks can be obtained, for example, using the classical approach proposed by Harten for his ENO subcell resolution (ENO-SR) algorithm [32, 33]. In [32], the authors introduce the ENO-SR algorithm and propose its implementation for the point-values discretization of the data.

Once the location of the singularity has been approximated, it is also possible to approximate the jump in the function and its derivatives using one-sided interpolation. See for example [31] for a detailed explanation about this point.

We will verify the previous assertions in the Numerical experiments section.

Corollary 2

(Error of the second order derivative of the corrected cubic Hermite interpolation)

Under the assumptions of Theorem (4), the error for the second order derivative of the corrected cubic Hermite interpolation in the interval that contains the singularity is given by

$$\begin{aligned} |f''(x)-g''_j(x)|\le \left\{ \begin{array}{ll} \displaystyle \max \left\{ \frac{1}{3}h^2 |[f^{(4)}]|, \frac{1}{2} h^2\sup _{\xi \in [x_{j},x^*]}|f^{(4)}(\xi )|\right\} , &{} x\in [x_{j},x^*],\\ \displaystyle \max \left\{ \frac{5}{12}h^2 |[f^{(4)}]|, \frac{1}{2} h^2\sup _{\xi \in (x^*,x_{j+1}]}|f^{(4)}(\xi )|\right\} , &{} x\in (x^*,x_{j+1}],\\ \displaystyle \frac{1}{2} h^2\sup _{\xi \in [x_j,x_{j+1}]}|f^{(4)}(\xi )|, &{} \text {in any other case,}\\ \end{array} \right. \end{aligned}$$

(56)

for all $h>0$, independent of f.

Proof

If we differentiate the error terms in (46) and (48) of Lemma (2), it is not difficult to arrive to the expression for the error of the second order derivative of the corrected spline in the intervals $[x_j,x^*)$ or $(x^*,x_{j+1}]$:

$$\begin{aligned} \begin{aligned} (T_4^-(x))''&=\left( -\frac{(\alpha - h)^3 (-3 \alpha - h + 6 \alpha h^2 (x - x_j) + 6 h^3 (x - x_j))}{12 h^2}\right) [f^{(4)}]+\cdots ,\\ (T_4^+(x))''&=\left( \frac{3 \alpha ^4 (h -2(x-xj))-4 h\alpha ^3(2h-3 (x - x_j))}{12 h^3}\right) [f^{(4)}]+\cdots \end{aligned} \end{aligned}$$

(57)

Taking into account these results and Corollary 1, we finish the proof:

$$\begin{aligned} |f''(x)-g''_j(x)|\le \left\{ \begin{array}{ll} \displaystyle \max \left\{ \frac{1}{3}h^2 |[f^{(4)}]|, \frac{1}{2} h^2\sup _{\xi \in [x_{j},x^*]}|f^{(4)}(\xi )|\right\} , &{} x\in [x_{j},x^*],\\ \displaystyle \max \left\{ \frac{5}{12}h^2 |[f^{(4)}]|, \frac{1}{2} h^2\sup _{\xi \in (x^*,x_{j+1}]}|f^{(4)}(\xi )|\right\} , &{} x\in (x^*,x_{j+1}],\\ \displaystyle \frac{1}{2} h^2\sup _{\xi \in [x_j,x_{j+1}]}|f^{(4)}(\xi )|, &{} \text {in any other case,}\\ \end{array} \right. \end{aligned}$$

(58)

$\square$

Corollary 3

(Smoothness of the classical Hermite spline with corrected first order derivatives)

The correction proposed in Theorem 4 leads to a piecewise $C^2$ spline that preserves the singularity of the function at $x^*$. Consequently, the use of the adapted first order derivatives obtained in Sect. 3 without the correction proposed in Theorem 4 produces a jump in the second order derivative of the resultant spline $[S''(x)]=S''(x^+)-S''(x^-)$, at $x_j$ and $x_{j+1}$ if the singularity is placed in the interval $[x_j, x_{j+1}]$. This jump is equal to

$$\begin{aligned} {[}S''(x_j)]= & {} \frac{6 }{h^2}[f] + \frac{(-6 \alpha + 4 h) }{h^2}[f'] + \frac{(3 \alpha ^2 +4 \alpha h) }{h^2}[f''] + \frac{(-\alpha ^3 + 2 \alpha ^2 h) }{h^2}[f'''],\nonumber \\ \end{aligned}$$

(59)

$$\begin{aligned} {[}S''(x_{j+1})]= & {} -\frac{6 }{h^2}[f] + \frac{(6 \alpha - 2 h) }{h^2}[f'] + \frac{(-3 \alpha ^2 + 10 \alpha h) }{h^2}[f''] + \frac{(\alpha ^3 - \alpha ^2 h) }{h^2}[f'''].\nonumber \\ \end{aligned}$$

(60)

Thus, in this case, the resultant spline is $C^1$. More specifically, it is $C^2$ except at $x_j$ and $x_{j+1}$, if $[x_j, x_{j+1}]$ is the interval that contains the singularity.

Proof

Let us consider that we apply the correction terms proposed in Theorem 4. In this case, the modification of the right hand side of the system in (33) does not affect the regularity conditions imposed at all the nodes except at $x_j$ and $x_{j+1}$. This is very easy to verify by computing the jumps $[S''(x_{i})]$ and $[S'(x_{i})]$ for $x_i, i=1, \ldots , j-1,j+2,\ldots m-1$, using the expression of the spline in (4), and confirming that these jumps are zero. Thus, we only have to consider the jump in the first and second order derivatives at these points. By differentiating once or twice the expression in (4) for the spline at $x_j^-$ and the expression in (51) for the spline at $x_j^+$, we obtain $[S'(x_j)]=S'(x_j^+)-S'(x_j^-)=0$ and $[S''(x_j)]=S''(x_j^+)-S''(x_j^-)=0$. In the same way, differentiating once or twice (4) for the spline at $x_{j+1}^+$ and (52) for the spline at $x_{j+1}^-$ we obtain $[S''(x_{j+1})]=S''(x_{j+1}^+)-S''(x_{j+1}^-)=0$ and $[S'(x_{j+1})]=S'(x_{j+1}^+)-S'(x_{j+1}^-)=0$. This means that the addition of the correction terms shown in Theorem 4 assures piecewise $C^2$ regularity of the spline everywhere except at the singularity placed at $x^*$.

By repeating the same calculations but eliminating the correction term $C^-(x)$ in (51) or $C^+(x)$ in (52) we obtain that the jumps in the first order derivative are still zero at $x_j$ and $x_{j+1}$, and we obtain the jumps in the second order derivative shown in (59) and (60) respectively. $\square$

The conclusions reached in Corollary 3 will be verified in the Numerical experiments section (Subsect. 7.3).

6 The Adaption as a Post-processing Procedure

In this section we will analyse the possibility of correcting the spline as a post-processing. We have two possibilities: the first one is when we want to obtain a smooth function in the interval that contains the singularity; the second one is when we want to obtain a piecewise continuous function in that interval. In the first case we only have to incorporate the correction of the first order derivatives exposed in Sect. 3. In the second case we need to include also the correction shown in Sect. 5. Let us start by the first case.

6.1 The Cubic Hermite Spline with Corrected First Order Derivatives as a post-processing procedure

Looking at the right hand side of the system in (33), we can see that the correction of the first order derivatives of the spline can be done as a post-processing, just adding to the solution of the uncorrected spline the solution of the system considering only as right hand side the sparse vector that contains the correction terms. It is clear that the system can be expressed as

$$\begin{aligned} A\bar{D}=d+C, \end{aligned}$$

(61)

and it is clear that the term that will add the desired accuracy to the vector is the solution of the system

$$\begin{aligned} A\tilde{D}=C. \end{aligned}$$

(62)

The solution $\tilde{D}$ of this system will be a vector which, added to the vector of first order derivatives obtained from the classical system of the spline

$$\begin{aligned} AD=d, \end{aligned}$$

(63)

will provide $O(h^3)$ accuracy for the first order derivatives even close to the singularities. Looking at the expression of the spline in (38), if we express each particular first order derivative resulting from (61)

$$\begin{aligned} \bar{D}=D+\tilde{D}, \end{aligned}$$

then the adapted spline can be expressed as

$$\begin{aligned} \begin{aligned} \bar{g}_i(x)=&\, \left( \frac{hD_{i+1}+D_i h+2y_i-2y_{i+1}}{h^3}(x-x_i)^3-\frac{hD_{i+1}+2D_i h+3y_i-3y_{i+1}}{h^2}(x-x_i)^2 \right. \\&\left. +D_i(x-x_i)+y_i\right) +\frac{h\tilde{D}_{i+1}+\tilde{D}_i h}{h^3}(x-x_i)^3-\frac{h\tilde{D}_{i+1}+2\tilde{D}_i h}{h^2}(x-x_i)^2\\&+\tilde{D}_i(x-x_i)=g_i(x)+\tilde{g}_i(x). \end{aligned} \end{aligned}$$

(64)

The term $\tilde{g}_i(x)$ can be interpreted as a spline where the function values are zero and the first order derivatives are different from zero. It can be added as a post-processing of the classical Hermite spline in order to attain $O(h^3)$ order of accuracy for the first order derivatives at the nodes.

In [34] p. 22, the authors provide a very efficient algorithm for solving tridiagonal systems with dominant main diagonal, like the ones in (25) and (33). If we have the equations

$$\begin{aligned} \begin{aligned} b_1x_1+c_1x_2&=d_1,\\ a_2x_1+b_2x_2+c_2x_3&=d_2,\\ \cdots&\\ a_{m-1}x_{m-2}+b_{m-1}x_{m-1}+c_{m-1}x_m&=d_{m-1},\\ a_{m}x_{m-1}+b_{m}x_{m}&=d_{m},\\ \end{aligned} \end{aligned}$$

(65)

then we can compute for $j=1,2, \ldots , m$,

$$\begin{aligned} \begin{aligned} p_j&=a_jq_{j-1}+b_j, \quad (q_0=0),\\ q_j&=-c_j/p_j,\\ u_j&=(d_j-a_ju_{j-1})/p_j,\quad (u_0=0). \end{aligned} \end{aligned}$$

The successive elimination of $x_1, x_2, \ldots , x_{m-1}$ from the second, third, $\ldots , m$th equations yields the equivalent system of equations

$$\begin{aligned} \begin{aligned} x_j&=q_{j}x_{j+1}+u_j, \quad (j=1, \ldots ,m-1),\\ x_m&=u_m, \end{aligned} \end{aligned}$$

where $x_m, x_{m-1}, \ldots , x_1$ are evaluated successively. We note that the quantities $p_j$ and $q_j$ depend on the mesh but not on the ordinates at the mesh nodes (that only appear in the $d_j$ in (65)). Thus, the correction can be calculated very efficiently once the $p_j$’s and $q_j$’s of the classical spline have been computed, as we are only changing the right hand side of the system of equations of the classical spline, i.e. the $d_j$’s.

6.2 The Fully Accurate Cubic Hermite Spline as a Post-processing Procedure

The fully accurate spline presents differences with the one presented in the previous subsection only at the interval $[x_j, x_{j+1}]$ that contains the singularity. Thus, it is possible to express the corrected spline with corrected first order derivatives inside that interval $g_j(x)$ as

$$\begin{aligned} g^-_j(x)= & {} g_j(x)+\tilde{g}_j(x)-C^-(x)+O(h^4),\quad x_j<x<x^*, \end{aligned}$$

(66)

$$\begin{aligned} g^+_j(x)= & {} g_j(x)+\tilde{g}_j(x)-C^+(x)+O(h^4),\quad x^*<x<x_{j+1}. \end{aligned}$$

(67)

7 Numerical Experiments

In this section we analyse the accuracy obtained when computing the corrected first order derivatives proposed in Sect. 3 and in the interpolation obtained through the corrected Hermite spline described in Sect. 5. We will compare the results of the corrected Hermite spline with those obtained by the classical non-corrected spline.

7.1 Accuracy Analysis of the Interpolation, First and Second Order Derivatives Obtained Through the Classical and Corrected Splines

Let us try to analyse the accuracy of the interpolation, first and second order derivatives attained by the corrected and non-corrected spline. We will consider that we start from discretized data that comes from the sampling of a function. The location of the discontinuity will be considered known in the case of jump discontinuities in the function. When the jumps are in the first order derivative, we will obtain approximations of its location using the algorithm in [32] (See, for example, the plot to the right in Fig. 3).

In order to check the accuracy of the spline between the nodes or the first and second order derivatives at the nodes, we can perform a grid refinement analysis. For the approximation of the derivatives, we can just obtain analytically the corresponding derivative of the function and compare it with the approximation of the derivative obtained through the spline for each step of the refinement analysis. The process for this grid refinement analysis will be the following:

• 1. Sample the function at a high resolution with the mesh size $h_0$ and keep this data to check the error.

• 2. Subsample the high resolution data using a mesh size $h_1$ to obtain the nodes of the spline.

• 3. Obtain the analytical first and second order derivatives of the function at the nodes obtained in the previous step.

• 4. Compute the spline with the low resolution nodes.

• 5. Compute the errors for the first order derivative at the low resolution nodes and for the function and the second order derivative at the high resolution data.

• 6. Reduce the mesh size for the high resolution data and between nodes so that $h_0=\frac{h_0}{2}, h_1=\frac{h_1}{2}$ and go back to step 1.

Once we have the errors for the data, the first, and the second order derivative for all the refinement steps, we can just apply the classical formula for approximating the numerical accuracy,

$$\begin{aligned} O^l=\log _2\left( \frac{E^l}{E^{l+1}}\right) , \end{aligned}$$

using the infinity norm $E^l=||f^l-\tilde{f^l}||_{L^\infty }$, where we have denoted the grid refinement step with the super index l.

7.1.1 Experiment 1: A Function with a Jump Discontinuity

For this first experiment we have chosen the function in (68), which corresponds to a function with jumps in the function and the first four derivatives:

$$\begin{aligned} f(x)=\left\{ \begin{array}{ll} (x - \frac{\pi }{6})(x - \frac{\pi }{6} - 3)^3 + 8\sin \left( \frac{\pi x}{8}\right) + 10, &{} \text { if } 0\le x<\frac{\pi }{6},\\ 8\sin \left( \frac{\pi x}{8}\right) , &{} \text { if } \frac{\pi }{6}\le x\le 1. \end{array} \right. \end{aligned}$$

(68)

Fig. 3

Functions in (68) and (69) used in the numerical experiments

In this case, we assume that we know the location of the discontinuity, but we approximate the jumps in the function and the derivatives using one-sided interpolation [31].

In this experiment we have set the relation between the high resolution data and the nodes to $\frac{h_1}{h_0}=32$, which means that we take one point of every 32 points of the high resolution data to obtain the low resolution data, i.e. the nodes of the spline. Thus, the low resolution nodes have been sampled with $m=2^l$ nodes and the high resolution data with 32m points. The results are shown in Table 1 for the interpolation of the spline at the high resolution data obtained from the low resolution data, in Table 2 for the first order derivative at the nodes and in Table 3 for the second order derivative at the high resolution data points. We can see that the accuracy obtained by the corrected spline between the nodes tends to $O(h^4)$ when h goes to zero and to O(1) for the classical spline. The accuracy of the first order derivative at the nodes tends to $O(h^3)$ for the corrected spline and to O(1/h) for the classical one. The accuracy of the second order derivative at the high resolution data tends to $O(h^2)$ for the corrected spline and to $O(1/h^2)$ for the classical one. These results have been also presented in Fig. 4, in a semilogarithmic scale in order to show the numerical decreasing of the error and to compare it with the theoretical one, which has been represented with dashed blue lines. In the plot to the left of this figure we can see the errors obtained by the classical spline (stars in red) and by the corrected spline (stars in blue) at the high resolution nodes. The slope of the dashed lines represents the decreasing of the error of the classical spline with order of accuracy O(1) and with $O(h^4)$ for the corrected spline. At the center we present the errors in the first order derivatives at the low resolution nodes (resulting from the solution of the system of equations of the spline) for each spline with the same markers. The slope of the dashed lines represent again the decreasing of the error for the classical spline with $O\left( \frac{1}{h}\right)$ accuracy and with $O(h^3)$ for the corrected spline. The plot to the right shows the errors obtained in the second order derivative at the high resolution nodes. In this case, the slope of the dashed lines represents the decreasing of the error for the classical spline with $O\left( \frac{1}{h^2}\right)$ accuracy and with $O(h^2)$ for the corrected spline. In Subsect. 7.2 we show that the loss of accuracy of the classical approach is due, among other things, to Gibbs-like oscillations close to the discontinuity.

Table 1 Grid refinement analysis for the accuracy of the corrected and classical splines using the infinity norm at the high resolution nodes

Table 2 Grid refinement analysis for the accuracy of the first order derivative at the low resolution nodes of the corrected and classical splines using the infinity norm

Table 3 Grid refinement analysis for the accuracy of the second order derivatives at the high resolution nodes of the corrected and classical splines using the infinity norm

Fig. 4

In this figure we represent the results of the grid refinement analysis for the accuracy of the corrected and classical splines and their first and second order derivatives using the infinity norm shown in Tables 1, 2 and 3. The original data has been sampled from the function in (68). In the first row we show the error $E^l$ for the corrected and classical spline, in the second one we show the error for the first order derivative, and finally, the error for the second order derivative. The x axis represents de resolution level l: the increase of l implies the reduction of the grid-spacing of the low resolution data. Thus, the slope of the dashed lines represents the accuracy of each method

7.1.2 Experiment 2: A Function with a Discontinuity in the Derivatives

For the second experiment we have chosen the function in (69)

$$\begin{aligned} { f(x)=\left| \sin \left( x-\frac{\pi }{6}\right) \right| ,\quad 0\le x<1, } \end{aligned}$$

(69)

which corresponds to a function with jumps in the odd derivatives. In this case, we approximate the location of the singularity with $O(h^4)$ accuracy using the algorithm in [32], and we also approximate the jumps in the function (which in this case is zero) and the derivatives [31]. Following the same setting as in the previous experiment, the results are shown in Table 4 for the interpolation of the spline, in Table 5 for the first order derivatives at the nodes and in Table 6 for the second order derivative of the spline at the high resolution data points. It can be observed that with the corrected spline we are capable of obtaining high resolution approximations of the function, the first and the second order derivatives, while with the classical spline the approximations are affected by the presence of the singularity. These results have been presented in Fig. 5, where the slope of the dashed blue lines represent again the theoretical decreasing of the error for each case. In this figure we represent the results of the grid refinement analysis for the accuracy of the corrected and classical splines and their first and second order derivatives using the infinity norm shown in Tables 4, 5 and 6.

Table 4 Grid refinement analysis for the accuracy of the corrected and classical splines using the infinity norm at the high resolution nodes

Table 5 Grid refinement analysis for the accuracy of the first order derivative at the nodes of the corrected and classical splines using the infinity norm

Table 6 Grid refinement analysis for the accuracy of the second order derivatives at the high resolution nodes of the corrected and classical splines using the infinity norm

Fig. 5

In this figure we represent the results of the grid refinement analysis for the accuracy of the corrected and classical splines and their first and second order derivatives using the infinity norm shown in Tables 4, 5 and 6. The original data has been sampled from the function in (69). As in Fig. 4, in the first row we show the error $E^l$ for the corrected and classical spline, in the second one we show the error for the first order derivative, and finally, the error for the second order derivative. The x axis represents de resolution level l: the increase of l implies the reduction of the grid-spacing of the low resolution data. Thus, the slope of the dashed lines represents the accuracy of each method

7.2 Elimination of the Gibbs Phenomenon Using Adapted First Order Derivatives

In this section we will analyse the behaviour of the corrected spline, the classical spline with corrected first order derivatives and the classical spline close to discontinuities. In order to do so, we will use the function in (68), which presents a jump discontinuity. As in Subsect. 7.1.1, we assume that we know the location of the discontinuity for this function, but we approximate the jumps in the function and the derivatives using one-sided interpolation [31]. The results of the experiments are shown in Fig. 6. In the first row, we present the approximation (red circles) of the function in (68) sampled with a high resolution of 256 points (blue crosses) and then subsampled to obtain 16 points (solid black circles) in order to obtain the low resolution data. We can see that the classical spline (first row, plot to the left) shows Gibbs oscillations close to the discontinuity. If we correct the first order derivatives of the spline and we do not use any correction for the spline itself (first row, central plot), we obtain an approximation that is free of oscillations close to the discontinuity, but that presents diffusion in the interval that contains the discontinuity. If we use the corrected first order derivatives and we also correct the approximation of the spline close to the discontinuity (first row, plot to the right), we can obtain a piecewise continuous function without oscillations close to the discontinuity. The second row of Fig. 6 shows the error obtained in the computation of the first order derivatives at the low resolution nodes for each spline: to the left we show the result of the classical spline and at the center and to the right, the error for corrected first order derivatives (which are the same both for the classical spline with corrected first order derivatives and for the corrected spline). The third row shows the error at the high resolution nodes attained by each of the splines in the same order. From these plots we can see that the use of corrected first order derivatives leads to the elimination of the Gibbs oscillations close to the discontinuities. If we also correct the interpolation of the spline, it is possible to obtain a piecewise continuous function with high accuracy even close to the discontinuities.

Fig. 6

In the first row, we present the reconstruction or approximation g(x) of the function in (68) with red circles. Initially we sample the function at a high resolution of 256 points and we represent it with blue crosses. Then we subsample the data to obtain a low resolution version of 16 points, and we represent it with solid black circles. The second and third row are dedicated to present the error in the first and second order derivatives, respectively. The first column presents the results of the classical cubic Hermite spline, the second one presents those obtained by the classical cubic Hermite spline with corrected first order derivatives. Finally, the third one presents the results attained by the corrected cubic Hermite spline with corrected first order derivatives (Color figure online)

7.3 Regularity

In this subsection we will analyse the behaviour of the first and second order derivatives between the nodes for the classical spline, the classical spline with corrected first order derivatives and the corrected spline with corrected first order derivatives. As before, for the function in (68), we assume that we know the location of the discontinuity, and, for the function in (69), we approximate the location of the singularity with $O(h^4)$ [32]. In both cases we approximate the jumps in the function and the derivatives using one-sided interpolation [31]. Figures 7 and 8 show the results obtained for the functions in (68) and (69) respectively. In both figures, the approximation of the first order derivative of the functions is shown in the first column, and the approximation for the second order derivatives is shown in the second column. Functions in (68) and (69) are sampled with 16 points to obtain the low resolution data. Then, each polynomial piece of the spline is differentiated once or twice to obtain the results shown in Figs. 7 and 8. We have represented some points of the first or the second order derivative of the function with blue dots. The stars represent a sampling of the first or the second order derivatives of each one of the polynomial pieces of the spline g(x). Each polynomial piece is represented with a different colour. Our objective is just to compare qualitatively the first and second order derivatives of the functions, and the approximations obtained through each spline. In each of the two figures, the first row presents the result for the classical spline, and the second row shows the result for the classical spline with corrected first order derivatives. Finally, the third row shows the results for the corrected spline with corrected first order derivatives. It is clear that the classical spline introduces oscillations in the first and the second order derivatives close to the singularities, but maintains the $C^2$ regularity in the whole domain. The results also show that the modification of the right hand side of the system of equations of the spline in (33) eliminates the oscillations in the function outside the interval that contains the singularity, but leads to the loss of the regularity in the second order derivative if the spline is not corrected, as discussed in Corollary 3. Only the corrected spline with corrected first order derivatives, discussed in Sect. 5, allows for the elimination of the oscillations and keeps the correct piecewise regularity of the function from which the data has been obtained.

Fig. 7

Approximation of the first order derivative (first column), and second order derivative (second column) of the function in (68) sampled with 16 points to obtain the low resolution data. The blue dots represent the first or the second order derivative of the function. The coloured stars represent the first or the second order derivatives of the polynomial pieces of the spline g(x). Each colour of the stars represents a different polynomial piece of the spline. The first row presents the result for the classical spline. The second row shows the result for the classical spline with corrected first order derivatives. The third row shows the results for the corrected spline with corrected first order derivatives (Color figure online)

Fig. 8

Approximation of the first order derivative (first column) and second order derivative (second column) of the function in (69) sampled with 16 points to obtain the low resolution data. The blue dots represent the first or the second order derivative of the function. The coloured stars represent the first or the second order derivatives of the polynomial pieces of the spline g(x). Each colour of the stars represents a different polynomial piece of the spline. The first row presents the result for the classical spline. The second row shows the result for the classical spline with corrected first order derivatives. The third row shows the results for the corrected spline with corrected first order derivatives (Color figure online)

8 Conclusion

In this article we have presented a new technique that allows for the correction of the first order derivatives and the approximation obtained through cubic Hermite splines close to singularities. The technique consists in taking into account the effect introduced by the singularity by adding correction terms, in order to assure a high order of accuracy close to the singularities. The correction can be included in the right hand side of the system of the spline, allowing for its application as a post-processing. It has been shown theoretically that the correction of the first order derivatives eliminates the Gibbs phenomenon introduced by the classical spline close to jump discontinuities. Through the correction introduced in the system of equations of the spline, we obtain accurate first order derivatives even close to the singularities. Thus we have also discussed the theoretical accuracy in the infinity norm obtained through the correction of the first order derivatives. Even so, diffusion is always present in the interval that contains a jump discontinuity if we only correct the first order derivatives of the spline. If the correction is also introduced in the equation of the spline (1) in the interval that contains the singularity, it is possible to approximate accurately piecewise continuous functions. Proofs for the accuracy attained through the corrected spline in the infinity norm have been provided. The new technique requires the knowledge of the location of the singularities and the jumps in the function and its derivatives at the singularities, or high order of accuracy approximations. We have also discussed the order of accuracy needed in the approximated location of the singularity to keep the order of the spline. The numerical experiments confirm the theoretical results obtained.