PCHIP_INTERPOLATE

Overview

The PCHIP_INTERPOLATE function performs Piecewise Cubic Hermite Interpolating Polynomial (PCHIP) interpolation, a shape-preserving method that maintains monotonicity in the interpolated data. Unlike standard cubic spline interpolation, PCHIP guarantees that the interpolant does not overshoot or oscillate between data points, making it ideal for datasets where preserving the original data’s monotonic behavior is critical.

This implementation uses SciPy’s pchip_interpolate function from the scipy.interpolate module. The underlying algorithm is based on the work of Fritsch and Butland, which constructs local monotone piecewise cubic interpolants.

The PCHIP algorithm computes derivatives at each data point x_k using a weighted harmonic mean approach. Let h_k = x_{k+1} - x_k represent the interval width and d_k = (y_{k+1} - y_k) / h_k represent the slope. The derivative f'_k at interior points is determined as follows:

If d_k and d_{k-1} have different signs or either equals zero, then f'_k = 0
Otherwise, the derivative is computed using the weighted harmonic mean:

\frac{w_1 + w_2}{f'_k} = \frac{w_1}{d_{k-1}} + \frac{w_2}{d_k}

where w_1 = 2h_k + h_{k-1} and w_2 = h_k + 2h_{k-1}.

The resulting interpolant is C1 continuous (first derivatives are continuous), though second derivatives may have discontinuities at the knot points. This property makes PCHIP particularly useful for scientific and engineering applications where smooth, monotonic interpolation is required—such as thermodynamic property tables, dose-response curves, or any dataset where overshooting would be physically meaningless. The function also supports computing first and second derivatives at the interpolation points via the der parameter.

For more technical details, see the PchipInterpolator documentation and the original reference: F. N. Fritsch and J. Butland, “A method for constructing local monotone piecewise cubic interpolants,” SIAM J. Sci. Comput., 5(2), 300-304 (1984). DOI:10.1137/0905021.

This example function is provided as-is without any representation of accuracy.

Excel Usage

=PCHIP_INTERPOLATE(xi, yi, x, der)

xi (list[list], required): X-coordinates of the data points (must be strictly increasing)
yi (list[list], required): Y-coordinates of the data points
x (list[list], required): X-coordinates at which to evaluate the interpolated values
der (int, optional, default: 0): Order of derivative to compute (0=values, 1=first derivative, 2=second derivative)

Returns (list[list]): A 2D list (column vector) of interpolated values or derivatives, or an error message (str) if invalid.

Examples

Example 1: PCHIP interpolation at single point

Inputs:

xi	yi	x
0	0	0.5
1	1
2	4

Excel formula:

=PCHIP_INTERPOLATE({0;1;2}, {0;1;4}, 0.5)

Expected output:

Result
0.3125

Example 2: PCHIP at multiple evaluation points

Inputs:

xi	yi	x
0	1	0.5
1	2	1.5
2	0	2.5
3	3

Excel formula:

=PCHIP_INTERPOLATE({0;1;2;3}, {1;2;0;3}, {0.5;1.5;2.5})

Expected output:

Result
1.8125
1
0.8125

Example 3: First derivative computation

Inputs:

xi	yi	x	der
0	0	0	1
1	1	1
2	4	2

Excel formula:

=PCHIP_INTERPOLATE({0;1;2}, {0;1;4}, {0;1;2}, 1)

Expected output:

Result
0
1.5
4

Example 4: All parameters specified

Inputs:

xi	yi	x	der
1	2	1.5	0
2	4	2.5
3	3	3.5
4	5

Excel formula:

=PCHIP_INTERPOLATE({1;2;3;4}, {2;4;3;5}, {1.5;2.5;3.5}, 0)

Expected output:

Result
3.4375
3.5
3.5625

Python Code

import math
from scipy.interpolate import pchip_interpolate as scipy_pchip_interpolate

def pchip_interpolate(xi, yi, x, der=0):
    """
    PCHIP 1-D monotonic cubic interpolation.

    See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.pchip_interpolate.html

    This example function is provided as-is without any representation of accuracy.

    Args:
        xi (list[list]): X-coordinates of the data points (must be strictly increasing)
        yi (list[list]): Y-coordinates of the data points
        x (list[list]): X-coordinates at which to evaluate the interpolated values
        der (int, optional): Order of derivative to compute (0=values, 1=first derivative, 2=second derivative) Default is 0.

    Returns:
        list[list]: A 2D list (column vector) of interpolated values or derivatives, or an error message (str) if invalid.
    """
    def to2d(val):
        return [[val]] if not isinstance(val, list) else val

    def flatten(arr):
        result = []
        for sublist in arr:
            if not isinstance(sublist, list):
                return "Error: arguments must be 2D lists."
            for item in sublist:
                result.append(item)
        return result

    def validate_numeric_list(arr, name):
        for i, sublist in enumerate(arr):
            if not isinstance(sublist, list):
                return f"Error: {name} must be a 2D list."
            for j, item in enumerate(sublist):
                if not isinstance(item, (int, float)):
                    return f"Error: {name}[{i}][{j}] must be numeric."
                if math.isnan(item) or math.isinf(item):
                    return f"Error: {name}[{i}][{j}] must be finite."
        return None

    try:
        # Normalize inputs to 2D lists
        xi = to2d(xi)
        yi = to2d(yi)
        x = to2d(x)

        # Validate der parameter
        if not isinstance(der, (int, float)):
            return "Error: der must be an integer."
        der = int(der)
        if der < 0:
            return "Error: der must be non-negative."
        if der > 2:
            return "Error: der must be 0, 1, or 2."

        # Validate inputs are 2D lists with numeric values
        error = validate_numeric_list(xi, "xi")
        if error:
            return error
        error = validate_numeric_list(yi, "yi")
        if error:
            return error
        error = validate_numeric_list(x, "x")
        if error:
            return error

        # Flatten the 2D lists
        xi_flat = flatten(xi)
        if isinstance(xi_flat, str):
            return xi_flat
        yi_flat = flatten(yi)
        if isinstance(yi_flat, str):
            return yi_flat
        x_flat = flatten(x)
        if isinstance(x_flat, str):
            return x_flat

        # Validate that xi and yi have the same length
        if len(xi_flat) != len(yi_flat):
            return "Error: xi and yi must have the same length."

        # Validate minimum length for interpolation
        if len(xi_flat) < 2:
            return "Error: xi and yi must have at least 2 data points."

        # Validate that x is not empty
        if len(x_flat) == 0:
            return "Error: x must not be empty."

        # Validate that xi is strictly increasing (required by scipy PCHIP)
        is_increasing = all(xi_flat[i] < xi_flat[i + 1] for i in range(len(xi_flat) - 1))
        if not is_increasing:
            return "Error: xi must be strictly increasing."

        # Call scipy.interpolate.pchip_interpolate
        result = scipy_pchip_interpolate(xi_flat, yi_flat, x_flat, der=der)

        # Validate result
        if not hasattr(result, "__iter__"):
            if isinstance(result, (int, float)):
                numeric_result = float(result)
                if math.isnan(numeric_result) or math.isinf(numeric_result):
                    return "Error: result contains non-finite values."
                return [[numeric_result]]
            return "Error: unexpected result type."

        # Convert result to 2D list (column vector)
        output = []
        for val in result:
            if not isinstance(val, (int, float)):
                return "Error: result contains non-numeric values."
            numeric_val = float(val)
            if math.isnan(numeric_val) or math.isinf(numeric_val):
                return "Error: result contains non-finite values."
            output.append([numeric_val])

        return output

    except Exception as exc:
        return f"Error: {exc}"

Overview

Excel Usage

Examples

Python Code

Online Calculator