I have a function evaluated at predefined grid points
$x_i=\frac{i}{N},\ i=0,...,N,\ f(x_i)=y_i$
I am looking for an interpolation scheme satisfying
$\|I(f)-f\|_{L^\infty} = \mathcal{O}((1/N)^3)\qquad \qquad (1)$
for smooth functions $f$. The interpolation should only be dependend on the values $y_i$. This can of course be accomplished by using (piecewise) quadratic splines, i.e.
$I(f)(x) = y_i(1-\left(\frac{x-x_i}{N}\right)^2) + y_{i-1}(\left(\frac{x-x_i}{N}\right)^2/2+\frac{x-x_i}{2N})+y_{i+1}(\left(\frac{x-x_i}{N}\right)^2/2-\frac{x-x_i}{2N}),\ x\in[x_i,x_{i+1}] $
However, those quadratic splines use negative weights. Therefore the following condition is not satisfied.
$\|I(f)\|_{L^\infty} \leq \|f\|_{L^\infty}\qquad \qquad (2)$
Is there a way to perform an interpolation which satisfies $(1)$ and only uses positive weights? A nonlinear dependence on the $y_i$ would probably work for me too ad long as $(1)$ and $(2)$ hold true.