Description
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