How to find the best interpolating function if we know $y(x_i)$ and $dy(x_i)/dx$

Imagine you are given a set of data points $\{x_i,y_i\}$, supplemented by a list of known first derivatives $\{y'_i\}$.

How would you construct an interpolating function $y(x)$ (which satisfies $y(x_i)=y_i$ and $y'(x_i)=y'_i$), such that the derivative, calculated from this function has the smallest error.

The function is expected to be well behaved, consisting of piece-wise power functions, smoothely linked.

What would you do if you knew second derivatives as well and wanted to work with them instead of the first ones?

-
I would first try a cubic spline. – Chris Taylor Feb 27 '12 at 15:22
I like cubic splines, but 1) This isn't extendable to the case when you know second derivatives as well as the first ones, 2) Even in the simpler case of having only first derivatives, if you want to work with y', y'' shall be non-smooth, which isn't a good feature. – Alexey Bobrick Feb 27 '12 at 15:33
Like a Hermite Interpolation problem? – mixedmath Feb 27 '12 at 16:01
The minimum degree polynomial interpolation can be found in "Newton form" easily using "divided differences". This lets you replace cubic splines with quartic splines, etc. The method easily allows you to have 5 derivatives at one point given and only 2 at another, and none (only the function) at another (per spline). To do even better than this: I think the more clever interpolation ideas require the algorithm to choose the xi, and the user to supply the yi and yi' etc. – Jack Schmidt Feb 27 '12 at 16:03
@Jack Schmidt, when going to higher order splines, one would expect them to be unstable and not really recovering the original function. Or would you suggest power functions being well recovered by, say a 5-th order spline for the case of 2-nd derivatives given? – Alexey Bobrick Feb 27 '12 at 22:04

To settle this question: as already mentioned in the comments: one of the simplest approaches to this problem is (piecewise) Hermite interpolation (of which splines are a special case, since they have extra restrictions that guarantee continuity of their higher derivatives). In particular, there is an easy modification of the usual divided difference scheme to produce the unique polynomial that satisfies the function and derivative values supplied.

On the other hand, (polynomial) Hermite interpolation is not the be all and end all of interpolating given function values and derivatives. For instance, there are circumstances where a rational Hermite interpolant (a rational function that satisfies a given set of function and derivative values) might be more appropriate, when the function being approximated is expected to have poles or similar behavior, and if the function whose values are given is expected to be periodic, then certainly one needs to look into trigonometric Hermite interpolation. The computations needed for these variants are slightly more intricate than in the polynomial case, but if they give better results ("better" being up to your criteria of course), then the additional bother is certainly justifiable.

-
Thank you, sir, especially for the good references to the theory! – Alexey Bobrick Mar 19 '12 at 15:01

This answer is a pitch for the 'Weak Solution Techniques' especially intended for the case where your sampled data might have noise (like in the real world). You could choose just about any basis: $[\phi_k(x)| k \in 0:N]$ (one of my favorites is B-splines). The great part about such a basis is that you know the derivatives of all the basis vectors in advance.

In this approach you would represent $y(x) \approx f(x;a) = \sum_{k=0}^{N} a_k \phi_k(x)$ and seek to find the best coefficients $a_k$ in the sense that they minimize an energy functional. I usually take a variational approach to this problem - this way we can accommodate any noise that might be in the data $y_k$. In particular, we can know exactly how well we are are approximating the data by setting up a cost functional such as:

$$J(a) = ||y-\sum_{k=0}^{N} a_k \phi_k(x) ||_2^2 + \beta_1||\frac{dy}{dx}-\sum_{k=0}^{N} a_k \frac{d\phi_k}{dx} ||_2^2 + \cdots$$

the dots are because you could keep on adding penalties for each kind of information that you have. At this point you need to enforce the necessary condition for extrema: namely, for each $k$, $$\frac{\partial J(a)}{\partial a_k} = 0$$

Lets call the matrix of sampled basis functions $[\phi_0(x_k), \phi_1(x_k), \phi_2(x_k),\cdots ] =\Phi$ and similarly $[\frac{d\phi_0}{dx},\frac{d\phi_1}{dx},\frac{d\phi_2}{dx},\cdots ] = d\Phi$ then the necessary condition yields the Least Squares style equation:

$(\Phi^*\Phi+\beta_1d\Phi^*d\Phi)a = \Phi^*y + \beta_1d\Phi^*\frac{dy}{dx}$ where $\Phi^*$ and $d\Phi^*$ are the adjoints and $y$ and $dy/dx$ are the sampled data.

There is a lot of literature on how to choose the regularization parameter $\beta_1$ as well as other kinds of penalties (Tikhonov, $L_1$, TV) that you might add (as Bayesian priors)

-
Thank you! Indeed, rather simple and natural an approach. Do you know, what ideas are there behind when chosing $\beta_1$ and alike? – Alexey Bobrick May 15 '14 at 19:47
In the case of a two-norm energy penalty (Tikhonov regularization - also called a Wiener filter) one typically uses the ratio of the energy in the noise to the energy in the signal (the square of the inverse of the signal to noise ratio). A very good reference on this is Curt Vogel's book on inverse problems. He also covers techniques such as Ordinary Cross Validation and Generalized Cross Validation and the L-curve method. You can also get a bit of this idea from chapter 4 of Luenberger's 'Optimization by vector space methods' – eigenjohnson May 15 '14 at 23:55
Big thank you! Very nice of you to share these references! – Alexey Bobrick May 16 '14 at 14:38