From 73ea895598a37c915baff1d66680f6dc298ca0f1 Mon Sep 17 00:00:00 2001 From: Etienne Date: Mon, 31 Mar 2025 17:12:58 +0200 Subject: [PATCH] fix pchip endpoints --- src/interpolation_utils.jl | 28 +++++++++++++++++++++++----- 1 file changed, 23 insertions(+), 5 deletions(-) diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index bf06a19e..2caa7f48 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -70,7 +70,7 @@ function quadratic_spline_params(t::AbstractVector, sc::AbstractVector) # Create linear system Ac = u, where: # - A consists of basis function evaulations in t - # - c are 1D control points + # - c are 1D control points n = length(t) dtype_sc = typeof(t[1] / t[1]) @@ -235,6 +235,19 @@ function du_PCHIP(u, t) δ = diff(u) ./ h s = sign.(δ) + # Special handling of the slope at the endpoints, see + # Cleve Moler, Numerical Computing with MATLAB, Chap 3.6 (file pchiptx.m, function pchipend()) + function _edge_case(h₁, h₂, δ₁, δ₂) + d = ((2 * h₁ + h₂) * δ₁ - h₁ * δ₂) / (h₁ + h₂) + if sign(d) != sign(δ₁) + zero(eltype(δ)) + elseif sign(δ₁) != sign(δ₂) && abs(d) > 3 * abs(δ₁) + 3 * δ₁ + else + d + end + end + function _du(k) sₖ₋₁, sₖ = if k == 1 s[1], s[2] @@ -248,17 +261,22 @@ function du_PCHIP(u, t) zero(eltype(δ)) elseif sₖ₋₁ == sₖ if k == 1 - ((2 * h[1] + h[2]) * δ[1] - h[1] * δ[2]) / (h[1] + h[2]) + _edge_case(h[1], h[2], δ[1], δ[2]) elseif k == lastindex(t) - ((2 * h[end] + h[end - 1]) * δ[end] - h[end] * δ[end - 1]) / - (h[end] + h[end - 1]) + _edge_case(h[end], h[end - 1], δ[end], δ[end - 1]) else w₁ = 2h[k] + h[k - 1] w₂ = h[k] + 2h[k - 1] (w₁ + w₂) / (w₁ / δ[k - 1] + w₂ / δ[k]) end else - zero(eltype(δ)) + if k == 1 + _edge_case(h[1], h[2], δ[1], δ[2]) + elseif k == lastindex(t) + _edge_case(h[end], h[end - 1], δ[end], δ[end - 1]) + else + zero(eltype(δ)) + end end end