Numerical Interpolation and Applications

Charbel Tannous

► To cite this version:

Charbel Tannous. Numerical Interpolation and Applications. Master. Analyse Numérique, UBO Brest, France. 2016, pp.23. hal-04961026

HAL Id: hal-04961026

https://hal.science/hal-04961026v1

Submitted on 21 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical Interpolation and Applications

C. Tannous

<figure><img src="image_2.png" /></figure>

Université de Brest, Lab-STICC, CNRS-UMR 6285, F-29200 Brest, FRANCE

Interpolation is used after signal sampling or after collecting series of measurements during a physics, engineering experiment or some statistical polling. Interpolation transforms data into a continuous function that can be treated with standard analytical or transform methods. Applications cover Digital Signal Processing, CAD (Computer Assisted Design), Computer Graphics, Font Design... Analytic interpolation made by Planck is also described in order to distinguish it from its numerical version. It allowed to discover Quantum Mechanics and introduce Bose-Einstein distribution in Statistical Physics.

PACS numbers: 02.60.Ed, 02.10.Ud, 02.60.Cb, 03.65.-w

Keywords: Interpolation methods, Linear Algebra, Numerical simulation, solution of equations, Quantum mechanics

Version: February 21, 2025

## Contents

<table><tr><td>I. Introduction</td><td>1</td></tr><tr><td>II. Lagrange Interpolation</td><td>2</td></tr><tr><td>III. Hermite Interpolation</td><td>3</td></tr><tr><td>IV. Newton Interpolation</td><td>3</td></tr><tr><td>A. Linear Interpolation</td><td>5</td></tr><tr><td>B. Quadratic Interpolation</td><td>6</td></tr><tr><td>V. Splines</td><td>7</td></tr><tr><td>A. Cubic Spline Interpolation</td><td>7</td></tr><tr><td>B. Natural Spline</td><td>8</td></tr><tr><td>C. Bézier Curves</td><td>9</td></tr><tr><td>D. B-Spline Interpolation</td><td>10</td></tr><tr><td>E. Hermite splines</td><td>13</td></tr><tr><td>VI. Applications</td><td>15</td></tr><tr><td>A. Vector Graphics in font design and Typography</td><td>16</td></tr><tr><td>B. Digital signal processing (DSP)</td><td>17</td></tr><tr><td>1. Reconstruction of Signal by Interpolation</td><td>18</td></tr><tr><td>2. Digital-to-analog conversion (D/A)</td><td>19</td></tr><tr><td>C. Computer Graphics and Surface interpolation</td><td>20</td></tr><tr><td>A. Analytic interpolation and Planck radiation law</td><td>21</td></tr><tr><td>References</td><td>23</td></tr></table>

## I. INTRODUCTION

When an experiment is performed and data points are collected or when a signal is sampled providing data of the form $(x_i, f(x_i))$ for $i = 0, 1, 2, ..., n$ where $f$ is a function representing the outcome or the experiment or signal, the issue is based on the following assumption: It is possible to use these data points to represent $f(x)$ at $x \neq x_i$ for $i = 0, 1, ..., n$. Some methods to address this procedure are based on using a piecewise polynomial or other representation of $f$ as explained further below.

However there are three possibilities for representing $f$: Interpolation, Approximation or Extrapolation. When interpolation is considered, the values $(x_i, f(x_i))$, $i = 1...n$ are considered an exact representation of $f$ whereas in
approximation, one considers that these values contains errors regarding $f$ representation [1, 2]. Finally, supposing that $p$ is a (piecewise) polynomial representation of $f$, if $x$ is inside the smallest interval containing the $x_i, i = 1...n$, then $p(x)$ is a polynomial interpolation of $f(x)$, otherwise $p(x)$ is a polynomial extrapolation of $f(x)$.

Numerical interpolation is different from Analytic interpolation as was done by Planck to reconcile the high-temperature and low-temperature behavior of Black body radiation (cf. Appendix).

## II. LAGRANGE INTERPOLATION

Given a number of nodes $x_i, i = 1...n$ it is possible to build a polynomial function going through all these distinct points, with Lagrange interpolation [1, 3, 4] formula:

$$L_j(x) = \prod_{\substack{i=1 \\ i \neq j}}^{n} \frac{(x - x_i)}{(x_j - x_i)} \quad (1)$$

such that $p_j(x_i) = \delta_{ij}$.

More precisely, if $f: [a, b] \to \mathbf{R}$ is a function and $a = x_0 < x_1 < ... < x_n = b$, then the polynomial $p$ of degree $n$ defined by

$$p(x) = \sum_{i=0}^{n} f(x_i) \prod_{\substack{j=0 \\ j \neq i}}^{n} \left( \frac{x - x_j}{x_i - x_j} \right) \quad (2)$$

is such that $p(x_i) = f(x_i)$ for $0 \le i \le n$. The polynomial $p$ is called the Lagrange Interpolating Polynomial of $f$ at $x_0, x_1, ..., x_n$.

The polynomial $p$ is often used as an approximation of $f$ on the interval $[x_0, x_n]$.

Since polynomials of degree $n$ have exactly $n$ complex roots counted with multiplicity, the Lagrange Interpolating Polynomial in (2) is the unique polynomial of degree at most $n$ satisfying $p(x_i) = f(x_i)$ for $0 \le i \le n$.

### Runge phenomenon

Runge phenomenon is a problem of oscillation (cf. Fig. 1) at the edges of an interval that occurs when using polynomial interpolation with polynomials of high degree over a set of periodic interpolation points. It is similar to Gibbs phenomenon when Fourier harmonics are summed at a function discontinuity.

<figure><img src="image_3.png" /><figcaption>**Fig.1:** The black curve is the exact Runge function $f(x) = \frac{1}{1+25x^2}$. Third, 5-th and 9-th-order interpolating polynomials are shown. Between some interpolating points (close to points 1 and -1), the error between the Runge function and the interpolating polynomial increases with higher order [4].</figcaption></figure>

Consider the case where one desires to interpolate through $n + 1$ equally spaced points of a function $f(x)$ using the $n$-degree polynomial $P_n(x)$ that passes through those points. Naturally, one might think that using more points would lead to a more accurate reconstruction of $f(x)$. However, this particular set of polynomial functions $P_n(x)$ is not guaranteed to have the property of uniform convergence.
The $P_n(x)$ produced in this manner may in fact diverge away from $f(x)$ as $n$ increases; this occurs in an oscillating pattern that is amplified near interpolation end points.

Consider the Runge function $f(x) = \frac{1}{1 + 25x^2}$ interpolated at equidistant points $x_i \in [-1, 1]$ such that: $x_i = \frac{2i}{n} - 1$, $i \in \{0, 1, \dots, n\}$ with a polynomial $P_n(x)$ of degree $\le n$, the resulting interpolation oscillates toward the end of the interval, i.e. close to -1 and 1 (cf. Fig. 1).

## Derivative estimation:

Let us reconstruct a continuous and noise-free function $y = f(x)$ that passes through a set of $n + 1$ sample points $x_i, y_i = f(x_i)$, around a point $x$ for which the derivative $f'(x)$ needs to be estimated by Lagrange interpolation:

$$p_n(x) = \sum_{i=-n/2}^{n/2} \prod_{j=-n/2(\ne i)}^{n/2} \frac{(x - x_j)}{(x_i - x_j)} y_i$$

$p_n(x)$ is a polynomial function of $x$ of order $n$ which passes through all $n + 1$ sample points $p_n(x_i) = y_i$, $i = -n/2, \dots, n/2$.

For example, when $n = 2$:

$$p_2(x) = \frac{(x - x_0)(x - x_1)}{(x_1 - x_0)(x_1 - x_1)} y_1 + \frac{(x - x_{-1})(x - x_1)}{(x_0 - x_{-1})(x_0 - x_1)} y_0 + \frac{(x - x_{-1})(x - x_0)}{(x_1 - x_{-1})(x_1 - x_0)} y_1$$

The derivative of function $f(x)$ at any point $x$ can be estimated by differentiating $P_n(x)$. For example, when $n = 2$:

$$p'_{2}(x) = \frac{2x - x_{0} - x_{1}}{(x_{1} - x_{0})(x_{1} - x_{1})} y_{1} + \frac{2x - x_{-1} - x_{1}}{(x_{0} - x_{-1})(x_{0} - x_{1})} y_{0} + \frac{2x - x_{-1} - x_{0}}{(x_{1} - x_{-1})(x_{1} - x_{0})} y_{1}$$

In particular, when the point $x = x_0$ is the middle point of the set of $n + 1$ points:

$$p'_{2}(x_{0}) = p'_{2}(x)|_{x=x_{0}} = \frac{x_{0} - x_{1}}{(x_{1} - x_{0})(x_{1} - x_{1})} y_{1} + \frac{(x_{0} - x_{-1}) + (x_{0} - x_{1})}{(x_{0} - x_{-1})(x_{0} - x_{1})} y_{0} + \frac{x_{0} - x_{-1}}{(x_{1} - x_{-1})(x_{1} - x_{0})} y_{1}$$

## III. HERMITE INTERPOLATION

Hermite method is an osculatory interpolation [3] requiring the function and its derivatives $f_i^{(k)}$ at the given points $x_i, i = 0, \dots, n$ with $a = x_0 < x_1 < \dots < x_n = b$ and $k = 0, \dots, n-1$ such that the interpolating polynomial writes: $p(x) = \sum_{i=0}^{n} \sum_{k=0}^{n_i-1} f_i^{(k)} L_{ik}(x)$ where $L_{ik}(x)$ are the generalized Lagrange polynomials and defined from the auxiliary polynomials:

$$\ell_{ik}(x) = \frac{(x - x_i)^k}{k!} \prod_{\substack{j=0 \\ j \ne i}}^{n} \frac{(x - x_j)}{x_i - x_j}, \quad 0 \le i \le n, \ 0 \le k \le n_i \qquad (3)$$

to get recursively, starting from $L_{i,n_i-1}(x) = \ell_{i,n_i-1}, i = 0, \dots, n$, and: $L_{ik}(x) = \ell_{ik}(x) - \sum_{\nu=k+1}^{n_i-1} \ell_{ik}^{(\nu)}(x)L_{i\nu}(x)$, $k = n_i - 2, n_i - 3, \dots, 0$ obeying the constraints:

$$p(x_i) = f_i, \quad p^{(1)}(x_i) = f_i^{(1)}, \quad p^{(2)}(x_i) = f_i^{(2)}, \dots, p^{(n)}(x_i) = f_i^{(n)}, \quad i = 0, \dots, k \qquad (4)$$

where the upper index $(i)$, $i = 0, \dots, k$ indicates derivation order.

As an example, a first order Hermite osculatory interpolation is treated in subsection VE.

## IV. NEWTON INTERPOLATION

The definition of interpolating polynomial of a function at $(n + 1)$ distinct points $x_0, x_1, \dots, x_n$ can be extended to the case where the $x_i$ are not all distinct.
Let $f : ]a, b[ \to \mathbf{R}$ and $g : ]a, b[ \to \mathbf{R}$ be two functions sufficiently differentiable. Suppose that $x_0, x_1, \dots, x_n$ are $(n+1)$ points in $]a, b[$ (not necessarily distinct). $f$ and $g$ agree at the points $x_0, x_1, \dots, x_n$ when:

$$ \frac{d^j f}{dx^j}(z) = \frac{d^j g}{dx^j}(z) $$

for $j = 0, 1, \dots, m-1$ whenever $z$ appears $m$ times in the list $x_0, x_1, \dots, x_n$. Then:

$$ \frac{d^j f}{dx^j}(z) = f(z) \quad \text{for } j = 0 $$

## Unicity theorem

Let $f : ]a, b[ \to \mathbf{R}$ be a function sufficiently differentiable. Suppose that $x_0, x_1, \dots, x_n$ are $(n+1)$ points in $]a, b[$ not necessarily distinct. Then there is a unique polynomial $p$ of degree at most $n$ such that $f$ and $p$ agree at $x_0, x_1, \dots, x_n$.

The polynomial $p$ in the Unicity theorem is called the interpolating polynomial of $f$ at the interpolatory points $x_0, x_1, \dots, x_n$.

Let $f : ]a, b[ \to \mathbf{R}$ be a function sufficiently differentiable. The $k^{th}$ divided difference of $f$ at $k+1$ not necessarily distinct points $x_0, x_1, \dots, x_k$ in $]a, b[$, denoted $f[x_0, x_1, \dots, x_k]$, is the coefficient of $x^k$ in the unique polynomial of degree at most $k$ that agrees with $f$ at $x_0, x_1, \dots, x_k$.

Let $f : ]a, b[ \to \mathbf{R}$ be a function sufficiently differentiable. Suppose that $x_0, x_1, \dots, x_n$ are $n+1$ not necessarily distinct points in $]a, b[$. Then the unique polynomial $p$ of degree at most $n$ that agrees with $f$ at $x_0, x_1, \dots, x_n$ is given by

$$ \begin{aligned} p(x) = f[x_0] &+ f[x_0, x_1](x - x_0) + f[x_0, x_1, x_2](x - x_0)(x - x_1) + \dots \\ &+ f[x_0, x_1, \dots, x_n](x - x_0)(x - x_1) \dots (x - x_{n-1}). \end{aligned} \quad (5) $$

Moreover, for $x_j \neq x_{j+k}$,

$$ f[x_j, x_{j+1}, \dots, x_{j+k}] = \frac{f[x_{j+1}, x_{j+2}, \dots, x_{j+k}] - f[x_j, x_{j+1}, \dots, x_{j+k-1}]}{x_{j+k} - x_j} \quad (6) $$

and, for $x_j = x_{j+1} = x_{j+2} = \dots = x_{j+k}$,

$$ f[x_j, x_{j+1}, \dots, x_{j+k}] = \frac{1}{k!} \frac{d^k f}{dx^k}(x_j). \quad (7) $$

Finally,

$$ f(x) = p(x) + f[x_0, x_1, \dots, x_n, x] (x - x_0)(x - x_1) \dots (x - x_{n-1})(x - x_n). \quad (8) $$

It is then easy to deduce the first divided difference of $f$.

The steps are given by:

* The interpolating polynomial $p$ of $f$ of degree 0 at $x_0$ is given by the constant function $p(x) = f(x_0)$ for all $x$. Hence, the coefficient $f[x_0]$ of $x^0$ is

$$ f[x_0] = f(x_0). $$

* The interpolating polynomial $p$ of $f$ of degree at most 1 at $x_0, x_1$ is given by

$$ p(x) = \begin{cases} f(x_0) + \frac{f(x_1) - f(x_0)}{x_1 - x_0} (x - x_0) & \text{if } x_0 \neq x_1 \\ f(x_0) + f'(x_0) (x - x_0) & \text{if } x_0 = x_1 \end{cases} $$

In the first case, it is the equation of the secant line through $(x_0, f(x_0))$ and $(x_1, f(x_1))$. In the second case, it is the equation of the tangent line at $x_0$ since $p(x_0) = f(x_0)$ and $p'(x_0) = f'(x_0)$. Thus, the coefficient $f[x_0, x_1]$ of $x^1$ is

$$ f[x_0, x_1] = \begin{cases} \frac{f(x_1) - f(x_0)}{x_1 - x_0} & \text{if } x_0 \neq x_1 \\ f'(x_0) & \text{if } x_0 = x_1 \end{cases} $$
* The interpolating polynomial $p$ of $f$ of degree at most 2 at $x_0, x_1, x_2$ is of the form

$$p(x) = A + B(x - x_0) + C(x - x_0)(x - x_1).$$

If the $x_i$ are distinct, the polynomial $p$ must satisfy

$$\begin{align*}
f(x_0) &= p(x_0) = A, \\
f(x_1) &= p(x_1) = A + B(x_1 - x_0), \\
f(x_2) &= p(x_2) = A + B(x_2 - x_0) + C(x_2 - x_0)(x_2 - x_1).
\end{align*}$$

Hence,

$$A = f(x_0) = f[x_0], \quad B = \frac{f(x_1) - f(x_0)}{x_1 - x_0} = f[x_0, x_1], \quad C = \frac{f[x_1, x_2] - f[x_0, x_1]}{x_2 - x_0}.$$

Thus, the coefficient $f[x_0, x_1, x_2]$ of $x^2$ is

$$f[x_0, x_1, x_2] = \frac{f[x_1, x_2] - f[x_0, x_1]}{x_2 - x_0}$$

if the $x_i$ are distinct. A similar formula is obtained if $x_0 = x_1 \neq x_2$ or $x_1 \neq x_2 = x_3$.

If $x_0 = x_1 = x_2$, the interpolating polynomial $p$ of $f$ of degree at most 2 is given by the Taylor polynomial of degree 2 at $x_0$; namely,

$$p(x) = f(x_0) + f'(x_0)(x - x_0) + \frac{1}{2!} f''(x_0)(x - x_0)^2.$$

It is easy to check that $p(x_0) = f(x_0)$, $f'(x_0) = p'(x_0)$ and $f''(x_0) = p''(x_0)$. Hence, the coefficient $f[x_0, x_1, x_2]$ of $x^2$ is

$$f[x_0, x_1, x_2] = \frac{1}{2!} f''(x_0) \quad \text{if} \quad x_0 = x_1 = x_2$$

1. There are two ways to represent the polynomial of degree at most $n$ that agrees with $f$ at the $n+1$ distinct points $x_0, x_1, \dots, x_n$,

2. Because the interpolating polynomial of degree at most $n$ of $f$ at $x_0, x_1, x_2, \dots, x_n$ is independent of the order in which the $x_i$ are listed, in particular the coefficient of $x^n$ is not going to change if the order of the $x_i$ is changed, then:

$$f[x_0, x_1, \dots, x_k] = f[x_{\sigma(0)}, x_{\sigma(1)}, \dots, x_{\sigma(k)}]$$

for any permutation $\sigma$ of $\{0, 1, 2, 3, \dots, k\}$.

3. To be able to use (7) and thus get simple divided difference formulae, assuming that if a value $z$ appears more than once in $\{x_0, x_1, x_2, \dots, x_n\}$, then all occurrences of $z$ are contiguous.

4. From a computational point of view, (5) is better than (2) because less operations are needed to evaluate $p(x)$ if the nested form is used to evaluate polynomials.

5. The form (5) of the interpolating polynomial of $f$ at $x_0, x_1, \dots, x_n$ can be easily extended to the form (5) of the interpolating polynomial of $f$ at the $n+1$ points $x_0, x_1, \dots, x_n$ and $x_{n+1}$, where $x_{n+1}$ is a new point. Only the divided difference $f[x_0, x_1, \dots, x_n, x_{n+1}]$ needs to be computed because we already have the divided differences $f[x_0], f[x_0, x_1], \dots, f[x_0, x_1, \dots, x_n]$ from the interpolating polynomial of $f$ at $x_0, x_1, \dots, x_n$.

## A. Linear Interpolation

Consider only two points $(x_0, f(x_0))$ and $(x_1, f(x_1))$ with $x_0 \neq x_1$. Then

$$p(x) = f[x_0] + f[x_0, x_1](x - x_0),$$
where

$$f[x_0] = f(x_0) \quad \text{and} \quad f[x_0, x_1] = \frac{f(x_1) - f(x_0)}{x_1 - x_0}.$$

If $(x_0, f(x_0)) = (2.2, 6.2)$ and $(x_1, f(x_1)) = (2.5, 6.7)$, find an approximation of $f(x)$ at $x = 2.35$.

Then $f[2.2] = f(2.2) = 6.2$ and

$$f[2.2, 2.5] = \frac{f(2.5) - f(2.2)}{2.5 - 2.2} = \frac{6.7 - 6.2}{2.5 - 2.2} = 1.66.$$

Thus $p(x) = 6.2 + 1.66(x - 2.2)$ and $f(2.35) \approx p(2.35) = 6.45$.

Suppose that only the values of $f(x) = \cos(x)$ at $x = 0$ and $x = \pi/6$ are known, find an approximation of $\cos(0.2)$.

Computing the interpolating polynomial at the points

$$(0, \cos(0)) = (0, 1) \quad \text{and} \quad (\pi/6, \cos(\pi/6)) = (\pi/6, \sqrt{3}/2).$$

Then $f[0] = f(0) = 1$ and

$$f[0, \pi/6] = \frac{f(\pi/6) - f(0)}{\pi/6 - 0} = \frac{\sqrt{3}/2 - 1}{\pi/6 - 0} = \frac{3\sqrt{3} - 6}{\pi}.$$

Thus

$$p(x) = 1 + \left( \frac{3\sqrt{3} - 6}{\pi} \right) x$$

and $\cos(0.2) \approx p(0.2) \approx 0.948825$. After six digit rounding, $\cos(0.2) = 0.980067$. Thus, the absolute error is about 0.031242.

Consider a function $f$ associated to data set

$$\{(x_0, f(x_0)), (x_1, f(x_1)), \dots, (x_n, f(x_n))\},$$

where $x_0 < x_1 < \dots < x_n$. It is possible to plot the graph of the piecewise linear function $p$ that passes through each point of the data set. Namely, $f$ is approximated by the piecewise polynomial $p$ (each piece is a polynomial of degree one) defined by

$$p(x) = f[x_i] + f[x_i, x_{i+1}] (x - x_i), \quad x_i \le x \le x_{i+1},$$

for $i = 0, 1, 2, \dots, n-1$.

## B. Quadratic Interpolation

Consider the points $(x_0, f(x_0))$, $(x_1, f(x_1))$ and $(x_2, f(x_2))$ with $x_i \neq x_j$ for $i \neq j$. Then

$$p(x) = f[x_0] + f[x_0, x_1] (x - x_0) + f[x_0, x_1, x_2] (x - x_0)(x - x_1),$$

where

$$f[x_0] = f(x_0), \quad f[x_0, x_1] = \frac{f(x_1) - f(x_0)}{x_1 - x_0}, \quad f[x_0, x_1, x_2] = \frac{f[x_1, x_2] - f[x_0, x_1]}{x_2 - x_0}.$$

Moreover, one must compute $f[x_1, x_2] = \frac{f(x_2) - f(x_1)}{x_2 - x_1}$ in the formula of $f[x_0, x_1, x_2]$.

Given the data $(x_0, f(x_0)) = (2.2, 6.2)$, $(x_1, f(x_1)) = (2.5, 6.7)$ and $(x_2, f(x_2)) = (2.7, 6.5)$. Find the approximation of $f(x)$ at $x = 2.35$.

Then $f[2.2] = f(2.2) = 6.2$,

$$f[2.2, 2.5] = \frac{f(2.5) - f(2.2)}{2.5 - 2.2} = \frac{6.7 - 6.2}{2.5 - 2.2} = 1.66,$$

$$f[2.5, 2.7] = \frac{f(2.7) - f(2.5)}{2.7 - 2.5} = \frac{6.5 - 6.7}{2.7 - 2.5} = -1.$$
and

$$f[2.2, 2.5, 2.7] = \frac{f[2.5, 2.7] - f[2.2, 2.5]}{2.7 - 2.2} = \frac{-1 - 1.66}{2.7 - 2.2} = -5.32.$$

Thus

$$\begin{aligned} p(x) &= 6.2 + 1.66(x - 2.2) - 5.32(x - 2.2)(x - 2.5) \\ &= 6.2 + (x - 2.2)(1.66 - 5.32(x - 2.5)) \end{aligned}$$

$$\text{and } f(2.35) \approx p(2.35) = 6.57.$$

# V. SPLINES

Several methods are described below to generate a piecewise [5] polynomial function $p$ that may provide a good representation of the function $f$ [6].

## A. Cubic Spline Interpolation

Let $f : [a, b] \to \mathbf{R}$ be a continuously differentiable function and $a = x_0 < x_1 < \dots < x_n = b$. Recall that one uses a piecewise linear function through the points $(x_i, f(x_i))$ for $0 \le i \le n$ to sketch the graph of $f$, i.e. the plot of a piecewise linear curve that looks like the graph of $f$. Instead of using linear interpolation to join the points $(x_{i-1}, f(x_{i-1}))$ and $(x_i, f(x_i))$, cubic polynomial interpolation is used on the intervals $[x_{i-1}, x_i]$. The function $p$ obtained is called a piecewise cubic polynomial [6] Using cubic polynomials, a better representation may be imposed between $f$ and $p$ than with linear polynomials. Cubic spline interpolation is ideal for approximating function with discontinuous derivatives.

Let $f : [a, b] \to \mathbf{R}$ be a continuously differentiable function and $a = x_0 < x_1 < \dots < x_n = b$. A free or natural spline interpolant for the function $f$ on the nodes $x_0, x_1, \dots, x_n$ is a piecewise cubic polynomial $p$ defined as follows.

$$p(x) = p_i(x) \quad \text{if} \quad x_i \le x \le x_{i+1},$$

where the $p_i$ are polynomials of degree three that satisfy

1. $p_i(x_i) = f(x_i)$ for $i = 0, 1, \dots, n-1$,

2. $p_i(x_{i+1}) = f(x_{i+1})$ for $i = 0, 1, \dots, n-1$,

3. $p'_i(x_{i+1}) = p'_{i+1}(x_{i+1})$ for $i = 0, 1, \dots, n-2$,

4. $p''_i(x_{i+1}) = p''_{i+1}(x_{i+1})$ for $i = 0, 1, \dots, n-2$,

5. $p''_0(x_0) = p''_{n-1}(x_n) = 0$.

If the third, fourth and fifth conditions are replaced by

1. $p'_i(x_i) = f'(x_i)$ for $i = 0, 1, \dots, n-1$,

2. $p'_i(x_{i+1}) = f'(x_{i+1})$ for $i = 0, 1, \dots, n-1$,

Then $p$ is called a piecewise cubic Hermite interpolant for the function $f$ on the nodes $x_0, x_1, \dots, x_n$.

It is possible to find a cubic polynomials $p_i$ needed for the spline interpolant [6]. Let $z_i = p''(x_i)$ for $i = 0, 1, \dots, n$. Values of the $z_i$ are needed to satisfy the natural cubic splines.

Assuming that the $p_i$ are cubic polynomials, $p_i''$ is a linear function through the points $(x_i, z_i)$ and $(x_{i+1}, z_{i+1})$. Recall that $\Delta x_i = x_{i+1} - x_i$. Hence,

$$p_i''(x) = \frac{z_{i+1} - z_i}{\Delta x_i} x + \frac{x_{i+1}z_i - x_iz_{i+1}}{\Delta x_i} = \left(\frac{z_i}{\Delta x_i}\right) (x_{i+1} - x) + \left(\frac{z_{i+1}}{\Delta x_i}\right) (x - x_i).$$
<figure><img src="image_4.png" /><figcaption>Fig.2: Cubic spline interpolating a series of arbitrarily spaced data points $x_i$, $f(x_i)$.</figcaption></figure>

Integrating twice gives

$$ 
\begin{aligned}
p_i(x) &= \left(\frac{z_i}{6\Delta x_i}\right) (x_{i+1} - x)^3 + \left(\frac{z_{i+1}}{6\Delta x_i}\right) (x - x_i)^3 + A_i x + B_i \\
&= \left(\frac{z_i}{6\Delta x_i}\right) (x_{i+1} - x)^3 + \left(\frac{z_{i+1}}{6\Delta x_i}\right) (x - x_i)^3 + C_i(x - x_i) + D_i(x_{i+1} - x),
\end{aligned} 
\quad (9) $$

where $C_i - D_i = A_i$ and $D_i x_{i+1} - C_i x_i = B_i$.

From $p_i(x_i) = f(x_i)$ and $p_i(x_{i+1}) = f(x_{i+1})$, then

$$ f(x_i) = \left(\frac{z_i}{6}\right) (\Delta x_i)^2 + D_i \Delta x_i \quad \text{and} \quad f(x_{i+1}) = \left(\frac{z_{i+1}}{6}\right) (\Delta x_i)^2 + C_i \Delta x_i . $$

Solving for $C_i$ and $D_i$, then

$$ C_i = \frac{f(x_{i+1})}{\Delta x_i} - \frac{z_{i+1} \Delta x_i}{6} \quad \text{and} \quad D_i = \frac{f(x_i)}{\Delta x_i} - \frac{z_i \Delta x_i}{6}. $$

Substituting these values of $C_i$ and $D_i$ in (9), we get

$$ 
\begin{aligned}
p_i(x) &= \left(\frac{z_i}{6\Delta x_i}\right) (x_{i+1} - x)^3 + \left(\frac{z_{i+1}}{6\Delta x_i}\right) (x - x_i)^3 \\
&\quad + \left(\frac{f(x_{i+1})}{\Delta x_i} - \frac{z_{i+1}\Delta x_i}{6}\right) (x - x_i) + \left(\frac{f(x_i)}{\Delta x_i} - \frac{z_i\Delta x_i}{6}\right) (x_{i+1} - x).
\end{aligned} 
\quad (10) $$

To determine the values of the $z_i$, the property $p'_i(x_i) = p'_{i-1}(x_i)$ for $1 \le i \le n-1$ is used. This gives $n-1$ equations to determine the $n+1$ variables $z_i$ for $i = 0, 1, \dots, n$.

## B. Natural Spline

For the natural spline [2] interpolant, one sets $z_0 = z_n = 0$ and determine the values of the other $z_i$ using $p'_i(x_i) = p'_{i-1}(x_i)$ for $1 \le i \le n-1$.

From (10), then

$$ 
\begin{aligned}
p'_i(x) &= -\left(\frac{z_i}{2\Delta x_i}\right) (x_{i+1} - x)^2 + \left(\frac{z_{i+1}}{2\Delta x_i}\right) (x - x_i)^2 \\
&\quad + \left(\frac{f(x_{i+1})}{\Delta x_i} - \frac{z_{i+1}\Delta x_i}{6}\right) - \left(\frac{f(x_i)}{\Delta x_i} - \frac{z_i\Delta x_i}{6}\right)
\end{aligned} 
\quad (11) $$

for $i = 0, 1, \dots, n-1$. Hence,

$$ 
\begin{aligned}
p'_i(x_i) &= -\left(\frac{z_i}{2}\right) \Delta x_i + \left(\frac{f(x_{i+1})}{\Delta x_i} - \frac{z_{i+1}\Delta x_i}{6}\right) - \left(\frac{f(x_i)}{\Delta x_i} - \frac{z_i\Delta x_i}{6}\right) \\
&= -\frac{z_{i+1}\Delta x_i}{6} - \frac{z_i\Delta x_i}{3} + \frac{f(x_{i+1}) - f(x_i)}{\Delta x_i}.
\end{aligned} 
$$
Similarly,

$$p'_{i-1}(x) = -\left(\frac{z_{i-1}}{2\Delta x_{i-1}}\right)(x_i - x)^2 + \left(\frac{z_i}{2\Delta x_{i-1}}\right)(x - x_{i-1})^2 + \left(\frac{f(x_i)}{\Delta x_{i-1}} - \frac{z_i \Delta x_{i-1}}{6}\right) - \left(\frac{f(x_{i-1})}{\Delta x_{i-1}} - \frac{z_{i-1} \Delta x_{i-1}}{6}\right)$$

for $i = 1, 2, \dots, n$. Hence,

$$p'_{i-1}(x_i) = \left(\frac{z_i \Delta x_{i-1}}{2}\right) + \left(\frac{f(x_i)}{\Delta x_{i-1}} - \frac{z_i \Delta x_{i-1}}{6}\right) - \left(\frac{f(x_{i-1})}{\Delta x_{i-1}} - \frac{z_{i-1} \Delta x_{i-1}}{6}\right) \\ = \frac{z_i \Delta x_{i-1}}{3} + \frac{z_{i-1} \Delta x_{i-1}}{6} + \frac{f(x_i) - f(x_{i-1})}{\Delta x_{i-1}}.$$

The relation $p'_i(x_i) = p'_{i-1}(x_i)$ yields

$$\begin{aligned} & z_{i+1} \Delta x_i + 2z_i (\Delta x_i + \Delta x_{i-1}) + z_{i-1} \Delta x_{i-1} \\ &= \frac{6}{\Delta x_i} (f(x_{i+1}) - f(x_i)) - \frac{6}{\Delta x_{i-1}} (f(x_i) - f(x_{i-1})) \end{aligned} \quad (12)$$

for $i = 1, 2, \dots, n-1$. In conclusion, the $z_i$ for $1 \le i \le n-1$ are given by the solution of the $n-1$ dimensional linear system $Az = b$.

To evaluate the polynomial $p_i$ defined in (10), the nested form is used. Expanding $p_i$ around $x = x_i$, yields:

$$p_i(x) = \frac{z_{i+1} - z_i}{6\Delta x_i}(x - x_i)^3 + \frac{z_i}{2}(x - x_i)^2 + \left(-\frac{z_i \Delta x_i}{3} - \frac{z_{i+1} \Delta x_i}{6} + \frac{f(x_{i+1}) - f(x_i)}{\Delta x_i}\right)(x - x_i) + f(x_i).$$

Thus,

$$p_i(x) = ((\alpha_i(x - x_i) + \beta_i)(x - x_i) + \gamma_i)(x - x_i) + \delta_i, \quad (13)$$

where

$$\alpha_i = \frac{z_{i+1} - z_i}{6\Delta x_i}, \quad \beta_i = \frac{z_i}{2}, \quad \gamma_i = -\frac{z_i \Delta x_i}{3} - \frac{z_{i+1} \Delta x_i}{6} + \frac{f(x_{i+1}) - f(x_i)}{\Delta x_i}, \quad \delta_i = f(x_i) \quad (14)$$

### C. Bézier Curves

Bézier Curves are interpolating functions based on Bernstein [7] polynomials.

A Bernstein polynomial of degree $k$ is defined as:

$$B_i^k(u) = \frac{k!}{i!(k-i)!} u^i (1-u)^{k-i}, \quad i = 0 \dots k \quad (15)$$

It originates from the binomial expansion:

$$([1+u]-u)^k = \sum_{i=0}^{k} \frac{k!}{i!(k-i)!} u^i (1-u)^{k-i} = 1 \quad \text{leading to} \quad \sum_{i=0}^{k} B_i^k(u) = 1 \quad (16)$$

A Bézier curve of degree $k$ defined over $[0, 1]$ can be written with Bernstein polynomials as:

$$C(u) = \sum_{i=0}^{n} B_{i}^{k}(u) P_{i}, \quad u \in [0, 1]$$
The quadratic Bézier curve of degree $k = 2$:

$$C(u) = \sum_{i=0}^{n} B_{i}^{2}(u) P_{i}, \quad u \in [0, 1]$$

can be written with Bernstein polynomials:

$$\begin{aligned}
B_{0}^{2}(u) &= u^{2} \\
B_{1}^{2}(u) &= 2u(1-u) \\
B_{2}^{2}(u) &= (1-u)^{2}
\end{aligned} \qquad (17)$$

A cubic Bézier curve of degree $k = 3$ defined from a set of points $P_i$, $i = 0, \dots, n$. can be written in terms of the Bernstein polynomial:

$$C(u) = \sum_{i=0}^{n} B_{i}^{3}(u) P_{i}$$

Bézier curves $B_i^3(u)$, $i = 0, 1, 2, 3$ defined over $[0, 1]$ and displayed in Fig. 3 are given by:

$$\begin{aligned}
B_{0}^{3}(u) &= u^{3} \\
B_{1}^{3}(u) &= 3u(1-u)^{2} \\
B_{2}^{3}(u) &= 3u^{2}(1-u) \\
B_{3}^{3}(u) &= (1-u)^{3}
\end{aligned} \qquad (18)$$

<figure><img src="image_5.png" /><figcaption>Fig.3: Cubic Bernstein curves defined over [0,1].</figcaption></figure>

## D. B-Spline Interpolation

B-Spline i.e. Basis-Spline is a generalization of the spline construction toward use of a set of basis functions underlying the interpolation process.

A Bézier curve is a particular case of a B-spline possessing no interior knots and belongs to a particular family of parametric curves using Bernstein polynomials.

The major difference between Bézier and B-Splines, is that Bézier curves are defined by all control points *globally*, whereas B-Splines are defined by control points *locally*, meaning changes to one control point affect only a portion of the curve.
B-Splines allow for non-uniform knot spacing and offer local control over curve segments, making them more flexible for complex shapes.

B-splines are commonly used in 3D modeling applications to define smooth transitions between data points.

Bézier curves are widely used in graphic design, animation, and font design due to their simplicity and ease of control.

**B-spline basis functions:**

B-splines of degree 0 are piecewise constant functions which are continuous from the right. For $k > 0$, the B-splines of degree $k$ are piecewise polynomials of degree $k$ and class $C^{k-1}$.

There are two possible methods:

1. Cox-de Boor [6]: Multiple curves in successive 1D intervals as in Fig. 4.

<figure><img src="image_6.png" /><figcaption>Fig.4: B-spline curves of degree 3 in different intervals.</figcaption></figure>

2. Riesenfeld [8]: Single curve in 2D (or 3D) as in Fig. 5.

<figure><img src="image_7.png" /><figcaption>Fig.5: B-spline curves of degrees 1,2,3 and 4 in 2D. Control points are at ends of the dashed linear segments. Degree 2 needs 3 control points, whereas degree 3 needs 4 and degree 4 needs 5.</figcaption></figure>

In the Cox-de Boor [6] method:

* In 1D, $u_i$ are a set of ordered and non-decreasing knots, i.e. $u_0 < ...u_{i-1} < u_i < u_{i+1}... < u_m$ and $N_{i,k}(u)$ starts at knot $u_i$ and is of order $k$.

$N_{i,0}(u) > 0$ for $0 \le i \le m-1$ and $N_{i,0}(u) = 0$ otherwise. For $k \ge 1$, $N_{i,k}(u) > 0$ and $0 \le i \le m-1$; $N_{i,k}(u) = 0$ otherwise.

More precisely, these functions [9] can be defined with degree 0 as $N_{i,0}(u)$ and from recursion expressions for $N_{i,k}(u)$ of arbitrary degree $k$ according to Cox-de Boor [6] formula:

$$N_{i,0}(u) = \begin{cases} 1 & \text{if } 0 \le i \le m-1 \\ 0 & \text{otherwise} \end{cases}$$ $$N_{i,k}(u) = \frac{u-u_i}{u_{i+k}-u_i} N_{i,k-1}(u) + \frac{u_{i+k+1}-u}{u_{i+k+1}-u_{i+1}} N_{i+1,k-1}(u) \quad \text{if } k \ge 1 \text{ and } 0 \le i \le m-1$$

where $u_i$ are a set of ordered and non-decreasing knots, i.e. $u_0 < ...u_{i-1} < u_i < u_{i+1}... < u_m$ and $N_{i,k}(u)$ starts at knot $u_i$ and is of order $k$.
* In the Riesenfeld [8] method, a single curve in 2D (or 3D) is defined from a set of points $P_i$, $i = 0, ..n$. B-spline curves of degree k:

$$C(u) = \sum_{i=0}^{n} R_{i,k}(u)P_i, \quad u \in [0, 1]$$

$R_{i,k}(u)$ are Riesenfeld polynomials [8] defined by:

$$R_{i,k}(u) = (k+1) \sum_{j=0}^{k-i} (-1)^j \frac{(u+k-i-j)^k}{j!(k+1-j)!} \quad (19)$$

Various degrees of 2D B-spline curves are displayed in Fig. 5.

### Quadratic B-spline example with Cox-de Boor [6] method:

A quadratic [10] B-spline is a parametric curve defined by a set of knots $u_0 < ...u_{i-1} < u_i < u_{i+1}... < u_n$. Assuming that the knots are equidistant and given by the sequence: $0 < 1 < 2 < 3...$, Cox-de Boor [6] formula yields the following:

$$N_{0,2}(u) = \begin{cases} u^2/2 & \text{if } 0 \le u < 1 \\ (-3+6u-2u^2)/2 & \text{if } 1 \le u < 2 \\ (3-u)^2/2 & \text{if } 2 \le u < 3 \\ 0 & \text{otherwise} \end{cases}$$

$$N_{1,2}(u) = \begin{cases} (u-1)^2/2 & \text{if } 1 \le u < 2 \\ (-11+10u-2u^2)/2 & \text{if } 2 \le u < 3 \\ (4-u)^2/2 & \text{if } 3 \le u < 4 \\ 0 & \text{otherwise} \end{cases}$$

### 2D Cubic B-spline example with Riesenfeld [8] method:

A cubic B-spline is a parametric curve defined by $u \in [0, 1]$ and four control points: $P_0, P_1, P_2, P_3$ as shown in Fig. 6.

<figure><img src="image_8.png" /><figcaption>Fig.6: Cubic B-spline curve with four control points $P_0, P_1, P_2, P_3$.</figcaption></figure>

With this definition and the associations [11] $u = 0 \to P_0$ and $u = 1 \to P_3$, the cubic curve is expressed as:

$$C(u) = \sum_{i=0}^{3} R_{i,3}(u)P_i, \quad u \in [0, 1] \quad (20)$$

$P_0$ et $P_3$ are at the end tips of the curve, with point $P_1$ on the tangent line to curve at $P_0$ and $P_2$ on the tangent line to curve at $P_3$.
Riesenfeld cubic polynomials $R_{i,3}(u)$ define B-spline curves for $u \in [0, 1]$ and are displayed in Fig. 7. Using eq. 19, they are expressed as:

$$
\begin{align*}
R_{0,3}(u) &= \frac{1}{6}(1-u)^3 \\
R_{1,3}(u) &= \frac{1}{6}(4-6u^2+3u^3) \\
R_{2,3}(u) &= \frac{1}{6}(1+3u+3u^2-3u^3) \\
R_{3,3}(u) &= \frac{1}{6}u^3 \tag{21}
\end{align*}
$$

<figure><img src="image_9.png" /><figcaption>Fig.7: Cubic B-spline Riesenfeld functions.</figcaption></figure>

E. Hermite splines

Hermite splines provide an alternative to Bézier for representing cubic curves. Hermite polynomials allow a curve to be defined in terms of its endpoints and its derivatives at its endpoints.

Given a number of nodes $x_i$, $i = 1...n$, the Hermite polynomial interpolant of a function $f(x)$, is denoted as $H_{2n+1}(x)$, since this polynomial must have degree $2n + 1$ in order to satisfy the $2n + 2$ constraints:
$H_{2n+1}(x_i) = f(x_i), H'_{2n+1}(x_i) = f(x_i), i = 0, 1, ..., n$.

In order to satisfy these constraints for $i = 0, 1, ..., n$, write:

$$
\begin{align}
H_i(x) &= |L_i(x)|^2 (1 - 2L'_i(x_i)(x - x_i)), \nonumber \\
K_i(x) &= |L_i(x)|^2 (x - x_i) \tag{22}
\end{align}
$$

where $L_i(x)$ are the Lagrange interpolant defined above for nodes $x_i$, $i = 1...n$.

One can verify that $H_i(x_j) = \delta_{ij}$, $H'_i(x_j) = 0$ and $K_i(x_j) = 0, H'_i(x_j) = \delta_{ij}$ where $\delta_{ij}$ is the Kroenecker symbol.

Consequently, the Hermite interpolant is:

$$
H_{2n+1}(x) = \sum_{i=0}^{n} f(x_i)H_i(x) + f'(x_i)K_i(x) \quad (23)
$$

As an application, consider Cubic Hermite polynomials $H_0(u)$, $H_1(u)$, $H_2(u)$, and $H_3(u)$ chosen so that
$H_0(0) = 1, H_1(0) = 0, H_2(0) = 0, H_3(0) = 0, H'_0(0) = 0, H'_1(0) = 1, H'_2(0) = 0, H'_3(0) = 0,$
$H'_0(1) = 0, H'_1(1) = 0, H'_2(1) = 1, H'_3(1) = 0, H_0(1) = 0, H_1(1) = 0, H_2(1) = 0, H_3(1) = 1.$
The advantage of Hermite polynomials is that when there is a need for a cubic curve $C(u)$ that has value equal to $a$ at $u = 0$ and equal to $d$ at $u = 1$ and has first derivative equal to $b$ at $u = 0$ and $c$ at $u = 1$, then one can simply define the curve as $C(u) = aH_0(u) + bH_1(u) + cH_2(u) + dH_3(u)$.

Since a cubic polynomial is uniquely determined by its values and first derivatives at points $u = 0$ and $u = 1$, there is only one way to define the Hermite polynomials $H_i(u)$ to satisfy the preceding conditions. Thus cubic Hermite polynomials are:

$$H_0(u) = (1+2u)(1-u)^2 = 2u^3 - 3u^2 + 1, H_1(u) = u(1-u)^2 = u^3 - 2u^2 + u$$

$$H_2(u) = -u^2(1-u) = u^3 - u^2, H_3(u) = u^2(3-2u) = -2u^3 + 3u^2$$

The Hermite polynomials are scalar functions that can be used to define curves in $\mathbb{R}^k$ by using vectors as coefficients. This allows any degree three Bézier curve to be expressed in a Hermite fashion. In fact, it is easy to convert a Bézier curve $C(u)$ with control points $P_0, P_1, P_2$ and $P_3$ in $\mathbb{R}^k$ into a Hermite representation: because the initial derivative is $C(0) = 3(P_1 - P_0)$ and the ending derivative is $C(1) = 3(P_3 - P_2)$.

The Hermite representation curve is:

$$C(u) = P_0H_0(u) + 3(P_1 - P_0)H_1(u) + 3(P_3 - P_2)H_2(u) + P_3H_3(u).$$

Unlike Bézier curves, the Hermite representation of a curve is not a weighted average since the sum $H_0+H_1+H_2+H_3$ is generally not equal to 1. The coefficients of $H_0$ and $H_3$ are points (the starting and end points of the curve), but the coefficients of $H_1$ and $H_2$ are vectors. Cubic Hermite spline functions are displayed in Fig. 8.

<figure><img src="image_10.png" /><figcaption>Fig.8: Cubic Hermite spline functions.</figcaption></figure>

As a consequence, Hermite polynomials lack many of the properties of Bézier curves; their advantage, however, is that sometimes it is more natural to define a curve in terms of its initial and ending positions and derivatives than with control points. For the opposite direction, converting a Hermite representation of a curve,

$C(u) = r_0H_0(u) + r_1H_1(u) + r_2H_2(u) + r_3H_3(u)$, into a Bézier representation of the curve is also simple. Let $P_0 = r_0, P_3 = r_3, P_1 = P_0 + (1/3)r_1$, and finally $P_2 = P_3 - (1/3)r_2$.

**Note:**

Hermite spline polynomials should not be confused with Hermite orthogonal polynomials that are well-known in Quantum Mechanics as the Harmonic Oscillator [12] eigenfunctions:

$$H_0(x) = 1, H_1(x) = 2x, H_2(x) = -2 + 4x^2, H_3(x) = -12x + 8x^3,$$

$$H_4(x) = 12 - 48x^2 + 16x^4, H_5(x) = 120x - 160x^3 + 32x^5$$

**Subdivision algorithm**

Discovery of Bézier splines was made by Pierre Bézier [13], who was an engineer for Renault and Paul de Casteljau [14], who was an engineer for Citroën, both working for automobile companies.

The most fundamental algorithm for dealing with Bézier curves is the subdivision algorithm. This was devised in 1959 by Paul de Casteljau [14] and is referred to as "de Casteljau algorithm".
Consider a Bézier curve defined over the parameter $t \in [0, 1]$. It is possible to subdivide such a curve into two new Bézier curves, each belonging to a different domain by introducing a parameter $\tau$ such that one has two distinct intervals $0 \le t \le \tau$ and $\tau \le t \le 1$.

These two new Bézier curves defined over two different intervals, considered together, are equivalent to the single original curve from which they were derived.

One immediate value of de Casteljau algorithm is providing a numerically stable means of computing the coordinates of any point along the curve. Its extension to curves of any degree is thus made obvious.

<figure><img src="image_11.png" /><figcaption>**Fig.9:** Subdivision process of a Bézier curve according to de Casteljau algorithm. Adapted from Joy [10].</figcaption></figure>

It can be shown that if a curve is repeatedly subdivided, the resulting collection of control points converges to the curve. Thus, one way of plotting a Bézier curve is to simply subdivide it an appropriate number of times and plot the control polygons.

De Casteljau algorithm [14] works even for parameter values outside the original interval. Fig. 9 shows a Bézier curve subdivided. It is numerically stable as long as the parameter subdivision parameter is within the parameter domain of the curve. Successive De Casteljau divisions are illustrated in Fig. 10 showing convergence progress toward the curve.

<figure><img src="image_12.png" /><figcaption>**Fig.10:** Successive subdivision of a Bézier curve with de Casteljau algorithm showing convergence progress toward the curve. One adds a superscript 0 to the original control point and compute: $P_i^j = (1-\tau)P_{i-1}^{j-1} + \tau P_{i+1}^{j-1}$; $j = 1, \dots, n$, $i = 0, \dots, n-1$. The curve over $0 \le t \le \tau$ is defined using control points $P_0^0, P_0^1, P_0^2, \dots$ and the curve over $\tau \le t \le 1$ is defined using control points $P_1^n, P_1^{n-1}, P_2^{n-2}, \dots$. As an example, the figure displays at left, a value of $\tau = 0.4$ whereas at right, the higher value of $\tau = 0.6$ shows better convergence to the curve. Adapted from Joy [10].</figcaption></figure>

# VI. APPLICATIONS

Splines, Bézier and B-splines have numerous applications [13] in many areas of Science and Engineering [6]. In Mathematics, for instance, they allow to solve ordinary differential (ODE), partial differential equations (PDE) [6], optimization problems...

In Engineering, they are used in Digital Signal Processing, Imaging [9], Control, Filtering, Target Tracking... Their most useful applications in Computer Aided Design (in Architecture, Tool, Machine and Car body designs...), in Computer Graphics (Animations and special effects in Cinema), Digital Image Processing [15, 16] Map and Cartography drawings, Typography font design...
## A. Vector Graphics in font design and Typography

Vector graphics [15] are often used in virtual reality [16] and 3D modelling, as well as in Macromedia Flash animations (SWF: Shockwave Flash). A simple example of a vector graphic at work is font design. Fonts, or images used to represent them, are vector based, meaning that size increase from e.g. point 10 to point 28, the image increases in size without any degradation. Vector data is also widely used in architecture and cartography, such as in computer aided design [16] (CAD) or in geographic information systems (GIS).

A significant benefit of vector images [16] is that zooming to any size does not lead to degradation due to pixelization since they are curve geometry based as opposed to pixel based raster images.

<figure><img src="image_13.png" /><figcaption>**Fig.11:** (Left) TrueType circle with quadratic splines and (Right) Postscript circle with cubic splines. Anchor points (special control points) are on the circle whereas endpoints on each segment represent ordinary control points. In TrueType there are 8 anchor and 8 ordinary control points (since a quadratic curve needs three control points), whereas in Postscript there are 4 anchor and 8 control points (since a cubic curve needs four control points).</figcaption></figure>

A font is simply a collection of glyphs. But computer fonts contain more information related for instance as how to map an input string of bytes into an output display of visual glyphs. This mapping is called an encoding. Many fonts will also contain rules for how adjacent glyphs should be arranged. For instance, in fonts for the latin alphabet, there is often a rule saying that if "f" and "i" are adjacent then a special glyph, a ligature of "f" and "i" should be used instead. Fonts for latin letters generally do not require such rules (but can be enhanced by them), fonts for other scripts like arabic require them.

A character is a concept and a glyph is a graphical representation of that concept. For instance, the letter "a" is a character, however several glyphs may represent an "S" as displayed in Fig.12. In latin fonts there is often a single glyph for every character and a single character for every glyph, but not always. A ligature is an example of a glyph that represents two (or more) characters. In arabic most of the letters (characters) have at least four different glyphs, the appropriate glyph being determined by the letters around it.

<figure><img src="image_14.png" /><figcaption>**Fig.12:** "S" letter font glyph varieties with CM (Computer Modern), CM Italic, CM Sans-Serif, Palatino, Utopia and Times-Roman.</figcaption></figure>

Each glyph is a font possessing a shape however, there are various ways of describing that shape on a computer. PostScript and TrueType fonts generally describe the outline of the shape and then color in the interior of that outline, this coloring process is called rasterizing. On the other hand a bitmap font describes each glyph by drawing a picture of that character on a rectangular grid of pixels.

PostScript fonts use cubic Bézier splines, TrueType fonts use quadratic Bézier splines. SVG fonts (Scalable Vector Graphics, a language for describing two-dimensional vector and mixed vector/raster graphics) can be either.

An outline is composed of a set of contours or paths. In Fig. 13, the letter "S" Times-Roman outline is shown with its spline representation consisting of segments and splines.

Each cubic spline is defined by four points, two of which mark the end points of the spline and two more to describe the slope of the spline at those end points.
<figure><img src="image_15.png" /><figcaption>**Fig.13:** Vector graphics Postscript representation of Times-Roman letter "S" glyph outline. The font outline is in black while font contours in dashed blue are defined by a sequence of points: circles that lie on the outline are the anchors, while circles at the end points of the segment tangent to the outline at the anchor are ordinary control points. Contours and control points define a Bézier or B-spline interpolation for the font representation.</figcaption></figure>

The points which describe the slope are called control points. Moving a control point will change the shape of the spline, but it will not alter where the spline begins or ends.

A series of splines and lines are joined together to make a path. The three points are joined by splines or lines, but the type of joins alters the appearance of the connections.

In the first case there are no control points and the result is two straight lines, in the second case there is a line and a curved spline which are tangent at their point of intersection, in the third case there are two splines which are again tangent at their point of intersection. The point of intersection in the first case is called a corner point, and the incoming and outgoing splines or lines may be in any orientation with respect to each other. In the second case the point of intersection is called a tangent point, it is represented by a triangle. The third case is a curved point and is represented by a circle.

## B. Digital signal processing (DSP)

Continuous signals need the following three steps in order to be processed by digital computers:

* Convert analog signal to digital signal (A/D) by sampling.

* Process the digital signal (DSP)

* Convert the processed digital signal back to analog signal (D/A)

The discrete signals are processed by the digital filters [17] or more generally by various data processing algorithms carried out in the digital computers. A processing algorithm can be generally considered as a discrete impulse (unit sample) response $h[n]$ in time domain

$$y_d[n] = h[n] * x_d[n]$$

or its spectrum $H(\omega)$ in frequency domain:

$$Y_d(\omega) = H(\omega)X_d(\omega)$$

where $X$, $Y$ and $H$ are, respectively, the discrete Fourier transforms of $x$, $y$ and $h$. In particular,

$$X_d(\omega) = \mathcal{F}[x_d[n]] = \sum_{n=-\infty}^{\infty} x_d[n]e^{jn\omega}$$

which can be compared with the continuous Fourier transform of the sampled signal $x_s(t)$:

$$X_s(\omega) = \mathcal{F}[x_s(t)] = \int_{-\infty}^{\infty} \left[ \sum_{m=-\infty}^{\infty} x_c(mT)\delta(t - mT) \right] e^{-j\omega t} dt = \sum_{m=-\infty}^{\infty} x_c(mT)e^{-j\omega mT}$$
By comparing the above two forms of Fourier transform:

$$X_d(\omega) = X_s(\omega/T)$$

Of course this is not a surprise as the definition of discrete Fourier transform is based on the assumption that the sampling period is unity $T = 1$, i.e., the unit used for the time scale is the sampling period $T$:

$$X_d(\omega) = X_s(\omega)|_{T=1}$$

## 1. Reconstruction of Signal by Interpolation

In time domain, the reconstruction [17] $x_r(t)$ of the continuous signal $x(t)$ from its sampled version $x_s(t)$ can be considered as an interpolation process of filling the gaps between neighboring samples. The interpolation can be considered as convolution of $x_s(t)$ with a certain function $h(t)$:

$$x_r(t) = h(t) * x_s(t)$$

In frequency domain, the interpolation can be considered as a filtering process:

$$X_r(\omega) = H(\omega)X_s(\omega)$$

with the general effect of reserving the central portion of the periodic spectrum $X_x(\omega)$ while suppressing all its replica at higher frequencies.

### • Zero-order hold

$$h_0(t) = \begin{cases} 1 & 0 \le t < T_s \\ 0 & \text{otherwise} \end{cases}$$

A continuous signal $x_0(t)$ can be recovered by

$$x_0(t) = h_0(t) * x_s(t)$$

which is a series of square pulses with their heights modulated by $x(mT_s)$. The interpolation corresponds a low-pass filtering in frequency domain by

$$H_0(\omega) = \mathcal{F}[h_0(t)] = \frac{2}{\omega} \sin\left(\frac{\omega T_s}{2}\right) e^{-j\omega T_s/2}$$

### • First-order hold

$$h_1(t) = \begin{cases} 1 - |t|/T_s & 0 \le |t| < T_s \\ 0 & \text{otherwise} \end{cases}$$

A continuous signal $x_1(t)$ can be recovered by

$$x_1(t) = h_1(t) * x_s(t)$$

which is the linear interpolation of the sample train $x(mT_s)$ (connecting every two consecutive samples by a straight line). This interpolation corresponds a low-pass filtering in frequency domain by

$$H_1(\omega) = \mathcal{F}[h_1(t)] = \frac{4}{\omega^2 T_s} \sin^2\left(\frac{\omega T_s}{2}\right)$$

### • Ideal reconstruction

The reconstructed signals $x_0(t)$ and $x_1(t)$ using 0th or 1st order hold interpolation [17] are certainly different from the original signal $x(t)$, for the reason that the low-pass filter is non-ideal. To find the interpolation function for a perfect reconstruction of the original signal $x(t)$, consider an ideal low-pass filter in frequency domain:

$$H_{lp}(\omega) = \begin{cases} T_s & -\omega_c < \omega < \omega_c \\ 0 & \text{elsewhere} \end{cases}$$
with time domain impulse response

$$h_{lp}(t) = \mathcal{F}[H_{lp}(\omega)] = T_s \frac{\sin(\omega_c t)}{\pi t}$$

The reconstruction of the continuous signal $X(\omega)$ can now be realized by applying this ideal low-pass filter to the sampled signal $X_s(\omega)$:

$$X_r(\omega) = H_{lp}(\omega) X_s(\omega)$$

If the cutoff frequency $\omega_c$ is higher than $\omega_m$ but lower than $\omega_s - \omega_m$, then the central portion of $X_s(\omega)$ can be extracted and scaled by factor $T$, while all other replicas beyond the cutoff frequency are suppressed, i.e., the original signal is perfectly reconstructed:

$$X_r(\omega) = X(\omega)$$

In time domain, this perfect reconstruction is

$$\begin{aligned}
x(t) &= h_{lp}(t) * x_s(t) = T_s \frac{\sin(\omega_c t)}{\pi t} * \sum_{m=-\infty}^{\infty} x(mT_s)\delta(t - mT_s) \\
&= T_s \sum_{m=-\infty}^{\infty} x(mT_s) \left[ \delta(t - mT_s) * \frac{\sin(\omega_c t)}{\pi t} \right] \\
&= T_s \sum_{m=-\infty}^{\infty} x(mT_s) \frac{\sin(\omega_c(t - mT_s))}{\pi(t - mT_s)}
\end{aligned}$$

## 2. Digital-to-analog conversion (D/A)

Sometimes it is necessary to convert the processed signals back to continuous-time domain by interpolating the discrete data points in time domain, or, equivalently, a low-pass filter [17] in frequency domain. To do so, the processed sequence of numbers $y_d[n]$ is converted back into a train of impulses $y_p(t)$ (with time interval $T_s$ between two consecutive impulses), and then the time signal of the processed signal can be obtained by an ideal low-pass filter $H_{lp}(\omega)$ with cut-off frequency $\omega_c = \omega_s/2$:

$$Y_c(\omega) = H_{lp}(\omega) Y_p(\omega)$$

or in time domain an interpolation by a sinc($t$) = $\sin(t)/t$ function $h_{lp}(t)$

$$y_c(t) = h_{lp}(t) * y_p(t)$$

**Example:** The system carries out the following five operations sequentially:

* Sampling of the continuous signal:

$$x_s(t) = x(t) p(t) = \sum_{m=-\infty}^{\infty} x(mT_s)\delta(t - mT_s)$$

or in frequency domain

$$X_s(\omega) = \frac{1}{T_s} \sum_{k=-\infty}^{\infty} X_c(\omega - k\omega_s)$$

where $T_s$ is the sampling period and $\omega_s = 2\pi/T_s$ is the sampling frequency;

* Conversion [17] of the impulse train to a sequence of numbers

$$x[m] = x(mT_s)$$
* Digital filtering:

$$y[m] = h[m] * x[m], \quad \text{or in frequency domain} \quad Y(e^{j\omega}) = H(e^{j\omega})X(e^{j\omega})$$

* Conversion from sequence of numbers back to impulse train

$$y_s(t) = \sum_{m=-\infty}^{\infty} y[m]\delta(t - mT_s)$$

* Interpolation of impulse train to obtain a continuous-time signal

$$y_c(t) = y_s(t) * h_{lp}(t) = T_s \sum_{m=-\infty}^{\infty} y[m] \frac{\sin(\omega_c(t - mT_s))}{\pi(t - mT_s)}$$

or in frequency domain

$$Y_c(\omega) = H_{lp}(\omega)Y_s(\omega)$$

where $h_{lp}(t)$ and $H_{lp}(\omega)$ are respectively the impulse and frequency response functions of an ideal low-pass filter with cut-off frequency $\omega_c = \omega_s/2$.

Note that the spectrum of a sampled signal $x_s(t) = \sum_m x(mT_s)\delta(t - mT_s)$ is

$$X_s(\omega) = \int_{-\infty}^{\infty} x_s(t)e^{-j\omega_s t}dt = \int_{-\infty}^{\infty} \sum_{m=-\infty}^{\infty} x(mT_s)\delta(t - mT_s)e^{-j\omega_s t}dt = \sum_{m=-\infty}^{\infty} x(mT_s)e^{-j\omega_s mT_s}$$

whose period is $\omega_s = 2\pi/T_s$. If the sampling period and rate could be assumed to be unity $T_s = f_s = 1$ and $\omega_s = 2\pi$, then an alternative definition of discrete Fourier transform can be used to yield:

$$X_s(e^{j\omega}) = \sum_{m=-\infty}^{\infty} x[m]e^{-j\omega_m m}$$

with period $\omega_s = 2\pi$. This definition [17] is no longer valid if the sampling rate is specifically determined according to the sampling theorem by the highest frequency component $\omega_m$ contained in the signal $x(t)$, instead of arbitrarily assumed to be unity. In this case, the second definition is only valid if the argument $\omega$ is replaced by $\omega T_s$. Or, equivalently, the argument $\omega$ of the spectrum defined in the first equation has to be replaced by $\omega/T_s$. However, if the first definition is used throughout the analysis, such conversion is not needed.

### C. Computer Graphics and Surface interpolation

When interpolation is generalized to the surface case, it becomes a very powerful tool for computer graphics. The representation of a surface that is most frequently used in computer graphics, was independently discovered by Pierre Bézier and Paul de Casteljau, both working for automobile companies interested in the design of automotive parts and car bodies. These engineers initially developed a curve representation scheme that is geometrical in construction, and based upon polynomial functions.

Going from a curve to a surface called a Bézier patch requires a Cartesian (or Tensor) product providing $n \times m$ control points $\{P_{ij} : 0 \le i \le n, 0 \le j \le m\}$ and yielding a surface parameterized by two variables:

$$P(u, v) = \sum_{i=0}^{n} \sum_{j=0}^{m} B_{i,j}^{n}(u) B_{j,m}^{m}(v) \qquad (24)$$

where $B_{i,j}^n, B_{j,m}^m$ may be Bézier (or B-spline) functions with control points $n, m$ as displayed in Fig. 14.

Bézier and de Casteljau work on patch methodology has become the default standard for surface generation in computer graphics [9].
<figure><img src="image_16.png" /><figcaption>Fig.14: Smooth surface (in grey) obtained from Bézier spline interpolation of a linearly interpolated surface. Adapted from Buss [9].</figcaption></figure>

## APPENDIX A: ANALYTIC INTERPOLATION AND PLANCK RADIATION LAW

In this section, analytic interpolation made by Planck [18] is described. This interpolation allowed to solve the Black-body riddle and discover Quantum Mechanics [12, 19].

Black-Body radiation was noticed for the first time by Gustav Kirchhoff [20, 21] who used to watch the color change of holes present in heated metals worked by blacksmiths in his village. He noticed a systematic color change as the metal is being heated following a typical sequence from red to orange, yellow, white and finally to blue.

Hole color change did not depend on the nature of the metal but rather on its size.

Black-Body emission spectrum also called Kirchoff function shows a maximum at a wavelength $λ_{max}(T)$ such that $λ_{max}(T) ∝ 1/T$ where T is the temperature of the Black-Body (cf. Fig. 15). $λ_{max}(T)$ is the dominant color corresponding to the Black body temperature T.

<figure><img src="image_17.png" /><figcaption>Fig.15: Black body radiation law displaying a maximum at wavelength $λ_{max}(T)$ for every temperature T. $λ_{max}(T)$ shown by the dotted magenta line is the dominant color indicator of the Black body temperature. In Quantum Physics language, we have a photon gas in a cavity obeying Planck radiation law [21] with the (thermal) average spectral energy density $ρ(λ, T)$ of the Bose-Einstein form (the chemical potential is zero since the number of photons is not fixed), thus: $ρ(λ, T) = \frac{8πhc}{λ^5} \frac{1}{[\exp(\frac{hc}{λk_B T}) - 1]}$ gives the Black-body emission spectrum with λ the wavelength, h Planck constant and c velocity of light. Note that the human visible spectrum is limited roughly to the [400 nm: 800 nm] interval.</figcaption></figure>

Deriving Black-Body emission spectrum is based on assuming thermal equilibrium between the metallic material (considered as oscillators) and electromagnetic radiation since Maxwell theory established that light is as such. Rayleigh derived a simple relation between the spectral energy distribution $ρ = ρ(ν, T)$ where $ν, T$ are frequency and Black-Body temperature respectively and the mean oscillator energy $U(ν, T)$:

$$ρ(ν, T) = \frac{8πν^2}{c^3} U(ν, T), \qquad (A1)$$
Rayleigh [12] used classical theory based on the equipartition theorem yielding $U(\nu, T) = k_B T$ for the mean oscillator energy assuming equality of mean kinetic and potential energy densities, each equal to $k_B T/2$ where $k_B$ is Boltzmann constant.

Rayleigh did not obtain a maximum $\nu_{max} = c/\lambda_{max}$ but a divergence of the total radiation at all temperatures, with an ultra-violet catastrophe at small wavelengths i.e. divergence of $\rho(\nu, T)$ as $\nu \to 0$.

On the other hand, Wien in 1896, using the analogy with the velocity distribution formula of Maxwell for gas molecules, suggested the formula:

$$ \rho(\nu, T) = a\nu^3 \exp(-b\nu/T) \qquad (A2) $$

The law given by eq. A2 was thought for several years to fit well with experiments: namely, a power increase (scale invariance) at low $\nu$ followed by an exponential damping typical of many physical processes.

Constants $a$ and $b$ should have an universal character, and will play an important role in Planck interpolation based on reconciling Rayleigh formula and Wien formula.

From the early pioneering experiments performed during the 1870s, it was apparent that the total amount radiated grows like the fourth power of the absolute temperature [12]. This is obtained by summing [22, 23] over all radiation frequencies in the metallic hole, yielding Stefan-Boltzmann law [12, 19]:

$$ \int_{0}^{\infty} \rho(\nu, T)d\nu = \int_{0}^{\infty} a\nu^3 \exp(-b\nu/T)d\nu = \sigma T^4. \qquad (A3) $$

$\sigma$ is Stefan-Boltzmann constant given [24] by: $\sigma = 48\pi a/b^4 c^3$.

Taking account of $\lambda_{max}T$ = constant and Wien formula [12, 19]; combined with Stefan-Boltzmann result yields the scaling law:

$$ \rho(\nu, T) = \nu^3 f(\nu/T), \quad f(\nu/T) = a \exp(-b\nu/T) \qquad (A4) $$

reducing the dependence on frequency and temperature to a single universal formula $\nu/T$ and a scaling function $f(\nu/T)$ requiring constants $a, b$, from dimensional analysis.

In order to obtain Kirchoff formula, Planck needed to derive the mean oscillator energy $U(\nu, T)$ in a thermal bath at $T$. Planck took an indirect route, resorting to thermodynamics. Instead of guessing $U(\nu, T)$, he derived the entropy dependence over $U(\nu, T)$, and obtained by working backwards from eq. A2, a relation between the oscillator entropy $S$ and $U$ through the reciprocal function $(\partial^2S/\partial U^2)^{-1}$ such that:

$$ (\partial^2 S / \partial U^2)^{-1} = -b\nu U \qquad (A5) $$

as the simplest possibility.

When $\eta > 0$, this ensures entropy should increase with time [24] according to Thermodynamics second law. Together with $dS = dU/T$, leads easily to the form $\rho(\nu, T) = a\nu^3 \exp(-b\nu/T)$ allowing to recover Wien law.

Assuming universal character of the Kirchoff function, implies Wien law is universal as well as constants $a$ and $b$.

Planck realized later after additional experiments [24] that Wien law is not valid at long wavelengths requiring a quadratic term $U^2$ such that $(\partial^2S/\partial U^2)^{-1} \propto -U^2$, while at short wavelengths $(\partial^2S/\partial U^2)^{-1} \propto -U$.

In order to reconcile both behaviors: Rayleigh (at long wavelength) and Wien (at short wavelength), Planck decided to consider an interpolation formula containing both $U$ and $U^2$ contributions.

Experiments carried out along 1900 [24] at the longest infrared waves available (up to 60 microns) showed unequivocally that $U(\nu, T)$ is linear in $T$ for $T$ high.

Consequently one may suppose that $U(\nu, T) = \eta T$, and use second law $dS = dU/T$, to obtain:

$$ (\partial^2 S / \partial U^2)^{-1} = -U^2 / \eta \qquad (A6) $$

Interpolation consisted of adding up the linear term given by eq. A5 and the quadratic term given by eq. A6 such that:

$$ (\partial^2 S / \partial U^2)^{-1} = -\frac{1}{\eta} U(\eta b \nu + U) \qquad (A7) $$

to yield, with the boundary condition [24] $U \to \infty$ when $T \to \infty$:

$$ \frac{1}{T} = \frac{1}{b\nu} \ln(1 + \eta b\nu / U) \qquad (A8) $$
Finally $U(\nu, T)$ formula becomes:

$$U(\nu, T) = \frac{\eta b \nu}{\exp(b \nu/T) - 1} \qquad (A9)$$

This interpolated formula was in perfect agreement with more accurate experiments [24] performed later on and allowed to reveal the existence of Planck constant as $h = bk_B$. Moreover, eq. A9 allowed the introduction of Bose-Einstein distribution in Statistical Physics [12, 19].

[1] J. R. Rice, *Numerical Methods, Software and Analysis*, Second Edition, Academic Press, New-York (1993).

[2] W. H. Press, W. T. Vetterling, S. A. Teukolsky and B. P. Flannery, *Numerical Recipes in C: The Art of Scientific Computing*, Second Edition, Cambridge University Press (New-York, 1992).

[3] C. D. Meyer, *Matrix Analysis and Applied Linear Algebra*, SIAM, Philadelphia (2000).

[4] S. D. Conte, C. de Boor, *Elementary Numerical Analysis*, McGraw-Hill, New York, (1980).

[5] J. Stoer R. and Bulirsch, *Introduction to Numerical Analysis*, Second Edition, Springer-Verlag, (New York, 1992).

[6] C. de Boor, *A Practical Guide to Splines*, Revised Edition, Springer-Verlag, New York (2001)

[7] B. Dionne, *MAT3380 and graduate courses in Numerical Analysis*, Department of Mathematics and Statistics, University of Ottawa (2023).

[8] R. F. Riesenfeld, *Applications of B-splines to geometric problems of computer-aided design*, PhD Dissertation; Syracuse University (1973). See also: R. E Barnhill and R. F. Riesenfeld eds, *Computer-aided geometric design*, Academic Press, New-York (1974).

[9] S. R. Buss, *3D Computer Graphics*, Cambridge University Press, New-York (2003).

[10] K. I. Joy, *Geometric Modelling Notes*, University of California, Davis (2000)

[11] G. Dahlquist and A. Björk, *Equidistant Interpolation and the Runge Phenomenon*, Numerical Methods, pp. 101-103 (1974)

[12] L. D. Landau and E. M. Lifshitz, *Quantum Mechanics*, Vol. 3, Pergamon, Oxford (1965)

[13] P. Bézier, *The mathematical basis of the UNISURF CAD system*, Butterworth & Co (Publishers) Ltd, (1986).

[14] P. de Casteljau, *Courbes et Surfaces à Pôles*, Citroën Automobiles, Paris SA (1963).

[15] W. K. Pratt, *Digital Image Processing*, Wiley, New-York (1991)

[16] R. C Gonzalez and R. E. Woods, *Digital Image Processing*, Third Edition, Pearson Prentice Hall, Upper Saddle River N. J. (2008)

[17] R. Wang, *Introduction to Orthogonal Transforms*, Cambridge University Press, New-York (2011).

[18] M. Planck, *Annalen der Physik (Leipzig)* **37**, 642 (1912).

[19] L. D. Landau and E. M. Lifshitz, *Statistical Physics*, Vol. 5, Pergamon, Oxford (1975).

[20] G. Kirchhoff, *On the relation between the radiating and absorbing powers of different bodies for light and heat*, Phil. Mag. Series 4, Volume 20, 1-21 (1860).

[21] F. Grum and R. J. Becherer, *Optical radiation measurements, Volume 1: Radiometry; Volume 2: Color Measurements*, Academic Press (New-York, 1979).

[22] I. S. Gradstein and I. M. Ryzhik, *Table of Integrals, Series and Products*, Fifth Edition, edited by Allan Jeffrey (Academic Press) (1994).

[23] M. Abramowitz and I.S Stegun, *Handbook of Mathematical Tables*, Dover, New-York (1960).

[24] M. Nauenberg, *Max Planck and the birth of the quantum hypothesis*, Am. J. Phys. 84 (9), 709 (2016).