Monotone (Shape-Preserving) Interpolation
Monotone interpolation preserves the monotonicity of the data being interpolated, and when the data points are also monotonic, the slopes of the interpolant should also be monotonic.
spinterp(x, y, xp)
x, y |
x- and y-coordinates of the points that shall be interpolated. |
xp |
points that should be interpolated. |
This implementation follows a cubic version of the method of Delbourgo and Gregory. It yields ‘shaplier’ curves than the Stineman method.
The calculation of the slopes is according to recommended practice:
- monotonic and convex –> harmonic
- monotonic and nonconvex –> geometric
- nonmonotonic and convex –> arithmetic
- nonmonotonic and nonconvex –> circles (Stineman) [not implemented]
The choice of supplementary coefficients r[i] depends on whether
the data are montonic or convex or both:
- monotonic, but not convex
- otherwise
and that can be detected from the data. The choice r[i]=3 for all
i results in the standard cubic Hermitean rational interpolation.
The interpolated values at all the points of xp.
At the moment, the data need to be monotonic and the case of convexity is not considered.
Stan Wagon (2010). Mathematica in Action. Third Edition, Springer-Verlag.
stinepack::stinterp, demography::cm.interp
data1 <- list(x = c(1,2,3,5,6,8,9,11,12,14,15),
y = c(rep(10,6), 10.5,15,50,60,95))
data2 <- list(x = c(0,1,4,6.5,9,10),
y = c(10,4,2,1,3,10))
data3 <- list(x = c(7.99,8.09,8.19,8.7,9.2,10,12,15,20),
y = c(0,0.000027629,0.00437498,0.169183,0.469428,
0.94374,0.998636,0.999919,0.999994))
data4 <- list(x = c(22,22.5,22.6,22.7,22.8,22.9,
23,23.1,23.2,23.3,23.4,23.5,24),
y = c(523,543,550,557,565,575,
590,620,860,915,944,958,986))
data5 <- list(x = c(0,1.1,1.31,2.5,3.9,4.4,5.5,6,8,10.1),
y = c(10.1,8,4.7,4.0,3.48,3.3,5.8,7,7.7,8.6))
data6 <- list(x = c(-0.8, -0.75, -0.3, 0.2, 0.5),
y = c(-0.9, 0.3, 0.4, 0.5, 0.6))
data7 <- list(x = c(-1, -0.96, -0.88, -0.62, 0.13, 1),
y = c(-1, -0.4, 0.3, 0.78, 0.91, 1))
data8 <- list(x = c(-1, -2/3, -1/3, 0.0, 1/3, 2/3, 1),
y = c(-1, -(2/3)^3, -(1/3)^3, -(1/3)^3, (1/3)^3, (1/3)^3, 1))
## Not run:
opr <- par(mfrow=c(2,2))
# These are well-known test cases:
D <- data1
plot(D, ylim=c(0, 100)); grid()
xp <- seq(1, 15, len=51); yp <- spinterp(D$x, D$y, xp)
lines(spline(D), col="blue")
lines(xp, yp, col="red")
D <- data3
plot(D, ylim=c(0, 1.2)); grid()
xp <- seq(8, 20, len=51); yp <- spinterp(D$x, D$y, xp)
lines(spline(D), col="blue")
lines(xp, yp, col="red")
D <- data4
plot(D); grid()
xp <- seq(22, 24, len=51); yp <- spinterp(D$x, D$y, xp)
lines(spline(D), col="blue")
lines(xp, yp, col="red")
# Fix a horizontal slope at the end points
D <- data8
x <- c(-1.05, D$x, 1.05); y <- c(-1, D$y, 1)
plot(D); grid()
xp <- seq(-1, 1, len=101); yp <- spinterp(x, y, xp)
lines(spline(D, n=101), col="blue")
lines(xp, yp, col="red")
par(opr)
## End(Not run)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.