diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index 00a62a02..61b4d51d 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -134,22 +134,30 @@ end function du_PCHIP(u, t) h = diff(u) δ = h ./ diff(t) + s = sign.(δ) function _du(k) - if k == 1 - return ((2 * h[1] + h[2]) * δ[1] - h[1] * δ[2]) / (h[1] + h[2]) - elseif k == length(t) - return ((2 * h[end] + h[end - 1]) * δ[end] - h[end] * δ[end - 1]) / - (h[end] + h[end - 1]) + sₖ₋₁, sₖ = if k == 1 + s[1], s[2] + elseif k == lastindex(t) + s[end - 1], s[end] + else + s[k - 1], s[k] end - sₖ₋₁ = sign(δ[k - 1]) - sₖ = sign(δ[k]) + if sₖ₋₁ == 0 && sₖ == 0 zero(eltype(δ)) elseif sₖ₋₁ == sₖ - w₁ = 2h[k] + h[k - 1] - w₂ = h[k] + 2h[k - 1] - δ[k - 1] * δ[k] * (w₁ + w₂) / (w₁ * δ[k] + w₂ * δ[k - 1]) + if k == 1 + ((2 * h[1] + h[2]) * δ[1] - h[1] * δ[2]) / (h[1] + h[2]) + elseif k == lastindex(t) + ((2 * h[end] + h[end - 1]) * δ[end] - h[end] * δ[end - 1]) / + (h[end] + h[end - 1]) + else + w₁ = 2h[k] + h[k - 1] + w₂ = h[k] + 2h[k - 1] + δ[k - 1] * δ[k] * (w₁ + w₂) / (w₁ * δ[k] + w₂ * δ[k - 1]) + end else zero(eltype(δ)) end