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: Demo case 1

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: Demo case 2

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: Demo case 3

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: Demo case 4

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 "Invalid input: 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"Invalid input: {name} must be a 2D list."
            for j, item in enumerate(sublist):
                if not isinstance(item, (int, float)):
                    return f"Invalid input: {name}[{i}][{j}] must be numeric."
                if math.isnan(item) or math.isinf(item):
                    return f"Invalid input: {name}[{i}][{j}] must be finite."
        return None

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

    # Validate der parameter
    if not isinstance(der, (int, float)):
        return "Invalid input: der must be an integer."
    der = int(der)
    if der < 0:
        return "Invalid input: der must be non-negative."
    if der > 2:
        return "Invalid input: 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 "Invalid input: xi and yi must have the same length."

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

    # Validate that x is not empty
    if len(x_flat) == 0:
        return "Invalid input: 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 "Invalid input: xi must be strictly increasing."

    # Call scipy.interpolate.pchip_interpolate
    try:
        result = scipy_pchip_interpolate(xi_flat, yi_flat, x_flat, der=der)
    except Exception as exc:
        return f"scipy.interpolate.pchip_interpolate error: {exc}"

    # 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 "scipy.interpolate.pchip_interpolate error: result contains non-finite values."
            return [[numeric_result]]
        return "scipy.interpolate.pchip_interpolate error: unexpected result type."

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

    return output

Overview

Excel Usage

Examples

Python Code

Online Calculator