|
| 1 | +# Interpolation |
| 2 | + |
| 3 | +The ***Numerics*** library provides interpolation methods for estimating values between known data points, from simple linear interpolation to smooth spline methods. |
| 4 | + |
| 5 | +## Overview |
| 6 | + |
| 7 | +| Method | Continuity | Best For | |
| 8 | +|--------|------------|----------| |
| 9 | +| Linear | C⁰ | Simple, fast | |
| 10 | +| Cubic Spline | C² | Smooth curves | |
| 11 | +| Akima Spline | C¹ | Data with outliers | |
| 12 | +| Polynomial | C∞ | Small datasets | |
| 13 | +| Bilinear | C⁰ | 2D gridded data | |
| 14 | +| Bicubic | C¹ | Smooth 2D surfaces | |
| 15 | + |
| 16 | +--- |
| 17 | + |
| 18 | +## One-Dimensional Interpolation |
| 19 | + |
| 20 | +### Linear Interpolation |
| 21 | + |
| 22 | +Connects points with straight line segments [[1]](#ref1): |
| 23 | + |
| 24 | +```math |
| 25 | +y = y_i + \frac{y_{i+1} - y_i}{x_{i+1} - x_i}(x - x_i) |
| 26 | +``` |
| 27 | + |
| 28 | +```cs |
| 29 | +using Numerics.Data.Interpolation; |
| 30 | + |
| 31 | +double[] x = { 0, 1, 2, 3, 4, 5 }; |
| 32 | +double[] y = { 0, 1, 4, 9, 16, 25 }; |
| 33 | + |
| 34 | +var linear = new LinearInterpolation(x, y); |
| 35 | + |
| 36 | +// Interpolate at a single point |
| 37 | +double value = linear.Interpolate(2.5); |
| 38 | +Console.WriteLine($"f(2.5) ≈ {value:F4}"); |
| 39 | + |
| 40 | +// Interpolate at multiple points |
| 41 | +double[] xNew = { 0.5, 1.5, 2.5, 3.5, 4.5 }; |
| 42 | +double[] yNew = linear.Interpolate(xNew); |
| 43 | +``` |
| 44 | + |
| 45 | +**Properties**: |
| 46 | +- Continuous (C⁰) but not smooth (discontinuous first derivative) |
| 47 | +- Never overshoots data |
| 48 | +- Fast computation |
| 49 | + |
| 50 | +### Cubic Spline Interpolation |
| 51 | + |
| 52 | +Fits piecewise cubic polynomials with continuous second derivatives [[1]](#ref1): |
| 53 | + |
| 54 | +```cs |
| 55 | +var spline = new CubicSplineInterpolation(x, y); |
| 56 | + |
| 57 | +double value = spline.Interpolate(2.5); |
| 58 | +Console.WriteLine($"f(2.5) ≈ {value:F4}"); |
| 59 | + |
| 60 | +// Get derivative at a point |
| 61 | +double derivative = spline.Derivative(2.5); |
| 62 | +Console.WriteLine($"f'(2.5) ≈ {derivative:F4}"); |
| 63 | + |
| 64 | +// Get second derivative |
| 65 | +double secondDeriv = spline.SecondDerivative(2.5); |
| 66 | +``` |
| 67 | + |
| 68 | +**Boundary Conditions**: |
| 69 | + |
| 70 | +```cs |
| 71 | +// Natural spline (second derivative = 0 at endpoints) |
| 72 | +var natural = new CubicSplineInterpolation(x, y, BoundaryCondition.Natural); |
| 73 | + |
| 74 | +// Clamped spline (specified first derivatives at endpoints) |
| 75 | +var clamped = new CubicSplineInterpolation(x, y, BoundaryCondition.Clamped, |
| 76 | + leftDerivative: 0.0, rightDerivative: 10.0); |
| 77 | + |
| 78 | +// Not-a-knot (default, continuous third derivative at second and second-to-last points) |
| 79 | +var notAKnot = new CubicSplineInterpolation(x, y, BoundaryCondition.NotAKnot); |
| 80 | +``` |
| 81 | + |
| 82 | +**Properties**: |
| 83 | +- C² continuous (smooth first and second derivatives) |
| 84 | +- Can overshoot/undershoot data |
| 85 | +- Good for smooth underlying functions |
| 86 | + |
| 87 | +### Akima Spline |
| 88 | + |
| 89 | +Locally-determined spline that reduces overshoot [[2]](#ref2): |
| 90 | + |
| 91 | +```cs |
| 92 | +var akima = new AkimaSplineInterpolation(x, y); |
| 93 | + |
| 94 | +double value = akima.Interpolate(2.5); |
| 95 | +``` |
| 96 | + |
| 97 | +**Properties**: |
| 98 | +- C¹ continuous (smooth first derivative only) |
| 99 | +- Less oscillation than cubic spline |
| 100 | +- Better for data with outliers or rapid changes |
| 101 | + |
| 102 | +### Polynomial Interpolation |
| 103 | + |
| 104 | +Fits a single polynomial through all points [[1]](#ref1): |
| 105 | + |
| 106 | +```cs |
| 107 | +var poly = new PolynomialInterpolation(x, y); |
| 108 | + |
| 109 | +double value = poly.Interpolate(2.5); |
| 110 | + |
| 111 | +// Get the polynomial coefficients |
| 112 | +double[] coeffs = poly.Coefficients; |
| 113 | +``` |
| 114 | + |
| 115 | +**Warning**: High-degree polynomials can oscillate wildly (Runge's phenomenon). Use only for small datasets (n ≤ 10). |
| 116 | + |
| 117 | +### Monotonic Interpolation |
| 118 | + |
| 119 | +Preserves monotonicity of the data (no spurious extrema): |
| 120 | + |
| 121 | +```cs |
| 122 | +var monotonic = new MonotonicInterpolation(x, y); |
| 123 | + |
| 124 | +double value = monotonic.Interpolate(2.5); |
| 125 | +``` |
| 126 | + |
| 127 | +**Properties**: |
| 128 | +- Guarantees interpolant is monotonic if data is monotonic |
| 129 | +- Useful for cumulative distribution functions |
| 130 | +- C¹ continuous |
| 131 | + |
| 132 | +--- |
| 133 | + |
| 134 | +## Extrapolation |
| 135 | + |
| 136 | +By default, interpolation methods only work within the data range. To extrapolate: |
| 137 | + |
| 138 | +```cs |
| 139 | +var linear = new LinearInterpolation(x, y); |
| 140 | +linear.AllowExtrapolation = true; |
| 141 | + |
| 142 | +// Now works outside [x_min, x_max] |
| 143 | +double extrapolated = linear.Interpolate(6.0); |
| 144 | +``` |
| 145 | + |
| 146 | +**Caution**: Extrapolation can be unreliable, especially for splines. |
| 147 | + |
| 148 | +--- |
| 149 | + |
| 150 | +## Two-Dimensional Interpolation |
| 151 | + |
| 152 | +### Bilinear Interpolation |
| 153 | + |
| 154 | +Linear interpolation on a 2D grid [[1]](#ref1): |
| 155 | + |
| 156 | +```cs |
| 157 | +double[] xGrid = { 0, 1, 2 }; |
| 158 | +double[] yGrid = { 0, 1, 2 }; |
| 159 | +double[,] zValues = { |
| 160 | + { 0, 1, 4 }, |
| 161 | + { 1, 2, 5 }, |
| 162 | + { 4, 5, 8 } |
| 163 | +}; |
| 164 | + |
| 165 | +var bilinear = new BilinearInterpolation(xGrid, yGrid, zValues); |
| 166 | + |
| 167 | +double value = bilinear.Interpolate(0.5, 0.5); |
| 168 | +Console.WriteLine($"f(0.5, 0.5) ≈ {value:F4}"); |
| 169 | +``` |
| 170 | + |
| 171 | +**Algorithm**: |
| 172 | +1. Linear interpolation in x-direction at two y-values |
| 173 | +2. Linear interpolation in y-direction between results |
| 174 | + |
| 175 | +### Bicubic Interpolation |
| 176 | + |
| 177 | +Cubic interpolation on a 2D grid for smoother surfaces: |
| 178 | + |
| 179 | +```cs |
| 180 | +var bicubic = new BicubicInterpolation(xGrid, yGrid, zValues); |
| 181 | + |
| 182 | +double value = bicubic.Interpolate(0.5, 0.5); |
| 183 | + |
| 184 | +// Partial derivatives |
| 185 | +double dfdx = bicubic.PartialDerivativeX(0.5, 0.5); |
| 186 | +double dfdy = bicubic.PartialDerivativeY(0.5, 0.5); |
| 187 | +``` |
| 188 | + |
| 189 | +--- |
| 190 | + |
| 191 | +## Interpolation from Irregular Data |
| 192 | + |
| 193 | +### Inverse Distance Weighting (IDW) |
| 194 | + |
| 195 | +For scattered (non-gridded) 2D data: |
| 196 | + |
| 197 | +```cs |
| 198 | +double[] xPoints = { 0, 1, 2, 0, 1, 2 }; |
| 199 | +double[] yPoints = { 0, 0, 0, 1, 1, 1 }; |
| 200 | +double[] zPoints = { 0, 1, 4, 1, 2, 5 }; |
| 201 | + |
| 202 | +var idw = new InverseDistanceWeighting(xPoints, yPoints, zPoints); |
| 203 | +idw.Power = 2; // Inverse square weighting |
| 204 | +
|
| 205 | +double value = idw.Interpolate(0.5, 0.5); |
| 206 | +``` |
| 207 | + |
| 208 | +**Formula**: |
| 209 | +```math |
| 210 | +z = \frac{\sum_{i=1}^{n} w_i z_i}{\sum_{i=1}^{n} w_i}, \quad w_i = \frac{1}{d_i^p} |
| 211 | +``` |
| 212 | + |
| 213 | +where $d_i$ is the distance to point $i$ and $p$ is the power parameter. |
| 214 | + |
| 215 | +--- |
| 216 | + |
| 217 | +## Applications |
| 218 | + |
| 219 | +### Resampling Time Series |
| 220 | + |
| 221 | +```cs |
| 222 | +// Original irregular time series |
| 223 | +double[] times = { 0, 1.5, 3.2, 4.0, 6.1, 8.0 }; |
| 224 | +double[] values = { 10, 15, 12, 18, 14, 20 }; |
| 225 | + |
| 226 | +var spline = new CubicSplineInterpolation(times, values); |
| 227 | + |
| 228 | +// Resample to regular intervals |
| 229 | +double dt = 0.5; |
| 230 | +int n = (int)((times.Last() - times.First()) / dt) + 1; |
| 231 | +double[] regularTimes = new double[n]; |
| 232 | +double[] regularValues = new double[n]; |
| 233 | + |
| 234 | +for (int i = 0; i < n; i++) |
| 235 | +{ |
| 236 | + regularTimes[i] = times.First() + i * dt; |
| 237 | + regularValues[i] = spline.Interpolate(regularTimes[i]); |
| 238 | +} |
| 239 | +``` |
| 240 | + |
| 241 | +### Rating Curve Interpolation |
| 242 | + |
| 243 | +```cs |
| 244 | +// Stage-discharge rating curve |
| 245 | +double[] stage = { 0, 1, 2, 3, 4, 5 }; // feet |
| 246 | +double[] discharge = { 0, 50, 200, 500, 1000, 1800 }; // cfs |
| 247 | +
|
| 248 | +// Use monotonic interpolation to ensure increasing discharge |
| 249 | +var rating = new MonotonicInterpolation(stage, discharge); |
| 250 | + |
| 251 | +// Convert stage reading to discharge |
| 252 | +double observedStage = 2.5; |
| 253 | +double estimatedQ = rating.Interpolate(observedStage); |
| 254 | +Console.WriteLine($"Stage {observedStage} ft → {estimatedQ:F0} cfs"); |
| 255 | +``` |
| 256 | + |
| 257 | +### Probability Plotting |
| 258 | + |
| 259 | +```cs |
| 260 | +// Empirical CDF interpolation |
| 261 | +double[] sortedData = data.OrderBy(x => x).ToArray(); |
| 262 | +int n = sortedData.Length; |
| 263 | +double[] plottingPositions = new double[n]; |
| 264 | + |
| 265 | +for (int i = 0; i < n; i++) |
| 266 | + plottingPositions[i] = (i + 1 - 0.4) / (n + 0.2); // Cunnane |
| 267 | +
|
| 268 | +// Interpolate to find quantiles |
| 269 | +var empiricalCDF = new MonotonicInterpolation(sortedData, plottingPositions); |
| 270 | +var inverseCDF = new MonotonicInterpolation(plottingPositions, sortedData); |
| 271 | + |
| 272 | +// Find value at 95th percentile |
| 273 | +double q95 = inverseCDF.Interpolate(0.95); |
| 274 | +Console.WriteLine($"95th percentile: {q95:F2}"); |
| 275 | +``` |
| 276 | + |
| 277 | +### DEM Interpolation |
| 278 | + |
| 279 | +```cs |
| 280 | +// Digital Elevation Model (gridded) |
| 281 | +double[] eastings = { 0, 100, 200, 300, 400 }; |
| 282 | +double[] northings = { 0, 100, 200, 300, 400 }; |
| 283 | +double[,] elevations = LoadDEM(); // 5x5 grid |
| 284 | +
|
| 285 | +var dem = new BicubicInterpolation(eastings, northings, elevations); |
| 286 | + |
| 287 | +// Get elevation at arbitrary point |
| 288 | +double x = 150, y = 175; |
| 289 | +double elevation = dem.Interpolate(x, y); |
| 290 | +Console.WriteLine($"Elevation at ({x}, {y}): {elevation:F1} m"); |
| 291 | + |
| 292 | +// Compute slope |
| 293 | +double dzdx = dem.PartialDerivativeX(x, y); |
| 294 | +double dzdy = dem.PartialDerivativeY(x, y); |
| 295 | +double slope = Math.Atan(Math.Sqrt(dzdx * dzdx + dzdy * dzdy)) * 180 / Math.PI; |
| 296 | +Console.WriteLine($"Slope: {slope:F1}°"); |
| 297 | +``` |
| 298 | + |
| 299 | +--- |
| 300 | + |
| 301 | +## Choosing an Interpolation Method |
| 302 | + |
| 303 | +| Scenario | Recommended Method | |
| 304 | +|----------|-------------------| |
| 305 | +| Fast, simple | Linear | |
| 306 | +| Smooth curves | Cubic Spline | |
| 307 | +| Data with outliers | Akima Spline | |
| 308 | +| Must preserve monotonicity | Monotonic | |
| 309 | +| CDF interpolation | Monotonic | |
| 310 | +| Small dataset, exact fit | Polynomial | |
| 311 | +| 2D gridded data | Bilinear or Bicubic | |
| 312 | +| 2D scattered data | IDW | |
| 313 | + |
| 314 | +--- |
| 315 | + |
| 316 | +## Performance Considerations |
| 317 | + |
| 318 | +1. **Precompute coefficients**: Spline setup is O(n), evaluation is O(log n) |
| 319 | +2. **Binary search for interval**: All methods use binary search to find the containing interval |
| 320 | +3. **Caching**: Reuse interpolation objects for multiple queries |
| 321 | + |
| 322 | +```cs |
| 323 | +// Efficient: Create once, use many times |
| 324 | +var spline = new CubicSplineInterpolation(x, y); |
| 325 | + |
| 326 | +for (int i = 0; i < 1000000; i++) |
| 327 | +{ |
| 328 | + double val = spline.Interpolate(queryPoints[i]); |
| 329 | +} |
| 330 | + |
| 331 | +// Inefficient: Recreating each time |
| 332 | +for (int i = 0; i < 1000000; i++) |
| 333 | +{ |
| 334 | + var spline = new CubicSplineInterpolation(x, y); // Don't do this! |
| 335 | + double val = spline.Interpolate(queryPoints[i]); |
| 336 | +} |
| 337 | +``` |
| 338 | + |
| 339 | +--- |
| 340 | + |
| 341 | +## Comparison of Methods |
| 342 | + |
| 343 | +```cs |
| 344 | +double[] x = { 0, 1, 2, 3, 4, 5 }; |
| 345 | +double[] y = { 0, 0.8, 0.9, 0.1, -0.8, -1.0 }; |
| 346 | + |
| 347 | +var linear = new LinearInterpolation(x, y); |
| 348 | +var cubic = new CubicSplineInterpolation(x, y); |
| 349 | +var akima = new AkimaSplineInterpolation(x, y); |
| 350 | + |
| 351 | +Console.WriteLine("x Linear Cubic Akima"); |
| 352 | +Console.WriteLine("--- ------ ----- -----"); |
| 353 | + |
| 354 | +for (double xi = 0; xi <= 5; xi += 0.5) |
| 355 | +{ |
| 356 | + Console.WriteLine($"{xi:F1} {linear.Interpolate(xi),6:F3} {cubic.Interpolate(xi),6:F3} {akima.Interpolate(xi),6:F3}"); |
| 357 | +} |
| 358 | +``` |
| 359 | + |
| 360 | +--- |
| 361 | + |
| 362 | +## References |
| 363 | + |
| 364 | +<a id="ref1">[1]</a> Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. (2007). *Numerical Recipes: The Art of Scientific Computing* (3rd ed.). Cambridge University Press. |
| 365 | + |
| 366 | +<a id="ref2">[2]</a> Akima, H. (1970). A new method of interpolation and smooth curve fitting based on local procedures. *Journal of the ACM*, 17(4), 589-602. |
| 367 | + |
| 368 | +<a id="ref3">[3]</a> de Boor, C. (2001). *A Practical Guide to Splines* (Rev. ed.). Springer. |
0 commit comments