kbarbary / Dierckx.jl

Julia package for 1-d and 2-d splines
Other
157 stars 30 forks source link

`roots` sometimes cannot find the correct root(s) #89

Open islent opened 2 years ago

islent commented 2 years ago

MWE:

using Dierckx
x = [
    -7.800000000000001,
    -7.7,
    -7.6000000000000005,
    -7.5,
    -7.4,
    -7.300000000000001,
];

y = [
    -0.02749549159389436,
    -0.01838389014507258,
    -0.00921755334192842,
     1.1102230246251565e-16,
     0.009265146320801443,
     0.018574158954026365,
];

spline = Spline1D(x,y)
roots(spline) #! No root

spline = Spline1D(x,y,s=1.0)
roots(spline) # correct

y2 = [
    -0.02749549159389436,
    -0.01838389014507258,
    -0.00921755334192842,
     0.0,
     0.009265146320801443,
     0.018574158954026365,
];

spline = Spline1D(x,y2)
roots(spline) # correct

y3 = [
    -0.02749549159389436,
    -0.01838389014507258,
    -0.00921755334192842,
    1.1102230246251565e-14,
     0.009265146320801443,
     0.018574158954026365,
];

spline = Spline1D(x,y3)
roots(spline) #! No root

y4 = [
    -0.02749549159389436,
    -0.01838389014507258,
    -0.00921755334192842,
    1.1102230246251565e-13,
     0.009265146320801443,
     0.018574158954026365,
];

spline = Spline1D(x,y4)
roots(spline) # correct

Output:

julia> using Dierckx

julia> x = [                                                                                                                                              
           -7.800000000000001,                                                                                                                            
           -7.7,                                                                                                                                          
           -7.6000000000000005,                                                                                                                           
           -7.5,                                                                                                                                          
           -7.4,                                                                                                                                          
           -7.300000000000001,                                                                                                                            
       ];

julia> 

julia> y = [                                                                                                                                              
           -0.02749549159389436,                                                                                                                          
               -0.01838389014507258,                                                                                                                      
           -0.00921755334192842,                                                                                                                          
            1.1102230246251565e-16,                                                                                                                       
            0.009265146320801443,                                                                                                                         
            0.018574158954026365,                                                                                                                         
       ];

julia> 

julia> spline = Spline1D(x,y)
Spline1D(knots=[-7.800000000000001, -7.6000000000000005, -7.5, -7.300000000000001], k=3, extrapolation="nearest", residual=0.0)

julia> roots(spline) #! No root
Float64[]

julia> 

julia> spline = Spline1D(x,y,s=1.0)
Spline1D(knots=[-7.800000000000001, -7.300000000000001], k=3, extrapolation="nearest", residual=1.5429392612813467e-15)

julia> y4 = [
           -0.02749549159389436,
           -0.01838389014507258,
           -0.00921755334192842,
           1.1102230246251565e-13,
            0.009265146320801443,
            0.018574158954026365,
       ];

julia>

julia> spline = Spline1D(x,y4)
Spline1D(knots=[-7.800000000000001, -7.6000000000000005, -7.5, -7.300000000000001], k=3, extrapolation="nearest", residual=0.0)

julia> roots(spline) # correct
1-element Vector{Float64}:
 -7.500000000000101

Environment:

julia> versioninfo()
Julia Version 1.7.2
Commit bf53498635 (2022-02-06 15:21 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Xeon(R) W-10885M CPU @ 2.40GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, skylake)
Environment:
  JULIA_DEPOT_PATH = E:\.julia
  JULIA_NUM_THREADS = 8
  JULIA_EDITOR = code

(@v1.7) pkg> st Dierckx
      Status `E:\.julia\environments\v1.7\Project.toml`
  [39dd38d3] Dierckx v0.5.2