Hermite-Birkhoff interpolation on scattered data on the sphere and other manifolds

TL;DR

This paper introduces a novel Hermite-Birkhoff interpolation method for scattered data on spheres and manifolds, using basis functions dependent on geodesic distance that do not require solving linear systems.

Contribution

It proposes a new interpolation approach that constructs basis functions with specific derivative properties, avoiding linear system solutions and applicable to arbitrary manifolds.

Findings

Effective interpolation on spheres demonstrated through numerical tests.

Basis functions depend on geodesic distance and have zero derivatives at data points.

Method belongs to partition of unity class, ensuring flexibility and efficiency.

Abstract

The Hermite-Birkhoff interpolation problem of a function given on arbitrarily distributed points on the sphere and other manifolds is considered. Each proposed interpolant is expressed as a linear combination of basis functions, the combination coefficients being incomplete Taylor expansions of the interpolated function at the interpolation points. The basis functions depend on the geodesic distance, are orthonormal with respect to the point-evaluation functionals, and have all derivatives equal zero up to a certain order at the interpolation points. A remarkable feature of such interpolants, which belong to the class of partition of unity methods, is that their construction does not require solving linear systems. Numerical tests are given to show the interpolation performance.

Tables

Table 1. Table 1: MAEs and RMSEs computed on the sphere for f 1 subscript 𝑓 1 f_{1} .

	T0		T1		T2
$n$	MAE	RMSE	MAE	RMSE	MAE	RMSE
$500$	$2.87E-2$	$6.56E-3$	$2.96E-3$	$1.12E-3$	$1.89E-4$	$2.38E-5$
$1000$	$2.07E-2$	$4.09E-3$	$1.45E-3$	$5.44E-4$	$5.87E-5$	$6.83E-6$
$2000$	$1.14E-2$	$2.83E-3$	$5.88E-4$	$2.67E-4$	$1.73E-5$	$2.24E-6$
$4000$	$8.41E-3$	$1.94E-3$	$3.43E-4$	$1.34E-4$	$6.14E-6$	$7.27E-7$
$8000$	$6.30E-3$	$1.35E-3$	$1.83E-4$	$6.65E-5$	$1.98E-6$	$2.25E-7$
$16000$	$4.83E-3$	$9.08E-4$	$7.70E-5$	$3.28E-5$	$5.66E-7$	$7.66E-8$

Table 2. Table 2: MAEs and RMSEs computed on the sphere for f 2 subscript 𝑓 2 f_{2} .

	T0		T1		T2
$n$	MAE	RMSE	MAE	RMSE	MAE	RMSE
$500$	$1.44E-2$	$3.72E-3$	$4.98E-3$	$1.30E-3$	$2.00E-4$	$3.29E-5$
$1000$	$8.47E-3$	$2.37E-3$	$2.90E-3$	$6.36E-4$	$4.09E-5$	$1.08E-5$
$2000$	$6.78E-3$	$1.70E-3$	$1.22E-3$	$3.13E-4$	$1.72E-5$	$3.76E-6$
$4000$	$5.16E-3$	$1.20E-3$	$6.70E-4$	$1.55E-4$	$6.35E-6$	$1.28E-6$
$8000$	$3.36E-3$	$8.70E-4$	$2.52E-4$	$7.66E-5$	$2.16E-6$	$4.60E-7$
$16000$	$3.35E-3$	$6.20E-4$	$1.15E-4$	$3.78E-5$	$5.48E-7$	$1.57E-7$

Table 3. Table 3: MAEs and RMSEs computed on the sphere with lacunary data for f 1 subscript 𝑓 1 f_{1} . Left: missing a half of the 1st derivatives, right: missing a half of the 2nd derivatives.

	missing 1st der.		missing 2nd der.
$n$	MAE	RMSE	MAE	RMSE
$500$	$3.16E-2$	$5.68E-3$	$3.25E-3$	$9.17E-4$
$1000$	$2.20E-2$	$3.35E-3$	$1.80E-3$	$4.47E-4$
$2000$	$1.19E-2$	$2.13E-3$	$7.86E-4$	$2.23E-4$
$4000$	$8.28E-3$	$1.43E-3$	$4.02E-4$	$1.13E-4$
$8000$	$6.40E-3$	$9.38E-4$	$1.90E-4$	$5.73E-5$
$16000$	$4.90E-3$	$6.63E-4$	$9.30E-5$	$2.84E-5$

Table 4. Table 4: MAEs and RMSEs computed on the cylinder for f 1 subscript 𝑓 1 f_{1} .

	T0		T1		T2
$n$	MAE	RMSE	MAE	RMSE	MAE	RMSE
$500$	$4.01E-2$	$1.15E-2$	$4.22E-3$	$1.36E-3$	$1.38E-4$	$3.18E-5$
$1000$	$2.37E-2$	$6.19E-3$	$1.29E-3$	$5.22E-4$	$1.69E-5$	$5.67E-6$
$2000$	$1.27E-2$	$3.63E-3$	$8.71E-4$	$2.99E-4$	$1.06E-5$	$2.48E-6$
$4000$	$5.09E-3$	$1.87E-3$	$3.66E-4$	$1.34E-4$	$2.74E-6$	$6.79E-7$
$8000$	$4.74E-3$	$1.19E-3$	$1.27E-4$	$6.56E-5$	$4.79E-7$	$1.44E-7$
$16000$	$2.93E-3$	$1.07E-3$	$7.87E-5$	$3.55E-5$	$1.82E-7$	$6.81E-8$

Table 5. Table 5: MAEs and RMSEs computed on the cylinder for f 2 subscript 𝑓 2 f_{2} .

	T0		T1		T2
$n$	MAE	RMSE	MAE	RMSE	MAE	RMSE
$500$	$1.91E-2$	$5.13E-3$	$2.59E-3$	$9.23E-4$	$7.94E-5$	$2.44E-5$
$1000$	$6.16E-3$	$2.85E-3$	$1.09E-3$	$3.85E-4$	$1.74E-5$	$6.65E-6$
$2000$	$3.96E-3$	$1.55E-3$	$4.99E-4$	$1.79E-4$	$5.41E-6$	$1.75E-6$
$4000$	$2.73E-3$	$1.10E-3$	$2.21E-4$	$9.47E-5$	$1.59E-6$	$6.09E-7$
$8000$	$2.15E-3$	$7.33E-4$	$8.83E-5$	$4.28E-5$	$5.43E-7$	$1.83E-7$
$16000$	$1.23E-3$	$5.20E-4$	$5.66E-5$	$2.51E-5$	$1.93E-7$	$6.33E-8$

Table 6. Table 6: MAEs and RMSEs computed on the cone for f 1 subscript 𝑓 1 f_{1} .

	T0		T1		T2
$n$	MAE	RMSE	MAE	RMSE	MAE	RMSE
$500$	$2.29E-2$	$9.64E-3$	$3.20E-3$	$1.62E-3$	$7.80E-5$	$2.41E-5$
$1000$	$1.16E-2$	$5.64E-3$	$1.52E-3$	$7.35E-4$	$3.95E-5$	$9.93E-6$
$2000$	$6.60E-3$	$3.17E-3$	$1.01E-3$	$3.76E-4$	$1.06E-5$	$2.97E-6$
$4000$	$7.81E-3$	$2.51E-3$	$3.50E-4$	$1.72E-4$	$3.89E-6$	$9.76E-7$
$8000$	$4.11E-3$	$1.47E-3$	$1.68E-4$	$8.78E-5$	$8.80E-7$	$2.64E-7$
$16000$	$2.74E-3$	$9.02E-4$	$8.58E-5$	$4.50E-5$	$2.43E-7$	$8.61E-8$

Table 7. Table 7: MAEs and RMSEs computed on the cone for f 2 subscript 𝑓 2 f_{2} .

	T0		T1		T2
$n$	MAE	RMSE	MAE	RMSE	MAE	RMSE
$500$	$1.18E-2$	$3.90E-3$	$3.52E-3$	$1.35E-3$	$1.58E-4$	$4.32E-5$
$1000$	$9.06E-3$	$2.39E-3$	$1.93E-3$	$7.41E-4$	$4.31E-5$	$1.29E-5$
$2000$	$4.72E-3$	$1.31E-3$	$8.12E-4$	$3.24E-4$	$1.57E-5$	$5.55E-6$
$4000$	$3.47E-3$	$9.21E-4$	$3.36E-4$	$1.61E-4$	$6.62E-6$	$1.59E-6$
$8000$	$2.29E-3$	$6.63E-4$	$2.11E-4$	$9.30E-5$	$1.85E-6$	$4.63E-7$
$16000$	$1.10E-3$	$3.92E-4$	$1.22E-4$	$4.40E-5$	$5.27E-7$	$1.56E-7$

Equations

$$D^{\beta }H(z_i):=\frac{{\partial }^{|\beta |}H(z_i)}{\partial v_1^{{\beta }_1}\cdots \partial v_m^{{\beta }_m}}=f_{i\beta },\beta \in {\Delta }_i,i=1,\dots ,n,$$

$$D^{\beta }H(z_i)=D^{\beta }f(z_i),\beta \in {\Delta }_i,i=1,\dots ,n.$$

$$g_j\in C^k(U),g_j(u)\ge 0,\sum _{j=1}^n{g}_j(u)=1,g_j(z_i)={\delta }_{ji},$$

$$D^{\beta }{g}_j(z_i)=0,\beta \in {\Delta }_i,|\beta |\ne 0,i=1,\dots ,n.$$

$$\displaystyle H(u)=\sum_{i=1}^{n}~{}T\bigl{(}u;f,z_{i},\Delta_{i})~{}g_{i}(u),$$

$$\displaystyle T\bigl{(}u;f,z_{i},\Delta_{i}):=\sum_{\beta\in\Delta_{i}}\frac{D^{\beta}f(z_{i})}{\beta_{1}!\cdots\beta_{m}!}\big{(}v-v(z_{i})\big{)}^{\beta}=\sum_{\beta\in\Delta_{i}}\frac{D^{\beta}f(z_{i})}{\beta_{1}!\cdots\beta_{m}!}\big{(}v_{1}-v_{1}(z_{i})\big{)}^{\beta_{1}}\cdots\big{(}v_{m}-v_{m}(z_{i})\big{)}^{\beta_{m}},$$

$$H(u)=\sum _{i=1}^n\sum _{\beta \in {\Delta }_i}{D}^{\beta }f(z_i)g_{i\beta }(u),$$

$$g_{i\beta}(u):=\frac{\big{(}v-v(z_{i})\big{)}^{\beta}}{\beta_{1}!\cdots\beta_{m}!}g_{i}(u)=\frac{\big{(}v_{1}-v_{1}(z_{i})\big{)}^{\beta_{1}}\cdots\big{(}v_{m}-v_{m}(z_{i})\big{)}^{\beta_{m}}}{\beta_{1}!\cdots\beta_{m}!}g_{i}(u),$$

$$\displaystyle a)\quad\|H(u)\|\leq\max_{i}\|T\bigl{(}u;f,z_{i},\Delta_{i})\|,$$

$$\displaystyle b)\quad\|f(u)-H(u)\|\leq\sum_{i=1}^{n}g_{i}(u)\|f(u)-T\bigl{(}u;f,z_{i},\Delta_{i})\|\leq\max_{i}\|f(u)-T\bigl{(}u;f,z_{i},\Delta_{i})\|,$$

$$[D^{\beta }\alpha (u,z_i){]}_{u=z_i}=0,i=1,\dots ,n,0<|\beta |\le k.$$

$$g_i(u)=\frac{\prod _{j=1,j\ne i}^n\alpha (u,z_j)}{\sum _{k=1}^n\prod _{j=1,j\ne k}^n\alpha (u,z_j)},i=1,\dots ,n,$$

$$D^{\beta }{g}_i(z_j)=0,0<|\beta |\le k,i,j=1,\dots ,n.$$

$$g_i(u)=\frac{1/\alpha (u,z_i)}{\sum _{k=1}^{n}1/\alpha (u,z_k)},g_i(z_i)=1,i=1,\dots ,n,$$

$$\alpha (u,w)=\vartheta (d_g(u,w)),$$

$$\alpha (u,w)=(d_g(u,w){)}^{\mu },\mu \in {\mathbb{R}}^{+},\mu \ge k,u,w\in U,$$

$${\alpha }_{\gamma }(u,w)=\frac{\exp (\gamma (d_g(u,w){)}^{\mu })}{(d_g(u,w){)}^{\mu }},\gamma >0,\mu \ge k,$$

$${\alpha }_{\delta }(x,y)=\exp (\delta d_g(u,w{)}^{\mu }),\delta \ge 0,\mu \ge k,$$

$$g_i(u)=\frac{(d_g(u,z_i){)}^{-\mu }}{{\sum }_{k=1}^n(d(u,z_k){)}^{-\mu }},i=i,\dots ,n,$$

$${\tilde{g}}_i(u)=\frac{{\tau }_i(u)(d_g(u,z_i){)}^{-\mu }}{{\sum }_{k=1}^n{\tau }_k(u)(d_g(u,z_k){)}^{-\mu }},$$

$${\tau }_i(u)=\{\begin{array}{l}>0,foru:d_g(u,z_i)<\delta ,\\ =0,foru:d_g(u,z_i)\ge \delta ,\end{array}$$

$$\tau_{i}(u)=\Big{(}1-\frac{d_{g}(u,z_{i})}{\delta}\Big{)}^{k+1}_{+},\quad u\in U,\quad i=1,\ldots,n.$$

$${\tau }_i(u)=\{\begin{array}{l}1,foru:d_g(u,z_i)<\delta ,\\ 0,foru:d_g(u,z_i)\ge \delta .\end{array}$$

$$\displaystyle\tilde{H}(u)=\sum_{i=1}^{n}~{}T\bigl{(}u;f,z_{i},\Delta_{i})~{}\tilde{g}_{i}(u),$$

$$\displaystyle|f(u)-\tilde{H}(u)|=\Bigg{|}\sum_{i=1}^{n}~{}f(u)\,\tilde{g}_{i}(u)-\sum_{i=1}^{n}~{}T(u;f,z_{i},\Delta_{i})\,\,\tilde{g}_{i}(u)\Bigg{|}=$$

$$\displaystyle\Bigg{|}\sum_{i=1}^{n}~{}\big{[}f(u)-T(u;f,z_{i},\Delta_{i})\big{]}\,\,\tilde{g}_{i}(u)\Bigg{|}\leq\sum_{i=1}^{n}~{}\big{|}f(u)-T(u;f,z_{i},\Delta_{i})\big{|}\,\,\tilde{g}_{i}(u),$$

$$|f(u)-T(u;f,z_i,{\Delta }_i)|\le C_i\Vert v-v(z_i){\Vert }^{q+1},$$

$$|f(u)-T(u;f,z_i,{\Delta }_i)|\le C_i{d}_g^{q+1}(u,z_i).$$

$$|f(u)-\tilde{H}(u)|\le \sum _{i=1}^n{C}_i{d}_g^{q+1}(u,z_i){\tilde{g}}_i(u)\le {\delta }^{q+1}\sum _{i=1}^n{C}_i{\tilde{g}}_i(u)\le C{\delta }^{q+1},$$

$$h_{U,\mathcal{X}}:=\underset{u\in U}{\sup }\underset{z_i\in \mathcal{X}}{\inf }{d}_g(u,z_i).$$

$$|f(u)-\tilde{H}(u)|\le CK{h}_{U,\mathcal{X}}^{q+1}.$$

$$u_1=u_1(v_1,\dots ,v_d),\dots ,u_{d+1}=u_{d+1}(v_1,\dots ,v_d),$$

Full text

††journal: Applied Mathematics and Computation

1 Introduction

In previous papers we concerned with Hermite-Birkhoff interpolation of a function given on arbitrarily distributed points on Euclidean spaces [4] (see also e.g. [18, 19]), and with Lagrange and Hermite-Birkhoff interpolation of a function given on arbitrarily distributed points on Banach spaces [3, 5]. In particular, in [4] we considered data on a general domain $\Omega\subset\mathbb{R}^{d+1},~{}d\geq 1$, without using any topological information about $\Omega$. Nevertheless many applications provide us with additional information on the underlying domain. For example, problems coming from geology, geophysics, metereology, oceanography, satellite-based techniques, etc., often relate to the entire earth or a large part of it, so that the unit sphere would be an appropriate model and the additional information would lead to a better approximant. The considered data typically represent some physical phenomena, such as temperature, rainfall, pressure, and gravitational forces, measured at various points on the surface of the earth, possibly at different times.

The sphere is a particularly interesting example of a connected compact smooth manifold, even because considerations developed about interpolation on the sphere can be extended to other manifolds. We think convenient to discuss, as long as possible, a general framework, though in practice the most interesting manifolds are smooth two-dimensional manifolds, i.e. surfaces in $\mathbb{R}^{3}$. In fact, a problem that occurs frequently in science and engineering is to recover from a surface in three dimensions a real valued function that interpolates to a given set of data (see e.g. [37, 56]).

The generalized Hermite interpolations (or Hermite-Birkhoff interpolation, see e.g. [27], Chapter 19) on scattered data by means of basis functions depending on the distance have been developed since 1992, when appeared the pioneering paper by Wu [55]. Since then, the interest in this topic seems to have increased significantly (see e.g. the pertinent chapters in [26, 27, 53]). A number of authors have also considered the Hermite interpolation setting on scattered data on the sphere (see e.g. [30, 49, 25, 44, 38]) or even general Riemannian manifolds [23, 43].

The point of view we follow in this work is deeply different from that in the just quoted papers. The matter is that we do not consider a radial basis function method but a cardinal (radial) basis function method. In this way, to get the Hermite-Birkhoff interpolant we have not to solve a system of linear equations. In fact, the interpolant is directly expressed as a linear combination of basis functions, which depend on the geodesic distance, are orthonormal with respect to the point-evaluation functionals, and have all derivatives equal zero up to a certain order at the interpolation points. The coefficients of the linear combination are incomplete Taylor expansions of the interpolated function at the interpolation points. On the other hand, our interpolation method is strictly linked up with papers which discuss Lagrange interpolation by partition of unity methods, namely Shepard-like methods, on the sphere (see, in particular, [9, 10, 11, 12, 20, 21]). Our method also shows an interesting analogy with articles by Franssens [29, 28].

This paper is organized as follows. In Section 2 explicit expressions of Hermite-Birkhoff interpolants on manifolds are given. Since these definitions are based on a suitable class of cardinal basis functions, Section 3 describes a general way to construct such basis functions, which depend on geodesic distances on Riemannian manifolds. Moreover, upper bounds for errors in terms of the fill distance are showed. In Section 4 some basic considerations are pointed out on geometric objects and analytic tools to deal with. Beside sketching the general background, the main goal of these considerations is to motivate the strategies to be adopted in numerical computation. Sections 5 and 6 discuss the application of the proposed interpolants to numerical recovering of functions only known on scattered data on the sphere and other Riemannian manifolds. The reported numerical tests are restricted to surfaces, which are the most interesting and handy cases, but the adopted point of view is more general. Considering the sphere the well-known expression of the geodesic distance is available, but for other manifolds it is necessary to consider approximations of their geodesic distances.

2 Hermite-Birkhoff Interpolation on Manifolds

Now, we define Hermite-Birkhoff interpolation on Riemannian manifolds, in particular on the higher-dimensional sphere. The manifolds we consider are supposed to enjoy suitable properties, as specified in the following (see Sections 4–6 for geometric details).

Definition 2.1

Let us consider a $m$-dimensional Riemannian manifold ${\cal M}\subset\mathbb{R}^{d+1},~{}d\geq 2$, an open set $U=\{u\equiv(u_{1},\ldots,u_{d+1})\in\mathbb{R}^{d+1}\}\subset{\cal M}$, a function $\varphi:U\to\mathbb{R}^{d}$ which maps $U$ homeomorphically to the open set $V=\{v\equiv(v_{1},\ldots,v_{m})\in\mathbb{R}^{m}\}:=\varphi(U)$ so that $v=\varphi(u):=v(u)$ and $u=\varphi^{-1}(v):=u(v)$. Let ${\cal X}=\{z_{1},\ldots,z_{n}\}\subset U$ be a set of distinct points, possibly scattered, with associated finite sets $\Delta_{1},\ldots,\Delta_{n}\subset\mathbb{N}_{0}^{m}$. The Hermite-Birkhoff interpolation problem from $U$ to $\mathbb{R}$ consists in finding a function $H:U\rightarrow\mathbb{R}$ which satisfies the interpolation conditions

[TABLE]

where $\beta=(\beta_{1},\ldots,\beta_{m})$, $|\beta|=\beta_{1}+\cdots+\beta_{m}$, and the $f_{i\beta}\in\mathbb{R}$ are given values to be interpolated. It is assumed that $H\in C^{k}(U)$, where $k=\max\{|\beta|:\beta\in\Delta_{i},\,\hbox{for some }i,\,1\leq i\leq n\}$.

In the following, we will also think that the $f_{i\beta}$ are values assumed by an underlying function $f:U\rightarrow\mathbb{R}$, $f\in{\cal C}^{k}(U)$, so that the conditions (1) take the form

[TABLE]

In general, the values of $f$ and of some its derivatives are known only at the points of ${\cal X}$.

A constructive solution to the interpolation problem (1) can be given by introducing a suitable class of cardinal basis functions, which can be defined as follows.

Definition 2.2

Given a set of distinct points ${\cal X}=\{z_{i},~{}1\leq i\leq n\}$, arbitrarily distributed in the open set $U\subset{\cal M}$, the functions $g_{j}:U\to\mathbb{R},~{}1\leq j\leq n$, are cardinal basis functions with respect to ${\cal X}$ if they satisfy for all $u\in U$ the conditions

[TABLE]

where $\delta_{ji}$ is the Kronecker delta, and also satisfy the additional property

[TABLE]

It is clear that an interpolant based on these weights must be considered as a partition of unity method.

It can be easily checked that the following result holds.

Property 2.3

The interpolation conditions (2) are satisfied by the interpolant

[TABLE]

where

[TABLE]

is formally an incomplete Taylor expansion of $f$ at $z_{i}$, in the sense that it only includes the partial derivatives whose orders belong to $\Delta_{i}$. The interpolant (4) can also be expressed in the form

[TABLE]

where

[TABLE]

with $\beta\in\Delta_{i},~{}1\leq i\leq n$.

The formula (5) highlights that the interpolant is essentially constructed by using the $g_{i\beta}$ as basis functions.

The interpolant (4) enjoys the usual properties of cardinal basis interpolants.

Property 2.4

There hold the following inequalities:

[TABLE]

where the $i$-th term in the sum may be interpreted as the local error at the point $z_{i}$.

3 Cardinal Basis Functions on Manifolds

A classical way to construct cardinal basis functions defined on $\mathbb{R}^{d+1}$ is Cheney’s method (see [15] and [16], pp. 67–68), which can be used for manifolds as well, if we adopt a suitable distance.

Theorem 3.1

Let us consider $U\subset{\cal M}$ as in Definition 2.1 and let $\alpha:U\times U\rightarrow\mathbb{R}^{+}_{0}$ be a continuous and bounded function, such that $\alpha(u,z_{i})>0$ for all $u\in U$, $u\neq z_{i}$, and $\alpha(z_{i},z_{i})=0$ for all $z_{i}\in{\cal X}$. Moreover, let each $\alpha(u,z_{i})$ be $k$-times continuously differentiable on $U$ such that

[TABLE]

The corresponding cardinal basis functions

[TABLE]

are continuous and satisfy

[TABLE]

Proof: This result is essentially the $d-$dimensional case of the main theorem in [4]. $\square$

The $g_{i}$ in (7) can also be represented in the barycentric form

[TABLE]

which is more convenient from a computational point of view [2].

A natural choice is defining $\alpha$ using the distance between points. Since we are considering points on the manifold ${\cal M}$, we take the geodesic distance $d_{g}$ (see Definition 4.4 below) and define $\alpha$ in the general form

[TABLE]

which obviously must satisfy the assumptions of Theorem 3.1. In particular, in (9) we may consider the choice

[TABLE]

which ensures both the vanishing of the derivatives at the nodes and the regularity assumptions, and among the possible choices is the most direct. Other interesting choices are (see e.g. [9]):

[TABLE]

both of them being rapidly decreasing.

As a result of the choice (10), we obtain the cardinal basis functions

[TABLE]

but, for computational reasons, in many cases it may be preferable to use a localized version of the cardinal basis functions (11), that is,

[TABLE]

where $\tau_{i}:U\to\mathbb{R}^{+}_{0},~{}\tau_{i}\in{C}^{k}(U)$, such that

[TABLE]

and $\delta>0$ is a suitably chosen value. Hence, each function $\tilde{g}_{i}$ vanishes outside the neighborhood of $z_{i}$ such that $d_{g}(u,z_{i})\geq\delta,~{}u\in U$. It can be easily proved that the functions $\tilde{g}_{i}$, $i=1,\ldots,n$, are cardinal and enjoy the vanishing property on derivatives (3). The localization can be obtained by taking, for instance,

[TABLE]

A simpler, but a little naive, choice is

[TABLE]

For the Hermite-Birkhoff interpolant with cardinal basis functions (12)

[TABLE]

we can give more significant error estimates than for the basic case (4). Let $q\in\mathbb{N}$ be defined such that each Taylor-type expansion $T(u;f,z_{j},\Delta_{j})$ is a complete Taylor expansion up to order $q$, plus other terms of higher degree. For any $f:U\to\mathbb{R}$ with $f\in{C}^{q}(U)$ and for any $u\in U$, we have, since the cardinal basis functions $\tilde{g}_{i}$ are a partition of unity,

[TABLE]

each $\tilde{g}_{i}$ being non-zero only inside the ball of radius $\delta$ centered at $z_{i}$. Now, since each $T(u;f,z_{i},\Delta_{i})$ is a Taylor expansion complete up to order $q$, we can use the estimate

[TABLE]

where $C_{i}\in\mathbb{R}^{+}$ is a suitable constant and $\|\cdot\|$ is the Euclidean norm. Since $\|v-v(z_{i})\|$ is less than or equal to the geodesic distance $d_{g}(u,z_{i})$, it follows

[TABLE]

Inserting (17) in (15) and exploiting again the partition of unity property, since $|d_{g}(u,z_{i})|<\delta$, we get

[TABLE]

with $C=\max_{i}C_{i}$. Moreover, if we set the localization radius $\delta=Kh_{U,{\cal X}}$, where $h_{U,{\cal X}}$ is the so-called fill distance, that is,

[TABLE]

and $K\geq 1$, we obtain the estimate

[TABLE]

4 Some Facts from Geometry

Considering the Hermite-Birkhoff interpolation problem on the sphere $\mathbb{S}^{d}$ and Riemannian manifolds in $\mathbb{R}^{d+1}$, some basic considerations are to be pointed out on the geometric objects and analytic tools to deal with. In fact, besides giving the general background, they affect deeply numerical computation strategies.

Let us consider first the interpolation problem on the sphere. It is well known that the spherical earth model is navigated using flat maps, collected in an atlas, and no single flat map can represent the entire earth. Similarly, the sphere $\mathbb{S}^{d}\subset\mathbb{R}^{d+1}$ can be described using an atlas of charts, each chart mapping part of the sphere to a subset of $\mathbb{R}^{d}$. Precisely, for every $u\in\mathbb{S}^{d}$ there exist an open set $U\subset\mathbb{S}^{d}$ with $u\in U$ and a mapping $\varphi:U\to\mathbb{R}^{d}$ that maps $U$ homeomorphically to the open set $V:=\varphi(U)$. The pair $(U,\varphi)$ is a chart of $U$ and a collection ${\cal A}=\{(U_{\alpha},\varphi_{\alpha})\}$ of charts, which covers the sphere, is an atlas. Charts in an atlas may overlap and a single point of the sphere may be represented in different charts. Given two overlaping charts $(U_{\alpha},\varphi_{\alpha})$ and $(U_{\beta},\varphi_{\beta})$, a transition function, that is, a coordinate transformation, can be defined which goes from $\varphi_{\alpha}(U_{\alpha}\cap U_{\beta})\subset\mathbb{R}^{d}$ to the sphere and then back to $\varphi_{\beta}(U_{\alpha}\cap U_{\beta})\subset\mathbb{R}^{d}$. A chart $(U,\varphi)$ is of class $C^{k}$ if $\varphi^{-1}\in C^{k}(\varphi(U))$, whereas a $C^{k}$-atlas consists of $C^{k}$-charts and $C^{k}$-transition functions.

Through the chart $(U,\varphi)$ the neighborhood $U$ inherits the coordinates given on the Euclidean space $\mathbb{R}^{d}$ and the homeomorphism $\varphi$ leads us to describe $U$ as a locally Euclidean space. In fact, considering $\varphi^{-1}:V\to U$ we have that the coordinates $u_{1},\ldots,u_{d+1}$ of a point $u\in U$ can be given by $d+1$ parametric equations

[TABLE]

where the parameters $v_{1},\ldots,v_{d}$ identify the point $v=\varphi(u)\in V\subset\mathbb{R}^{d}$. Then, the map $\varphi^{-1}$ can be written in terms of its components as \varphi^{-1}(v_{1},\ldots,v_{d})=\big{(}u_{1}(v_{1},\ldots,v_{d}),\ldots, u_{d+1}(v_{1},\ldots,v_{d})\big{)}.

An atlas is not unique as the sphere can be covered in multiple ways using different combinations of charts. To describe a possible atlas for the sphere ${S}^{d}$, we consider for any $u^{*}\in{S}^{d}$ the open neighborhood of $u^{*}$ in ${S}^{d}$ given by $U^{+}_{u^{*}}=\{u\in{S}^{d}:(u,u^{*})>0\}$, where $(\cdot,\cdot)$ is the usual inner product in $\mathbb{R}^{d+1}$. Choosing the coordinate system in $\mathbb{R}^{d+1}$ so that the vector $u^{*}$ has, to say, components $(0,0,\ldots,1)$ we have that the neighborhood $U^{+}_{u^{*}}$ can be homeomorphically projected on an open set in $\mathbb{R}^{d}$. Similarly, we may consider the neighborhood $U^{-}_{u^{*}}=\{u\in{S}^{d}:(u,u^{*})<0\}$. Since ${S}^{d}$ can be thought as the union of a suitable number of charts, we get an atlas and the considered parametrization holds. In particular, considering the two-dimensional sphere $\mathbb{S}^{2}$ an atlas of six charts is obtained which covers the entire sphere. Otherwise for the sphere $\mathbb{S}^{2}$, as well as for $\mathbb{S}^{d}$, it may be sometimes convenient to use spherical coordinates (see e.g. [25, 43]). Choosing one or another atlas has significant effects, especially for actual applications and their numerical treatment.

Referring to a real function $\psi$ defined on the sphere $\mathbb{S}^{d}$, we say that $\psi$ is $k$-times differentiable on $\mathbb{S}^{d}$, or $\psi\in C^{k}(\mathbb{S}^{d})$, if $\psi\circ\varphi^{-1}\in C^{k}(\varphi(U))$ for every chart $(U,\varphi)$ of ${\mathbb{S}^{d}}$. The function $\psi:{\mathbb{S}^{d}}\to\mathbb{R}$ inherits the local coordinates $v_{1},\ldots,v_{d}$ of the chart $U$, since $\psi(u)=\psi\circ\varphi^{-1}(v)$.

The considerations on the sphere, just seen, can be extended to other manifolds. It is useful to recall a formal definition of a topological manifold (see e.g. [53, 7]).

Definition 4.1

*A set ${\cal M}\subset\mathbb{R}^{d+1}$ is called a topological manifold of dimension $m$, if it is a Hausdorff space with a countable basis of open sets such that for every $u\in{\cal M}$ there exist an open set $U\subset{\cal M}$ with $u\in U$ and a mapping $\varphi:U\to\mathbb{R}^{m}$ that maps $U$ homeomorphically to the open set $V:=\varphi(U)\subset\mathbb{R}^{m}$. The pair $(U,\varphi)$ is called a coordinate neighborhood of $u$ or a chart and for every $u\in U$ the vector $\varphi(u)=(v_{1}(u),\ldots,v_{m}(u))\in\mathbb{R}^{m}$ represents the local coordinates of $u$ in $V$. A chart is of class $C^{k}$ if $\varphi^{-1}\in C^{k}(\varphi(U))$. A collection ${\cal A}=\{(U_{\alpha},\varphi_{\alpha})\}$ of *$C^{k}-$*charts is called a $C^{k}-$atlas of ${\cal M}$ if the sets $U_{\alpha}$ cover ${\cal M}$ and, moreover, for any $U_{\alpha},U_{\beta}$ with $U_{\alpha}\cap U_{\beta}\neq\emptyset$ the transition maps $\varphi_{\beta}\circ\varphi_{\alpha}^{-1}$ and $\varphi_{\alpha}\circ\varphi_{\beta}^{-1}$ are in $C^{k}$ on $\varphi_{\alpha}(U_{\alpha}\cap U_{\beta})$ and $\varphi_{\beta}(U_{\alpha}\cap U_{\beta})$, respectively. Finally, a manifold ${\cal M}$ is called a $C^{k}-$manifold if it possesses a *$C^{k}-$atlas.

The smoothness of a function $f:{\cal M}\to\mathbb{R}$ is defined by the smoothness of $f\circ\varphi^{-1}$ with a chart $(U,\varphi)$, as it is pointed out in the following:

Definition 4.2

We say that $f:{\cal M}\to\mathbb{R}$ is $k$-times differentiable on ${\cal M}$ or $f\in C^{k}({\cal M})$ provided that for every chart $(U,\varphi)$ of ${\cal M}$ the composition

[TABLE]

is $k$-times differentiable.

It is important to realize that the definition of differentiability of a real-valued function on a $C^{k}$-manifold does not depend on the choice of the chart.

To introduce on a $C^{k}$-manifold ${\cal M}\subset\mathbb{R}^{d+1}$, $k\geq 1$, the notions of length and distance, each tangent space must be equipped with an inner product, so that it varies smoothly from point to point. The tangent space $T_{u}{\cal M}$ for a point $u\in{\cal M}$ is the space formed by the tangent vectors to all the curves in ${\cal M}$ passing through $u$. Here, a vector $\tau$ is a tangent vector in $u\in{\cal M}$ if there exists a differentiable curve $\gamma(t)$ on ${\cal M}$, depending on a parameter $t$ with $-\varepsilon\leq t\leq\varepsilon$ for some $\varepsilon>0$, such that $\gamma(0)=u$ and $\gamma^{\prime}(0)=\tau$, where

[TABLE]

It turns out that $T_{u}({\cal M})$ is a $m-$dimensional vector suspace of $\mathbb{R}^{d+1}$ and that a basis is given by

[TABLE]

It is interesting to note that the tangent space can be thought as the best linear approximation to ${\cal M}$ in $u$.

More explicity, let us consider a chart $(U,\varphi)$ of a manifold ${\cal M}$ with parametric equations

[TABLE]

and a curve $\gamma:[-\epsilon,+\epsilon]\to{\cal M}$ on $U$ whose equations in $\varphi(U)$ are

[TABLE]

Substituting the latter equations into those of the chart, we get the equations of the curve $\gamma(t)$ on $U$ as a function of $t$, that is, $\gamma(t)=u(v(t))$. Then, differentiating we obtain the tangent vector to the curve at the point $u=\gamma(0)$

[TABLE]

and the basis vectors are

[TABLE]

namely (21).

To operationalize the Hermite-Birkhoff interpolation technique we start considering the geometric problem involving the computation of lengths of curves lying on a surface ${\cal M}\subset\mathbb{R}^{3}$. The key idea is essentially based on replacing an infinitesimal element of a smooth curve by the corresponding element of its tangent plane. As a significant example, which concerns the sphere and other major surfaces, let us consider a surface ${\cal M}$, a curve $\gamma$ on it and a point $u=(u_{1}(v_{1},v_{2}),~{}u_{2}(v_{1},v_{2}),u_{3}(v_{1},v_{2}))$ on $\gamma$. Taking the arclength $s$ of the curve as a parameter, the vector $du/ds$ is of unit length and we have from (22) for $m=2$

[TABLE]

where $(\cdot,\cdot)$ is the scalar product. Making use of the notations

[TABLE]

we obtain the first fundamental quadratic form of the surface

[TABLE]

The components $g_{ij},(i,j=1,2),$ of the metric form the entries of a $2\times 2$ symmetric matrix, namely $ds^{2}$ is a positive quadratic form related to this matrix. The first fundamental quadratic form of a surface provides the expression for the length $ds$ of an infinitesimal arc and the length of a finite curve lying on the surface is obtained from it by integration. More precisely, if a curve on the surface is given by the equation $\gamma(t)=u(v(t)),~{}t_{1}\leq t\leq t_{2},$ its length is

[TABLE]

To handle this idea in a more general situation, we recall the concept of Riemannian manifold.

Definition 4.3

A $C^{k}$-manifold is called a $C^{k}-$Riemannian manifold if for every $u\in{\cal M}$ there exists a positive definite inner product $g_{u}:T_{u}({\cal M})\times T_{u}({\cal M})\to\mathbb{R}$ such that for every chart $(U,\varphi)$ the $m^{2}$ functions

[TABLE]

are in $C^{k}(V)$ with $V=\varphi(U)$. The family of $g_{ij}^{\varphi}(u)$, assuming compatibility among different charts, is called a Riemannian metric on ${\cal M}$, the $g_{ij}^{\varphi}(u)$ are the component of the metric and form the entries of a $m\times m$ symmetric matrix. The first fundamental quadratic form associated to the metric is

[TABLE]

All differentiable manifolds (of constant dimension) can be given the structure of a Riemannian manifolds. The Euclidean space itself carries a natural structure of Riemannian manifold, where the tangent spaces are naturally identified with the Euclidean space itself and the scalar product of the space is the standard scalar product. Precisely, with ${\partial u}/{\partial v_{i}}$ identified with the $i$-th standard basis vector $e_{i}=(0,\ldots,1,\ldots,0)$, the (canonical) Euclidean metric over an open subset $U\subset\mathbb{R}^{d+1}$ is defined by $g^{\mathrm{can}}_{ij}=(e_{i},e_{j})=\delta_{ij}$. Many familiar curves and surfaces, including for example all $d-$spheres, are specified as subspaces of a Euclidean space and inherit a metric from their imbedding.

Finally, we use the Riemannian metric to define the length of a curve on ${\cal M}$.

Definition 4.4

*Suppose that ${\cal M}$ is a connected *$C^{k}-$Riemannian manifold. Let $u,w\in{\cal M}$ be two distinct points and let $\gamma:[a,b]\to{\cal M}$ be a piecewise $C^{1}$ curve that connects these points, i.e., $\gamma(a)=u,\gamma(b)=w$. Then the length of $\gamma$ is expressible in one of the equivalent forms

[TABLE]

where the first integrand represents the length of an infinitesimal arc in the Riemannian metric and $\|\cdot\|$ denotes the norm induced by the inner product on the tangent space. Supposing ${\cal M}$ to be compact, we set $d_{g}(u,w)$ to be the infimum over the length of all such curves connecting $u$ and $w$. The shortest of such curves is called the shortest path for $u$ and $w$, and $d_{g}(u,w)$ is their geodesic or Riemannian distance.

If ${\cal M}=\mathbb{R}^{m}$ and if $g_{u}$ is the canonical inner product on $\mathbb{R}^{m}$ then our definition of the length of a curve in the Riemannian metric coincides with the classical definition. In this case, $d_{g}(u,w)=\|u-w\|_{2}$, i.e. the shortest curve between two points in $\mathbb{R}^{m}$ is the straight line between them. On the sphere, our definition of geodesic coincides with the old one, since both denote the length of the shorter portion of the great circle connecting the two points.

5 Numerical Results on the Sphere

In this section we discuss numerical calculation of Hermite-Birkhoff interpolation on the sphere. Referring to the framework in Sections 2 and 3, we develop our considerations on $\mathbb{S}^{d},~{}d\geq 2$, as long as possible.

5.1 Computation of the Geodesic Distance on the Sphere

The $d-$dimensional sphere represents a case where the concept of atlas is essential. As described in Section 4, suitable charts are given for $\mathbb{S}^{d}$ by the function $\varphi^{+}_{i}(u_{1},\ldots,u_{d+1})=(v_{1},\ldots,v_{i-1},v_{i+1},\ldots,v_{d})$, which projects the subset of $\mathbb{S}^{d}$ with $u_{i}>0$ in the subspace $\mathbb{R}^{d}$, and by the similar function $\varphi^{-}_{i}$ with $u_{i}<0$, being $i=1,\ldots,d+1$. The family of all these charts forms an atlas.

Referring to Definition 4.4, the geodesic distance of any two points $u$ and $w$ of the unit sphere $\mathbb{S}^{d}$, that is, the sphere of unit radius and centered at the origin, is

[TABLE]

The geodesic distance $d_{g}(u,w)$ is always expressible in terms of the Euclidean distance $d_{E}(u,w)$ in $\mathbb{R}^{d+1}$. In fact, we have

[TABLE]

and, conversely,

[TABLE]

It follows from the asymptotic expansion of $\arcsin x,~{}x\in\mathbb{R}$,

[TABLE]

that is, the difference between $d_{g}(u,w)$ and $d_{E}(u,w)$ may be very small if $u$ and $w$ are sufficiently close. Therefore, using our local interpolation method, the Euclidean distance may be considered as a good approximation of the geodesic distance.

In general, taking a radial function $\varphi(d_{E}(u,w))$ we obtain a zonal basis function $\psi(d_{g}(u,w))$ with $u,w\in\mathbb{S}^{d}$ by setting

[TABLE]

In particular, referring to the construction of cardinal basis functions in Section 3, we remark that the expressions of $\alpha$ in terms of Euclidean and geodesic distances are closely related, because from (9) we have for suitable functions $\vartheta$ and $\eta$

[TABLE]

Hence, known expressions of $\alpha$ in terms of the Euclidean distance can be used as well to get expressions of $\alpha$ in terms of the geodesic distance (see e.g. [31]).

5.2 Computation of the Interpolant on the Sphere

Implementing the interpolation formula $H:U\subset\mathbb{S}^{d}\to\mathbb{R}$ in (4), a problem is to optimize the nearest neighbour searching procedure for spherical data, that is finding in a convenient way the set of points of ${\cal X}$ closest to any point $u\in U$. It is possible to use a cell-technique, which consists in a space decomposition into hypercubic cells by overlaying a spatial grid on the sphere. The procedure, which is based on the optimized Renka’s algorithm for trivariate interpolation [48, 13, 14], has been successfully tested for $\mathbb{S}^{2}$ [20]. Referring only to $\mathbb{S}^{2}$, another procedure, as well successfully tested, makes a decomposition of $U\subset{\cal M}$ in strips or spherical zones [6, 10].

Generally speaking, it is convenient to use the localized version (12) of the cardinal basis functions, so that more distant points have less influence. Alternatively, one could consider exponential-type weights, as in (11), which are strongly decaying as distance increases (see e.g. [8], p. 46). The drawback is that at least one parameter is necessary for the localization and this implies the requirement of determining its optimal value. The choice of an appropriate value for the localization parameter $\delta$ in (13) determines the efficiency of the local scheme and is a nontrivial problem. In practice, the localization can be obtained using for interpolation only the nodes that belong to a convenient neighborhood of the point $u$ considered, i.e. the nodes $z_{k}$ whose distance from $u$ is smaller than $\delta$.

To test the performance of our interpolation method, we must consider in general a uniform distribution of nodes on the sphere. To generate uniformly (or quasi-uniformly) distributed random (or pseudo-random) points on the high-dimensional unit sphere one can use, in principle, anyone of the suitable algorithms proposed in the literature, but they are not quite equivalent and only some of them work for $d\geq 3$. The most efficient and fast algorithm is based on the fact that the normal distribution function for a point $u\in\mathbb{R}^{d+1}$ has a density that depends only on the distance of the point from the origin, so that the points of $\mathbb{S}^{d}$ have the uniform distribution (see e.g. [33, 40, 42, 51, 52]). Since our interpolation algorithm works locally and on suitable charts, we experimented another way to find interpolation points on the unit sphere $\mathbb{S}^{2}$ $\subset$ $\mathbb{R}^{3}$ (and in the following on manifolds), considering $n$ Halton points [54] on the spherical cap of $\mathbb{S}^{2}$ for $z>0.5$. Interpolation errors are instead evaluated on a nearly uniform distribution of $n_{eval}=50$ spiral points belonging to the considered chart (see e.g. [50, 10]). An example of interpolation and evaluation points defined on the chart of $\mathbb{S}^{2}$ is shown in Figure 1.

Therefore, to investigate accuracy of the Hermite-Birkhoff interpolant, we compute the Maximum Absolute Error (MAE) and the Root Mean Square Error (RMSE) given by

[TABLE]

In order to get an idea of the distribution of points in the set ${\cal X}\subset U$ (see Definition 2.1) and, in particular, of their uniformity and density, we consider two common indicators, that is, the separation distance

[TABLE]

and the fill distance (18). These parameters are crucial in order to investigate the accuracy of interpolation methods.

Considering a function $f$ to be recovered, it is possible, from a theoretical point of view, to take in account several combinations of function values and derivatives of $f$ and, moreover, this combinations may change from point to point of the set ${\cal X}$. In fact, the interpolant $H(u)$ in (4) is such as to offer a complete flexibility. On the other hand, in practice, the most interesting situation is when the values of $f$ and its first (and possibly second) derivatives are known at each point of ${\cal X}$. Hence, our numerical tests are restricted to this case and also to $\mathbb{S}^{2}$, as it is usually done using other interpolation methods, achieving the advantage of permitting a comparison among different schemes (see e.g. [25, 26]).

The test functions to be interpolated are taken from the restriction to $\mathbb{S}^{2}$ of the following trivariate functions

[TABLE]

Since the performance of the interpolant does not change significantly using other test functions (see e.g. [31, 47, 24]), for shortness we here report only the numerical results obtained considering $f_{1}$ and $f_{2}$. In Tables 1–2 we report the errors obtained for the interpolant (14) using a complete Taylor expansion up to order zero (T0), one (T1) and two (T2). From these tables we can observe the significant improvement (in terms of accuracy) of the interpolant (14) when making use of first and second derivatives in the Taylor expansion. Finally, to give an idea, in Table 3 we show the results obtained in case of lacunary data, that is when a half of the first and second derivatives respectively are missing. This study points out that the interpolation scheme results in an unavoidable loss of accuracy due to the lack of information, but in any case the method can be applicable.

The considered interpolation scheme for the sphere is suitable for parallel implementation (see [2, 17, 21]). In the implementation of the parallel algorithm for a distributed memory machine, the data are assigned to $p$ processors by breaking the set ${\cal X}$ into subsets ${\cal X}_{k},~{}k=1,\ldots,p$. In this way each processor proceeds to solve the interpolation problem on a subset ${\cal X}_{k}$. The parallel algorithm consists of three steps: a) partitioning and data distribution, so that each subset has an approximately equal number of points, b) local interpolation solving, after having determined the radius of influence for each point of the considered subset, c) data collection and evaluation phase, where each slave processor sends its partial results to the master processor. In the ideal case, when the algorithm is completely and efficiently parallelizable, the speed-up must be equal to the number of processors involved.

6 Numerical Results on Riemannian Manifolds

Moving from the considerations on the sphere to those on Riemannian manifolds in general, nothing changes for what concerns the structure of the interpolant (14) and of the cardinal basis functions to be used. Instead, what changes dramatically is the problem of calculating the geodesic distance between points, because on the sphere one has got a simple analytic expression of the geodesic distance while this does not happen for other manifolds.

6.1 Computation of the Geodesic Distance on Manifolds

In order to study properties of geodesics on a $m$-dimensional Riemannian manifold ${\cal M}$, it is convenient to consider a connected chart $(U,\varphi)$ and for each $u\in U$ the vector $\varphi(u)=(v_{1},\ldots,v_{m})\in V$ of local coordinates. Then, a geodesic on $U$ is a curve given by $m$ functions $v_{1}(s),\ldots,v_{m}(s)$ which satisfy the system of second order differential equations, called geodesic equations,

[TABLE]

where $s$ is the arclength parameter and $\Gamma_{ij}^{k}$ are the Christoffel symbols of the second kind. It is possible to express the Christoffel simbols in terms of the components $g_{ks}$ of the metric matrix and their derivatives as follows

[TABLE]

where $g^{ks}$ are the components of the matrix $(g^{ij})$, inverse of the matrix $(g_{ij})$. As the manifold has dimension $m$, the geodesic equations are a system of $m$ ordinary differential equations for the $m$ variables $v_{k}$. Thus, allied with initial conditions consisting of a point on the manifold and a tangent vector at the point, the system can theoretically be uniquely solved, at least locally, but actually there are serious difficulties.

Computing the geodesic distance is less prohibitive if we merely consider a particular, but important, subclass of Riemannian manifolds, namely the regular surfaces in $\mathbb{R}^{3}$ parametrically defined on an open set of $\mathbb{R}^{2}$ by a map of the type $u=(u_{1}(v_{1},v_{2}),u_{2}(v_{1},v_{2}),$ $u_{3}(v_{1},v_{2}))$. In this case the system (29) reduces to two equations.

A further simplification is achieved considering for the set $U$ on the surface ${\cal M}$ orthogonal local coordinates, so that $g_{12}=g_{21}=0$. For $\gamma=\gamma(s)$ to be a geodesic on $U$, then it is necessary and sufficient that the geodesic equations

[TABLE]

are satisfied, expressing now the Christoffel symbols by (30).

An even more favorable situation is achieved considering the Clairaut parametrizations. We say that an orthogonal parametrization is a Clairaut parametrization in $v_{1}$ if

[TABLE]

Similarly, we say that an orthogonal parametrization is a Clairaut parametrization in $v_{2}$ if

[TABLE]

The geodesic equations simplify in these cases to the $v_{1}$-Clairaut geodesic equations

[TABLE]

and to the $v_{2}$-Clairaut geodesic equations

[TABLE]

respectively. As an example, considering the torus

[TABLE]

the $v_{1}$-Clairaut geodesic equations are

[TABLE]

Actually the surfaces of revolution, which include many cases important for applications, appear to be the most manageable. Without loss of generalization, let us consider a plane $\pi\subset\mathbb{R}^{3}$ generated by the unit vectors $e_{1},e_{3}\in\mathbb{R}^{3}$, a straight line $l$ generated by the unit vector $e_{3}$, and a curve $c\in C^{k},~{}k\geq 1$, which is disjoint from $l$ and belongs to the positive halfplane with respect to $e_{1}$. Rotating the curve $c$ around the line $l$ we obtain a surface of revolution ${\cal M}$ with generating curve $c$ and revolution axis $l$. If $c$ is parametrically represented by

[TABLE]

the surface is given by

[TABLE]

For a fixed $v_{1}=v_{1}^{0}$ the curve $p={\cal M}(v_{1}^{0},v_{2})\subset\mathbb{R}^{3}$ is called a parallel of ${\cal M}$ and represents the circle of radius $\alpha(v_{1}^{0})$ obtained by rotating the point $(\alpha(v_{1}^{0}),\beta(v_{1}^{0}))\in c$ around the line $l$. Similarly, for a fixed $v_{2}=v_{2}^{0}$ the curve $m={\cal M}(v_{1},v_{2}^{0})\subset\mathbb{R}^{3}$ is called a meridian of ${\cal M}$ and is obtained by rotating $c$ of an angle $v_{2}^{0}$ around $l$.

If the curve is parameterized with respect to the arclength $s$, the differential equations of geodesics for surfaces of rotation are

[TABLE]

Important consequences of these equations are:

1. i)

the meridians of a surface of revolution are geodesic curves, 2. ii)

a parallel is a geodesic curve if and only if it is obtained by rotating a point on the generating curve whose tangent vector is parallel to the axis of revolution.

Note that the geodesic equations of surfaces of revolution parameterised with respect to the arclength are particular cases of $v_{1}$-parameterization of Clairaut.

We used Matlab to numerically solve equations of geodesics on a parametric surface and to picture the relevant graphs. If ${\cal M}\subset\mathbb{R}^{3}$ is the surface

[TABLE]

we built two programs, which are similar but meet different needs. One of them resolves the system (29) with $m=2$ (as well as its particular cases) and finds $(v_{1}(s),v_{2}(s))$ with the arclength parameter $s$ or a multiple of it, starting from the initial point $(v_{1}(0),v_{2}(0))$ and the derivatives $(dv_{1}(0)/ds,dv_{2}(0)/ds)$. Then the program draws the support of ${\cal M}$ by varying $(v_{1},v_{2})\in\varphi(U)$ and traces the geodesic curve leaving $s$ to vary in a given interval $[s_{i},s_{e}]\subset\mathbb{R}$ required in input. Seeing pictures of geodesics is obviously interesting and useful, but it is not our primary goal (see e.g. [45, 1]).

The other program resolves the system (29) with $m=2$ (as well as its particular cases) and finds $(v_{1}(s),v_{2}(s))$ with the arclength parameter $s$ or a multiple of it, given the initial point $(v_{1}(s_{i}),v_{2}(s_{i}))$ and the end point $(v_{1}(s_{e}),v_{2}(s_{e}))$. Then, the program traces the geodesic connecting the two points and, above all, compute the geodesic distance between them. This second program starts with an approximate path of the geodesic and improves the solution iteratively. Since the geodesic distance between the initial and the end points is very small, we can choose a segment as the initial guess. In general, the method works well, but in a few cases the convergence is not assured despite requiring compactness (see e.g. [39, 32]).

6.2 Computation of the Interpolant on Manifolds

To optimize the nearest neighbour searching procedure for data on a general Riemannian manifold, one can continue to use the techniques already described for the sphere. Of course, in individual cases, more difficulties than with the sphere may arise in the use of those procedures.

To test the performance of our interpolation method, we need to get a uniform (or quasi-uniform) distribution of nodes on the considered Riemannian manifolds. Unfortunately, finding a convenient distribution is another critical point. To face the problem of generating a uniform distribution of points on analytic surfaces, there appear to be interesting the results of some recent papers [36, 35, 34, 41, 46]. On the other hand, to get information on uniformity and density of the distribution of points in the set ${\cal X}\subset U$, the separation distance (28) and the fill distance (18), already considered for the sphere, continue to be crucial parameters to assess the accuracy of interpolation methods.

To test our interpolant on manifolds, we focus on cylinder and cone. As interpolation nodes, we take some sets of $n$ uniformly random Halton data points, originally contained in the unit square $[0,1]^{2}\subset\mathbb{R}^{2}$ and then mapped onto the surface of cylinder and cone via, respectively, the equations

[TABLE]

and

[TABLE]

where $(p,q)\in[0,1]^{2}$, $r$ is the radius and $h$ denotes the height of the cone. The computation of interpolation errors, using (27), is carried out mapping as earlier a set of $n_{eval}=50$ evaluation points generated by the rand command of Matlab. Note that in our tests we consider a chart of ${\cal M}$, taking all points belonging to the cylinder for $x<-0.5$ and to the cone for $x<0$, assuming $r=1$ and $h=2$. An example of interpolation and evaluation points defined on the charts of cylinder and cone is given in Figure 2.

In order to recover a function $f$, known on a set ${\cal X}\subset U\subset{\cal M}$ together with some of its derivatives, the interpolant $H(u)$ in (4) is quite efficient for any combination of derivatives. However, actually, the most interesting situation is when the values of $f$ and its first and second derivatives are known at each point of ${\cal X}$. Hence, our numerical tests are restricted to this case and also to a chart of the cylinder and a chart of the cone. As regard to Hermite-Birkhoff interpolation on manifolds it is difficult to find numerical tests in the literature, as far as we know, while theoretical considerations are not lacking (see e.g. [22, 23, 43]). The test functions to be interpolated are taken from the restriction to $U\subset{\cal M}$ of the trivariate functions already considered for the sphere, but we report only the numerical results obtained considering $f_{1}$ and $f_{2}$. Interpolation errors computed for the interpolant (14) using a complete Taylor expansion up to order zero (T0), one (T1) and two (T2) are shown in Tables 4–5 for the cylinder and Tables 6–7 for the cone. From these numerical experiments we obtain an error behavior similar to that observed in the previous section for the sphere.

The considered interpolation schemes for the Riemannian manifolds are suitable for parallel implementation as explained for the sphere.

Acknowledgements

This work was supported by the University of Turin via grant “Metodi numerici nelle scienze applicate”.