Skip to content

Smooth equivalent of ConstantInterpolation in terms of integral #364

Open
@SouthEndMusic

Description

@SouthEndMusic

Is your feature request related to a problem? Please describe.

A bit of context: In my ODE problem for hydrological modelling, I want to model certain forcing events like precipitation. The precipitation is given as a time series, where constant interpolation ('forward fill') is assumed so that the exact total volume is known. Constant interpolation however does not play nice with ODE solvers because of its lack of smoothness. Therefore I want to have an interpolation type which:

  • Is at least $C^1$ smooth;
  • Has for each interval between data points the same integral as forward fill would have.

An additional requirement I have is that it is easy to extend this interpolation with new incoming data.

Describe the solution you’d like

The naive solution of assuming constant interpolation on the first interval and then quadratic on the others to satisfy the above conditions is not stable:

using DataInterpolations
using CairoMakie
using Random
using FindFirstFunctions: searchsortedfirst
using QuadGK

fig = Figure()

ax_itp = Axis(fig[1,1]; title = "Interpolation")
ax_int = Axis(fig[2,1]; title = "Integral")

N = 6
Random.seed!(5)
u = rand(N)
t = cumsum(rand(N))
t_eval = range(first(t), last(t), length = 250)

A = ConstantInterpolation(u, t; cache_parameters = true)
lines!(ax_itp, t_eval, A.(t_eval); label = "ConstantInterpolation")
lines!(ax_int, t_eval, DataInterpolations.integral.(Ref(A), t_eval); label = "ConstantInterpolation")

a = zeros(N - 1)
b = zeros(N - 1)
c = zeros(N - 1)

# Assume constant interpolation on first interval
c[1] = u[1]

for i in 1:(N-2)
    Δt = t[i+1] - t[i]

    # value at t[i+1]
    u_ = @evalpoly Δt c[i] b[i] a[i]

    # derivative at t[i+1]
    du_ = @evalpoly Δt b[i] 2a[i]

    # integral over (t[i+1], t[i+2])
    Δt = t[i+2] - t[i+1]  
    I = Δt * u[i+1]

    # Next coefficients
    c[i+1] = u_
    b[i+1] = du_
    a[i+1] = 3 * (u[i+1] - u_ - (du_ * Δt) / 2) / Δt^2
end

function B(t_)
    i = clamp(searchsortedfirst(t, t_)-1, 1, N-1)
    @evalpoly (t_ - t[i]) c[i] b[i] a[i]
end

lines!(ax_itp, t_eval, B.(t_eval); label = "Naive smooth interpolation")
lines!(ax_int, t_eval, [quadgk(B, first(t), t_)[1] for t_ in t_eval]; label = "Naive smooth interpolation")

axislegend(ax_itp)
axislegend(ax_int)

fig

figure

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions