Adaptive Rational Splines

Robert Schaback

Abstract. The classical interpolation problems for cubic and rational splines are merged to get an “adaptive” rational interpolating spline which automatically uses cubic pieces to model unavoidable inflection points and retains convexity/concavity elsewhere. An existence proof, a numerical method, and a series of examples are presented. Furthermore, the two– dimensional case is discussed.

Keywords. Rational functions, convexity–preserving splines.

1 Introduction

For a partition

a = x_{0} < x_{1} < ... < x_{n+1} = b, n ≥ 0 (1)

of a real interval [a, b] into subintervals I_{j} = [x_{j}, x_{j+1}], 0 ≤ j ≤ n, we want to construct a twice differentiable function f on [a, b] that is either a cubic polynomial or a rational function of the form

f_{j}(x) = \frac{a_{j} + b_{j}x + c_{j}x^{2}}{1 + d_{j}x}  (2)

on I_{j} and which satisfies the interpolation conditions

f(x_{j}) = y_{j} 0 ≤ j ≤ n

for given data values y_{0}, ..., y_{n} at x_{0}, ..., x_{n}. In addition, we impose a boundary condition of the type

α) f′(x) = y′ or

β) f′′(x) = y′′

on each of the endpoints x_{0} = a, x_{n+1} = b. Note that the parameters d_{j} are allowed to vary freely. This distinguishes our approach to rational splines from others (see e.g. [5]).

Due to the special form of (2), f′′(x) can have zeros only on intervals I_{j} where f(x) = f_{j}(x) is a cubic polynomial. We call such intervals cubic, the others rational. To avoid unnecessary inflection points, we make the number of cubic intervals as small as possible and place them where inflection points are necessarily induced by the data.

Mixed type nonlinear spline interpolation problems were treated by [3], but an existence proof was given there only for the case of sufficiently large n, using convergence of rational splines to cubic splines. This paper gives a general proof and a series of numerical examples. Approximation by rational splines is treated in [1], [4], and [7].

2 Basic equations

Following [2] and [3] we write the interpolation problem as a system of equations in terms of the second derivatives M_{j} := f”(x_{j}), 0 ≤ j ≤ n + 1, at the knots, involving the given data

\begin{aligned}h_{j} & := & x_{j+1} − x_{j} & (0 ≤ j ≤ n) \\ Δ_{j}^{1} & := & \frac{1}{h_{j}} (y_{j+1} − y_{j}) & (0 ≤ j ≤ n) \\ D_{j} & := & Δ_{j}^{1} − Δ_{j−1}^{1} & (1 ≤ j ≤ n). \end{aligned}

The quantities D_{j}, 0 ≤ j ≤ n + 1, are positive multiples of f′′. Therefore we use them to determine the type of the local interpolant from properties of the data alone:

Definition 2.1 A subinterval I_{j} = [x_{j}, x_{j+1}] is called rational iff

D_{j} · D_{j+1} > 0

and cubic otherwise. A knot x_{j} is called rational if it is the boundary point of at least one rational subinterval. Furthermore, let

R := {x_{j} | 0 ≤ j ≤ n + 1, x_{j} is rational }

denote the set of rational knots; the other knots will be called cubic.

We require the interpolant to be rational in rational subintervals and cubic in cubic subintervals. Strict convexity or concavity of the interpolant will thus be conserved as far as possible; inflection points will occur only in cubic intervals which will be placed near to points where the second derivative necessarily must have a zero for any C^{2} interpolant of the data.

We now express first–order derivatives f′_{j}(x_{j}), f′_{j}(x_{j+1}) by the data and the unknowns M_{j} := f”(x_{j}) :

a) For rational intervals I_{j} = [x_{j}, x_{j+1}] we get

\begin{aligned}f′_{j}(x_{j}) & = Δ_{j}^{1} − \frac{1}{2} h_{j}M_{j}^{2/3}M_{j+1}^{1/3} \\ f′_{j}(x_{j+1}) & = Δ_{j}^{1} + \frac{1}{2} h_{j}M_{j}^{1/3}M_{j+1}^{2/3}, \end{aligned}  (3)

using (2) as a representation of f_{j} in I_{j} and after elimination of the parameters a_{j}, b_{j}, c_{j}, and d_{j} in favor of the data and M_{j}, M_{j+1} (see e.g. [2]).

b) For cubic intervals we get

\begin{aligned}f′_{j}(x_{j}) & = & Δ_{j}^{1} & − & \frac{1}{6} h_{j}M_{j+1} & − & \frac{1}{3} h_{j}M_{j} \\ f′_{j}(x_{j+1}) & = & Δ_{j}^{1} & + & \frac{1}{6} h_{j}M_{j} & + & \frac{1}{3} h_{j}M_{j+1}, \end{aligned}  (4)

if f_{j} is written as a cubic polynomial in I_{j}.

Reference [3] contains a general treatment of interpolating C^{2} spline functions composed of pieces from certain 4–parameter nonlinear families. Equations (3) and (4) are special cases of a general formula given in [3].

We have to enforce continuity of the first derivative on the boundary point x_{j} between intervals I_{j−1} = [x_{j−1}, x_{j}], I_{j} = [x_{j}, x_{j+1}], 1 ≤ j ≤ n. Combining two equations of the form (3) or (4) we get the following equations for the unknowns M_{j} := f”(x_{j}) :

Case 1 I_{j−1}, I_{j} rational

\frac{1}{2} h_{j−1}M_{j−1}^{1/3}M_{j}^{2/3} + \frac{1}{2} h_{j}M_{j+1}^{1/3}M_{j}^{2/3} = D_{j}  (5)

Case 2 I_{j−1} cubic, I_{j} rational (or vice versa, analogously)

\frac{1}{6} h_{j−1}M_{j−1} + \frac{1}{3} h_{j−1}M_{j} + \frac{1}{2} h_{j}M_{j+1}^{1/3}M_{j}^{2/3} = D_{j}  (6)

Case 3 I_{j−1}, I_{j} cubic

\frac{1}{6} h_{j−1}M_{j−1} + \frac{1}{3} (h_{j−1} + h_{j})M_{j} + \frac{1}{6} h_{j}M_{j+1} = D_{j}  (7)

Boundary conditions are treated as follows:

α) f′_{0}(x_{0}) = y′_{0} implies an equation

\frac{1}{6} h_{0}M_{1} + \frac{1}{3} h_{0}M_{0} = D_{0} := Δ_{1}^{1} − y′_{0}

for I_{0} cubic and

\frac{1}{2} h_{0}M_{0}^{2}M_{1} = D_{0} := Δ_{1}^{1} − y′_{0}

for I_{0} rational. These can easily be incorporated by setting h_{−1} := 0 and neglecting the corresponding terms in an additional equation at x = x_{0}. Obviously, the resulting equation has one of the above forms.

β) f′′_{0}(x_{0}) = y′′_{0} = M_{0} simply fixes M_{0} but leaves the rest as it was. We formally define D_{0} := M_{0} in this case.

In each case the form of the equations is not changed.

3 Cubic sections

For the proof of the main existence theorem in the next section we have to impose an additional restriction on “cubic sections”; i.e., maximal sequences of adjoining cubic intervals

[x_{i}, x_{i+1}], [x_{i+1}, x_{i+2}], ..., [x_{j−1}, x_{j}],  0 ≤ i ≤ j − 1 ≤ n.  (8)

On such a cubic section C := [x_{i}, x_{j}] the inequalities

D_{k} · D_{k+1} ≤ 0 (i ≤ k ≤ j − 1)  (9)

hold, and the “local boundary points” x_{k} ∈ ∂C are either boundary points of the original problem (i.e., i = k = 0 or j = k = n + 1) or rational points (i.e., boundary points to a neighbouring rational interval, x_{k} ∈ ∂C ∩ R).

Definition 3.1 Let s_{M} be a cubic spline on a cubic section C of the form (8) satisfying

1. s_{M}(x_{k}) = y_{k}, i ≤ k ≤ j,

2. standard boundary conditions on non–rational boundary points of C,

3. boundary conditions s′′_{M}(x_{k}) = M_{k} on rational boundary points x_{k} of C, using a boundary derivative value M_{k} with sgn M_{k} = sgn D_{k}.

A cubic section C of the form (8) is called nondegenerate if for any spline s_{M} defined as above the inequalities

M_{i} sgn s′′_{M}(x_{i+1}) ≤ 0 if x_{i} ∈ ∂C ∩ R
M_{j} sgn s′′_{M}(x_{j−1}) ≤ 0 if x_{j} ∈ ∂C ∩ R (10)

hold.

This makes sure that the second derivative actually changes sign when going from a rational boundary point to a cubic neighbor. A sufficient criterion for nondegeneracy is given by

Theorem 3.1 If the data in a cubic section (8) satisfy

(−1)^{k−i}D_{k}D_{i} ≥ 0 for i ≤ k ≤ j  (11)

(i.e., weak alternation of D_{k}), the cubic section is nondegenerate.

Proof: We consider the Gauss–Seidel iteration on the standard system of equations (7) for the cubic spline in terms of second derivatives. If started with M_{k} = D_{k}, it will satisfy

(−1)^{k−i}M_{k}M_{i} ≥ 0 for i ≤ k ≤ j

for all iterates, if (11) holds.

Sign distributions of D_{i}, ..., D_{j} like (+, 0, 0, +) are not covered by Theorem 3.1. They have to be removed by slight changes of the data in order to make our theory applicable, although they do not seem to cause problems in practical computations.

We shall need the following bound of the second derivatives of the solution at the knots in cubic sections:

Theorem 3.2 There is a constant C_{1} = C_{1}(D_{i},…, D_{j}) such that

|s′′_{M}(x_{k})| ≤ C_{1} + \frac{1}{2} max \{|M_{m}| | \begin{matrix}m = i or j, \\ x_{m} a rational \\ boundary point \end{matrix} \} , i < k < j,

for the cubic interpolant s_{M} from Definition 3.1.

Proof: Let s_{M} = g_{D} + g_{M}, where g_{D} interpolates for homogeneous boundary values on rational boundary points and inhomogeneous data values D_{i}, ..., D_{j}, while g_{M} interpolates for inhomogeneous boundary values and homogeneous data values. Then for C_{1} := max_{i≤k≤j} |g′′_{D}(x_{k})| we have

|s′′_{M}(x_{k})| ≤ C_{1} + |g′′_{M}(x_{k})| (i ≤ k ≤ j).

From the equations (see (7) )

\frac{1}{6} h_{k−1}g′′_{M}(x_{k−1}) + \frac{1}{3} (h_{k−1} + h_{k})g′′_{M}(x_{k}) + \frac{1}{6} h_{k}g′′_{M}(x_{k+1}) = 0

at interior points x_{k} we get

−g′′_{M}(x_{k}) = \frac{h_{k}}{h_{k−1} + h_{k}} \frac{1}{2} g′′_{M}(x_{k+1}) + \frac{h_{k−1}}{h_{k−1} + h_{k}} \frac{1}{2} g′′_{M}(x_{k−1})

and this convex combination implies

|g′′_{M}(x_{k})| ≤ \frac{1}{2} max(|g′′_{M}(x_{k+1})|, |g′′_{M}(x_{k−1})|).

Extending this inequality recursively to the local boundary yields the assertion.

4 Existence of adaptive rational splines

The main result of this paper is

Theorem 4.1 There exists a unique interpolating twice differentiable adaptive rational spline, if all cubic sections determined by the data via D_{i} are nondegenerate.

Proof: We solve the equation at x_{j} for M_{j} and thus get a nonlinear Gauss–Seidel iteration scheme as a mapping T : K ⟶ K which will be proven to satisfy the hypothesis of Brouwer’s fixed point theorem on a suitable compact convex domain K. To define the latter we need some notation and a series of a–priori estimates.

The iteration will ensure that

sgn s′′(x_{j}) = sgn M_{j} = sgn D_{j} ≠ 0

on rational points; i.e., for j ∈ R. Using constants satisfying

0 < d ≤ |D_{j}| ≤ D for all j ∈ R
0 < h ≤ h_{j} ≤ H 0 ≤ j ≤ n,
h ≤ 6D

we set

c := \frac{d^{2}h}{2DH^{2}} C = max \{\frac{4D^{2}H}{dh^{2}} , \frac{4}{h}D + \frac{2}{3}C_{1}, \frac{6}{h}D\}

choosing the constant C_{1} appearing in Theorem 3 large enough for all cubic sections.

Now let K be the compact convex set containing all vectors M = (M_{0}, ..., M_{n+1})^{T} with

c ≤ |M_{j}| ≤ C for j ∈ R,
sgn M_{j} = sgn D_{j} for j ∈ R,
|M_{j}| ≤ C_{1} + \frac{1}{2} C for j ∉ R.

The iteration mapping T : K → K to be defined will combine a Gauss–Seidel iteration of M_{j} on rational points with the standard solution method on cubic sections using second derivative local boundary values on the rational boundary points:

Case a) Let x_{j} be a rational point between two rational intervals or on the boundary (in case of a boundary condition of type α). Then we set in accordance with (5)

\begin{matrix}T_{j}(M)^{2/3} := & 2D_{j}(h_{j} M_{j+1}^{1/3} + h_{j−1}M_{j−1}^{1/3})^{−1} \\ sgn T_{j}(M) := & sgn D_{j}. \end{matrix}

Case b) Let [x_{j−1}, x_{j}] be a cubic, [x_{j}, x_{j+1}] be a rational interval, respectively. Then we use (6) to define T_{j}(M) as the (unique) sign–correct solution of

\frac{1}{6} h_{j−1}M_{j−1} + \frac{1}{3} h_{j−1}T_{j}(M) + \frac{1}{2} h_{j}M_{j+1}^{1/3}T_{j}(M)^{2/3} = D_{j}

(for details see below).

Case c) On non–rational points x_{k} of cubic sections we define T_{k}(M) := s′′(x_{k}) for the cubic spline interpolating the given data and the values M_{i} on rational boundary points x_{i} of the section.

In case a) the sign of T_{j}(M) will be correct since D_{j}, M_{j−1}, M_{j+1} are all of the same sign. The nondegeneracy condition for cubic sections guarantees that in case b) we have

D_{j}M_{j−1} ≤ 0, sgn D_{j} = sgn M_{j+1} ≠ 0.

Then (T_{j}(M))^{1/3} solves an equation of the form

x^{3} + b_{2}x^{2} = b_{0} = D_{j} − \frac{1}{6} h_{j−1}M_{j−1} ≠ 0  (12)

with sgn b_{0} = sgn b_{2} and should have the sign

sgn T_{j}(M) = sgn D_{j} = sgn b_{0} = sgn b_{2}.

Since b_{0}x+b_{0}b_{2} intersects b_{0}^{2}/x^{2} in exactly one point x with sgn x = sgn b_{0}, the cubic equation (12) has precisely one (simple) sign–correct solution depending smoothly on the data. At this point the effect of nondegeneracy of cubic sections is clearly visible: if D_{j} − \frac{1}{6} h_{j−1}M_{j−1} has the wrong sign, T_{j}(M) will get the wrong sign, too.

Now T is well–defined and continuous on K; we still have to prove that T (K) ⊆ K holds.

In case a) we have

dH^{−1}C^{−1/3} ≤ T_{j}(M)^{2/3} ≤ 2 D_{j}h_{j}^{−1}M_{j+1}^{−1/3} ≤ 2 Dh^{−1}c^{−1/3},

where we allowed a possible boundary point and assumed j < n without loss of generality. Then

c^{2} = d^{3}H^{−3}C^{−1} ≤ T_{j}(M)^{2} ≤ 8 D^{3}h^{−3}c^{−1} = C^{2},

as required. Case c) implies

|T_{k}(M)| ≤ C_{1} + \frac{1}{2} C for k ∉ R

via Theorem (3.2) and we are left with case b). To get upper bounds we rewrite T_{j}(M) as

T_{j}(M) = \frac{\frac{3}{h_{j−1}} D_{j} − \frac{1}{2} M_{j−1}}{1 + \frac{3h_{j}}{2h_{j−1}} (\frac{M_{j+1}}{T_{j}(M)})^{1/3}}

and get

|T_{j}(M)| ≤ \frac{3}{h} D + \frac{1}{2} |M_{j−1}| ≤ C.

If j − 1 ∉ R, then

\begin{aligned}|T_{j}(M)| & ≤ \frac{3}{h} D + \frac{1}{2} (C_{1} + \frac{1}{2} C) \\ & ≤ \frac{3}{h} D + \frac{1}{2} C_{1} + \frac{1}{4} C ≤ C, \end{aligned}

while for j − 1 ∈ R we conclude

|T_{j}(M)| ≤ \frac{3}{h} D + \frac{1}{2} C ≤ C.

To estimate |T_{j}(M)| from below we can assume |T_{j}(M)| ≤ 1 and write the equation as

\frac{h_{j−1}}{3} |T_{j}(M)| + \frac{h_{j}}{2} |M_{j+1}^{1/3}| |T_{j}(M)|^{2/3} = |D_{j} − \frac{h_{j−1}}{2} M_{j−1}| ≥ d,

taking the sign distribution into account. Then

\begin{aligned}d & ≤ \frac{H}{3} |T_{j}(M)| + \frac{H}{2} C^{1/3}T_{j}(M)^{2/3} \\ & ≤ (\frac{H}{3} + \frac{H}{2} C^{1/3})T_{j}(M)^{2/3}, \end{aligned}

\begin{aligned}T_{j}(M)^{2/3} & ≥ dH^{−1} (\frac{1}{3} + \frac{1}{2} C^{1/3})^{−1} ≥ dH^{−1}C^{−1/3} (\frac{1}{2} + \frac{1}{3} C^{−1/3})^{−1} \\ & ≥ dH^{−1}C^{−1/3}, \end{aligned}

and finally

T_{j}^{2}(M) ≥ d^{3}H^{−3}C^{−1} = c^{2}.

Uniqueness follows from a simple argument counting the zeros of the difference of two solutions (see [3]).

Remark. Under the hypothesis (11) of Theorem 3.1 a pure Gauss–Seidel iteration can be employed. This will streamline the proof of Theorem 4.1. However, the following example will show that (11) is not necessary for nondegeneracy.

Example. Consider a symmetric problem with n = 9, x_{i} = i − 5, 0 ≤ i ≤ 10, and data given in condensed form as D_{0}, D_{1}, 0, D_{3}, 0, D_{5}, 0, D_{7} = D_{3}, 0, D_{9} = D_{1}, D_{10} = D_{0} with sgn D_{0} = sgn D_{1} = 1. If a solution exists, it is unique (see [3]) and therefore it must be symmetric. The data for the cubic section [x_{1}, x_{9}] satisfy (11) in case D_{3} ≥ 0, D_{5} ≥ 0, but not if D_{3} < 0, D_{5} > 0. Direct elimination of unknowns M_{i} = s′′(x_{i}) yields the equation

−97M_{2} = 26M_{1} + 42D_{3} + 3D_{5}

and nongeneracy occurs if and only if

42D_{3} + 3D_{5} ≥ 0,

because then any choice of M_{1} > 0 implies M_{2} < 0. This shows that the condition (11) is not necessary for nondegeneracy.

Furthermore, nondegeneracy is not necessary for solvability of the interpolation problem, because the final equation for M_{4} is

\frac{168}{97}M_{4} + 3(2D_{0})^{1/6}\sqrt{M_{4}} = 6D_{1} + \frac{42}{97}D_{3} + \frac{3}{97}D_{5},

having a positive solution, iff

582D_{1} + 42D_{3} + 3D_{5} > 0.

5 Numerical methods

The iteration defined in the previous section can be used to calculate the solution numerically, but we propose to use a stabilized Newton–Raphson method for reasons of computational efficiency. Since results of [3] yield uniqueness of the interpolating adaptive spline, the nonlinear system to be solved does not have other solutions in the domain K. We have used Newton’s method with a simple stepsize control assuring that

• iterates stay in K and

• ∥F (M)∥_{2}^{2} decreases, if the equations are written as F (M) = 0.

The computational cost of solving the adaptive rational spline interpolation problem was observed to stay proportional to n. Thus it does not exceed the computational cost of cubic spline interpolation by more than a constant factor.

Each Newton step involves a tridiagonal system of linear equations. For single precision and the starting values

M_{j} := \frac{2}{h_{j−1} + h_{j}} D_{j} = s′′(ξ_{j}) ≈ s′′(x_{j})

normally 5 – 6 Newton iterations are sufficient. Stabilization is required for small n and near– singular problems, but even in those cases 15 iterations are rarely exceeded. The values of a cubic spline solution can not be used to start the Newton iteration, because they may have wrong signs due to the fact that cubic splines do not preserve convexity.

6 Numerical examples

Figures 1, 2, and 3 show plots of interpolating adaptive rational splines (dashed lines) to certain functions (dotted lines) over the interval [0, 1]. The cubic interpolating splines (and the cubic pieces within adaptive rational splines) are represented by solid lines.

The convexity preserving feature of the adaptive rational spline cuts off the parasitic wiggles of the cubic spline near singularities or at the foot of sharp peaks. For large numbers of knots it shows the same asymptotics as cubic splines (see [3] and [6] for the convergence of rational splines to cubic splines), but it is inferior to the cubic spline by a certain factor (about 10 in most examples). The adaptive rational spline is superior to the cubic spline for small numbers of knots and if the problem has nearby singularities or sharp bends.

7 Two–dimensional rational splines

On a rectangular grid defined by the partitions (1) and

c = z_{0} < z_{1} < ... < z_{n+1} = d m ≥ 0

we construct a function s ∈ C^{2}([a, b] × [c, d]) which interpolates the tensor product data

\begin{matrix}y_{ij} = f(x_{i}, z_{j}), & 0 ≤ i ≤ n + 1, & 0 ≤ j ≤ m + 1, \\ y′_{ij} = \frac{∂f}{∂x} (x_{i}, z_{j}), & i = 0, i = n + 1, & 0 ≤ j ≤ m + 1, \\ z′_{ij} = \frac{∂f}{∂z} (x_{i}, z_{j}), & 0 ≤ i ≤ n + 1, & j = 0, j = m + 1, \\ y′′_{ij} = \frac{∂^{2}f}{∂x∂z} (x_{i}, z_{j}), & i = 0, i = n + 1, & j = 0, j = m + 1. \end{matrix}

Though the formalism of tensor product interpolation fails here due to the nonlinearity, we want to use the one–dimensional rational interpolation process wherever possible.

A straighforward method to do this is the following:

1) Interpolate data y_{ij} at (x_{i}, z_{j}), 0 ≤ j ≤ m + 1, for each fixed value x_{i}, 0 ≤ i ≤ n + 1, using the one–dimensional process.

2) Interpolate derivatives y′′_{i0}, y′_{i0}, y′_{i1},…, y′_{i,m+1}, y′′_{i,m+1} for i = 0 and i = m + 1 using the one–dimensional process.

3) For each of these interpolations, take z–derivatives at grid points and store these only.

4) Then for a fixed z ∈ [z_{j}, z_{j+1}] ⊆ [c, d] the values w_{i} at (x_{i}, z) for interpolation along the line z = const can be calculated easily from the stored Hermite data of step 1 at (x_{i}, z_{j}) and (x_{i}, z_{j+1}). This Hermite interpolation is possible by solving equations of the form (3) for M_{j} and M_{j+1}.

5) The same local Hermite interpolation is applied for \frac{∂}{∂x}–values at (x_{i}, z) for i = 0 and i = m + 1, but using the derivative data of step 2.

6) Application of the one–dimensional process to the data generated by steps 4 and 5 for fixed z yields a C^{2} function s_{z}(x) which interpolates at (x_{i}, z) for 0 ≤ i ≤ n + 1.

7) Repeat steps 4 – 6 for the required number of z–values to evaluate the function s(x, z) := s_{z}(x) wherever needed. Each z–value requires the solution of a nonlinear one–dimensional interpolation problem in step 6. Therefore it is highly recommended to proceed along the lines with z = const.

This process leads to an interpolant s(x, z) with the following properties:

• If all one–dimensional interpolants are purely rational, the solution s_{z}(x) ∈ C^{2}[a, b] depends smoothly on the data which vary twice continuously differentiable with z. Then we have s(x, z) ∈ C^{2}([a, b] × [c, d]).

• If all one–dimensional interpolants are purely rational, s(x, z) = s_{z}(x) is a purely concave/convex function of x for each z.

• If adaptive rational splines are used for local one–dimensional interpolation, there will be small jumps in the solution where the type of the interpolant changes. These are eliminated when the type of interpolation does not vary for successive applications of step 6. This means that the local type of interpolation with respect to x for fixed values of z must be independent of z.

• So far no local representation of s(x, z) is known.

• Consequently, it is an open question whether s reproduces convexity/concavity.

The main difference to tensor product interpolation with bicubic splines is that here the tensor product structure is confined to the data and not shared by the interpolating family of functions. Consequently, successive rational interpolation processes along x and z grid lines do not commute and there is no linear local superposition of the one–dimensional interpolation processes, allowing one to construct values of the interpolant s together with \frac{∂s}{∂x}, \frac{∂s}{∂z}, and \frac{∂^{2}s}{∂x∂z} at each grid point for a subsequent local Hermite–type interpolation on each subrectangle. When interpolation is carried out along z first, the construction of the two–dimensional rational interpolant at a certain point requires calculation of a non–local interpolant along a complete x–line. If the interpolant has to be evaluated successively along complete x– or z–lines, the nonlinear scheme is not much inferior to the linear bicubic tensor product scheme.

Figure 4 shows the rational interpolant to the major part of a half–sphere over [0, 1]^{2} on a 5 × 5–grid, while the nonconvex cubic interpolant appears in Figure 5. Furthermore, a sharp (Runge–type rational) peak is interpolated by an adaptive rational spline on an 11 × 11–grid (see Fig. 6), and Fig. 7 shows the corresponding bicubic solution. No precautions to ensure C^{2} differentiability were taken for Fig. 6.

References

[1] D. Braess, H. Werner, Tschebyscheff–Approximation mit einer Klasse rationaler Spline–Funktionen II, J. Approx. Theory 10, 1974, pp. 379–399.

[2] R. Schaback, Spezielle rationale Splinefunktionen , J. Approx. Theory 7, No. 3, 1973, pp. 281–292.

[3] R. Schaback, Interpolation mit nichtlinearen Klassen von Spline–Funktionen, J. Approx. Theory 8, No. 2, 1973, pp. 173–188.

[4] R. Schaback, Calculation of best approximations by rational splines, in: G.G. Lorentz, C.K. Chui, and L.L. Schumaker (eds): Approximation Theory II, Academic Press, 1976, pp. 533–547.

[5] H. Späth, Spline–Algorithmen zur Konstruktion glatter Kurven und Flächen, 2. Auflage, Oldenbourg–Verlag, 1978.

[6] H. Werner, Interpolation and Integration of Initial Value Problems of Ordinary Differential Equations by Regular Splines, SIAM J. Numer. Anal. 12 (1975), pp. 255–271.

[7] H. Werner, H. Loeb, Tschebyscheff Approximation by Regular Splines with Free Knots, in: Lecture Notes in Mathematics, Bd. 556, Springer–Verlag, 1976, pp. 439–452.

Fig. 1. f(x) = arctan(100x − 50), 8 knots. Fig. 2. f(x) = (100x^{2} − 100x + 26)^{−1}, 6 knots.

Fig. 3. f(x) = −0.0001 + \sqrt{0.5001^{2} − (x − 0.5)^{2}}, 11 knots. Fig. 4. Purely rational bivariate spline interpolant, 5 × 5–grid.

Fig. 5. Bicubic spline interpolant, same data as Fig. 4. Fig. 6. Adaptive rational spline interpolant, 11 × 11–grid. Fig. 7. Bicubic spline interpolant, same data as Fig. 6.

Author’s address : Prof. Dr. R. Schaback

Institut für Numerische und Angewandte Mathematik

Georg–August–Universität Göttingen Lotzestraße 16–18 D–3400 Göttingen Fed. Rep. of Germany