Monotone cubic interpolation

From Wikipedia, the free encyclopedia
Jump to: navigation, search

In the mathematical subfield of numerical analysis, monotone cubic interpolation is a variant of cubic interpolation that preserves monotonicity of the data set being interpolated.

Monotonicity is preserved by linear interpolation but not guaranteed by cubic interpolation.

Monotone cubic Hermite interpolation[edit]

Example showing non-monotone cubic interpolation (in red) and monotone cubic interpolation (in blue) of a monotone data set.

Monotone interpolation can be accomplished using cubic Hermite spline with the tangents modified to ensure the monotonicity of the resulting Hermite spline.

An algorithm is also available for monotone quintic Hermite interpolation.

Interpolant selection[edit]

There are several ways of selecting interpolating tangents for each data point. This section will outline the use of the Fritsch–Carlson method.

Let the data points be for

  1. Compute the slopes of the secant lines between successive points:

    for .
  2. Initialize the tangents at every data point as the average of the secants,

    for ; if and have different sign, set . These may be updated in further steps. For the endpoints, use one-sided differences:

  3. For , if (if two successive are equal), then set as the spline connecting these points must be flat to preserve monotonicity. Ignore step 4 and 5 for those .
  4. Let and . If or are computed to be less than zero, then the input data points are not strictly monotone, and is a local extremum. In such cases, piecewise monotone curves can still be generated by choosing , although global strict monotonicity is not possible.
  5. To prevent overshoot and ensure monotonicity, at least one of the following conditions must be met:
    1. the function

      must have a value greater than or equal to zero;
    2. ; or
    3. .

If monotonicity must be strict then must have a value strictly greater than zero.

One simple way to satisfy this constraint is to restrict the magnitude of vector to a circle of radius 3. That is, if , then set and where .

Alternatively it is sufficient to restrict and . To accomplish this if , then set . Similarly for .

Note that only one pass of the algorithm is required.

Cubic interpolation[edit]

After the preprocessing, evaluation of the interpolated spline is equivalent to cubic Hermite spline, using the data , , and for .

To evaluate at , find the smallest value larger than , , and the largest value smaller than , , among such that . Calculate

and

then the interpolant is

where are the basis functions for the cubic Hermite spline.

Example implementation[edit]

The following JavaScript implementation takes a data set and produces a monotone cubic spline interpolant function:

/* Monotone cubic spline interpolation
   Usage example:
	var f = createInterpolant([0, 1, 2, 3, 4], [0, 1, 4, 9, 16]);
	var message = '';
	for (var x = 0; x <= 4; x += 0.5) {
		var xSquared = f(x);
		message += x + ' squared is about ' + xSquared + '\n';
	}
	alert(message);
*/
var createInterpolant = function(xs, ys) {
	var i, length = xs.length;
	
	// Deal with length issues
	if (length != ys.length) { throw 'Need an equal count of xs and ys.'; }
	if (length === 0) { return function(x) { return 0; }; }
	if (length === 1) {
		// Impl: Precomputing the result prevents problems if ys is mutated later and allows garbage collection of ys
		// Impl: Unary plus properly converts values to numbers
		var result = +ys[0];
		return function(x) { return result; };
	}
	
	// Rearrange xs and ys so that xs is sorted
	var indexes = [];
	for (i = 0; i < length; i++) { indexes.push(i); }
	indexes.sort(function(a, b) { return xs[a] < xs[b] ? -1 : 1; });
	var oldXs = xs, oldYs = ys;
	// Impl: Creating new arrays also prevents problems if the input arrays are mutated later
	xs = []; ys = [];
	// Impl: Unary plus properly converts values to numbers
	for (i = 0; i < length; i++) { xs.push(+oldXs[indexes[i]]); ys.push(+oldYs[indexes[i]]); }
	
	// Get consecutive differences and slopes
	var dys = [], dxs = [], ms = [];
	for (i = 0; i < length - 1; i++) {
		var dx = xs[i + 1] - xs[i], dy = ys[i + 1] - ys[i];
		dxs.push(dx); dys.push(dy); ms.push(dy/dx);
	}
	
	// Get degree-1 coefficients
	var c1s = [ms[0]];
	for (i = 0; i < dxs.length - 1; i++) {
		var m = ms[i], mNext = ms[i + 1];
		if (m*mNext <= 0) {
			c1s.push(0);
		} else {
			var dx_ = dxs[i], dxNext = dxs[i + 1], common = dx_ + dxNext;
			c1s.push(3*common/((common + dxNext)/m + (common + dx_)/mNext));
		}
	}
	c1s.push(ms[ms.length - 1]);
	
	// Get degree-2 and degree-3 coefficients
	var c2s = [], c3s = [];
	for (i = 0; i < c1s.length - 1; i++) {
		var c1 = c1s[i], m_ = ms[i], invDx = 1/dxs[i], common_ = c1 + c1s[i + 1] - m_ - m_;
		c2s.push((m_ - c1 - common_)*invDx); c3s.push(common_*invDx*invDx);
	}
	
	// Return interpolant function
	return function(x) {
		// The rightmost point in the dataset should give an exact result
		var i = xs.length - 1;
		if (x == xs[i]) { return ys[i]; }
		
		// Search for the interval x is in, returning the corresponding y if x is one of the original xs
		var low = 0, mid, high = c3s.length - 1;
		while (low <= high) {
			mid = Math.floor(0.5*(low + high));
			var xHere = xs[mid];
			if (xHere < x) { low = mid + 1; }
			else if (xHere > x) { high = mid - 1; }
			else { return ys[mid]; }
		}
		i = Math.max(0, high);
		
		// Interpolate
		var diff = x - xs[i], diffSq = diff*diff;
		return ys[i] + c1s[i]*diff + c2s[i]*diffSq + c3s[i]*diff*diffSq;
	};
};

References[edit]

External links[edit]