SciML / DataInterpolations.jl

A library of data interpolation and smoothing functions
MIT License
203 stars 43 forks source link

Wrong linear interpolation integral on interval `(t[1], t[2])` #261

Closed SouthEndMusic closed 1 month ago

SouthEndMusic commented 1 month ago

Describe the bug 🐞

I'm pretty sure the linear interpolation integral on the first interval is wrong.

Expected behavior

DataInterpolations.integral(itp, t, t) == 0 for any t, the integral is increasing when the integrand is positive.

Minimal Reproducible Example 👇

using DataInterpolations
using Random
using Plots

Random.seed!(9)
u = rand(5)
t = cumsum(rand(5))
itp = LinearInterpolation(u, t)

t_eval = range(t[1], t[end]; length = 200)
u_eval = itp.(t_eval)
u_int_eval = DataInterpolations.integral.(Ref(itp), t_eval)
u_int_numeric = cumsum((u_eval[1:(end - 1)] + u_eval[2:end]) .* t_eval.step.hi / 2)
pushfirst!(u_int_numeric, 0)

plot(t_eval, u_int_eval; label = "itp integral DataInterpolations")
plot!(t_eval, u_int_numeric; label = "itp integral numerical")
scatter!(t, DataInterpolations.integral.(Ref(itp), t); label = nothing)

plot

Suggested solution:

I don't completely understand the role of this part of the code:

https://github.com/SciML/DataInterpolations.jl/blob/5aed96dacf31344a5a73f4bbf054f571feebf96f/src/integrals.jl#L26

However, setting the first value in this tuple to -1 fixes the problem.

Environment:

  [82cc6244] DataInterpolations v5.1.0
  [64ca27bc] FindFirstFunctions v1.2.0
Julia Version 1.10.3
Commit 0b4590a550 (2024-04-30 10:59 UTC)
Build Info:
  Official https://julialang.org/ release        
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 24 × 12th Gen Intel(R) Core(TM) i7-12800HX
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, alderlake)
Threads: 16 default, 0 interactive, 8 GC (on 24 virtual cores)
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS = 16
  JULIA_PKG_PRESERVE_TIERED_INSTALLED = true
  JULIA_PKG_USE_CLI_GIT = true
sathvikbhagavan commented 1 month ago

Thanks for the report! I am able to reproduce it and looking into fixing it.