How can I find a non-negative interpolation function? In numerical mathematics I have learnt about some interpolation methods, however today I've come across some sort of interpolation problem which I don't know how to solve or even work with:
Let $(x_i, y_i) \in \mathbb{R}^2, i=0, \dots, n$ and $y_i \geq 0$ for all $0 \leq i \leq n$. I want to find a function (preferably a polynomial of arbitrary degree) $f \geq 0$ such that $f(x_i) = y_i$ for all $0 \leq i \leq n$.
How would I approach the general case? Is there a general solution? I know of interpolation techniques considering derivatives (Hermite) and also know of spline-interpolation, but those don't really help. Thanks a lot in advance.
 A: Positivity-preserving interpolation is hard, as searching for the italicized text will tell you. There are several things you could do, and the choice will depend on the practical context of your problem (not present in the question).
A cheap and dirty solution is to interpolate $(x_i,\sqrt{y_i})$ and square the interpolating function. If you interpolate with a polynomial, you'll have a nonnegative polynomial of twice the degree. The obvious drawback is that the squared function will have somewhat unnatural behavior in the regions where the interpolant was negative before squaring. Like on the right side of the graph here:  

Not to mention that interpolation by a high-degree polynomial is rarely desirable at all. 
If you use spline interpolation (e.g., cubic spline), the behavior of interpolant is much better; it takes a bit of effort to produce an example where the interpolant becomes negative. Here I started with the values 1 90 1 1 90 1 at equally spaced points, interpolated their square roots with a natural  cubic spline, and squared the spline:

The bump around 0.5 is unpleasant. It would be better to use logarithm instead of square root, but of course then we don't have a piecewise polynomial as a result: 

Still, in many cases the reason for data being positive is that it's naturally $\exp$ of something, so the above may be the best solution. 
A completely different approach is to turn to approximation instead of interpolation. Using positive compactly supported functions such as $B$-splines, you can approximate positive data by a positive function. 
A related question: Polynomial fitting where polynomial must be monotonically increasing
Related SE site: Computational Science
