JuliaApproximation / ApproxFun.jl

Julia package for function approximation
http://juliaapproximation.github.io/ApproxFun.jl/
Other
542 stars 70 forks source link

abs of higher order derivative fails. #643

Open macd opened 5 years ago

macd commented 5 years ago

The abs() fails on a forth order derivative where is succeeds on the third order. This looks like a regression since it did work in the November timeframe (this is on ApproxFun master). Details below.

   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.2.0-DEV.163 (2019-01-13)
 _/ |\__'_|_|_|\__'_|  |  Commit b56a9f0794 (0 days old master)
|__/                   |

julia (v1.2)> using ApproxFun

julia (v1.2)> x = Fun()
Fun(Chebyshev(),[0.0, 1.0])

julia (v1.2)> f = abs(sin(5x))^3
Fun(ChebyshevDirichlet{1,1,Segment{Float64},Float64}(the segment [-1.0,-0.6283185307179585])⨄ChebyshevDirichlet{1,1,Segment{Float64},Float64}(the segment [-0.6283185307179585,-5.070867905569396e-16])⨄ChebyshevDirichlet{1,1,Segment{Float64},Float64}(the segment [-5.070867905569396e-16,0.6283185307179586])⨄ChebyshevDirichlet{1,1,Segment{Float64},Float64}(the segment [0.6283185307179586,1.0]),[0.440883, 8.36157e-45, 8.90973e-45, 0.440883, -0.440883, 9.54702e-45, -8.91397e-45, 0.440883, -0.0527034, -0.287537  …  -2.16105e-42, 0.0, 8.99137e-50, 9.92937e-47, -5.65047e-48, 0.0, -1.29898e-49, 1.44048e-45, 1.89168e-45, 0.0])

julia (v1.2)> fd3 = f'''
Fun(Ultraspherical(3,the segment [-1.0,-0.6283185307179585])⨄Ultraspherical(3,the segment [-0.6283185307179585,-5.070867905569396e-16])⨄Ultraspherical(3,the segment [-5.070867905569396e-16,0.6283185307179586])⨄Ultraspherical(3,the segment [0.6283185307179586,1.0]),[527.348, -7.99886e-13, 2.97565e-13, -527.348, -79.7517, 212.82, 212.82, -79.7517, -92.9299, 5.55978e-13  …  -2.56711e-38, 0.0, 5.26734e-45, 1.20409e-42, -6.85209e-44, 0.0, -7.7716e-45, 1.78398e-41, 2.34277e-41, 0.0])

julia (v1.2)> fd4 = f''''
Fun(Ultraspherical(4,the segment [-1.0,-0.6283185307179585])⨄Ultraspherical(4,the segment [-0.6283185307179585,-5.070867905569396e-16])⨄Ultraspherical(4,the segment [-5.070867905569396e-16,0.6283185307179586])⨄Ultraspherical(4,the segment [0.6283185307179586,1.0]),[-2574.84, 4064.57, 4064.57, -2574.84, -3000.31, 1.06184e-11, -2.03649e-12, 3000.31, 317.319, -1533.12  …  -4.90281e-37, 0.0, 1.7006e-43, 2.29965e-41, -1.30865e-42, 0.0, -2.50912e-43, 3.40715e-40, 4.47437e-40, 0.0])

julia (v1.2)> abs(fd3)
Fun(ApproxFun.ContinuousSpace{Float64,Float64}(ApproxFun.PiecewiseSegment{Float64}([-1.0, -0.942478, -0.726495, -0.628319]))⨄ApproxFun.ContinuousSpace{Float64,Float64}(ApproxFun.PiecewiseSegment{Float64}([-0.628319, -0.530142, -0.314159, -0.0981765, -5.07087e-16]))⨄ApproxFun.ContinuousSpace{Float64,Float64}(ApproxFun.PiecewiseSegment{Float64}([-5.07087e-16, 0.0981765, 0.314159, 0.530142, 0.628319]))⨄ApproxFun.ContinuousSpace{Float64,Float64}(ApproxFun.PiecewiseSegment{Float64}([0.628319, 0.726495, 0.942478, 1.0])),[667.58, 1.36424e-12, 4.583e-13, 750.0, -1.13687e-12, 5.74651e-13, 1.00986e-12, 750.0, 1.36424e-12, -1.10134e-13  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.60548e-13])

julia (v1.2)> abs(fd4)
ERROR: ArgumentError: Domain cannot be empty
Stacktrace:
 [1] Type at /home/macd/.julia/dev/ApproxFun/src/Spaces/Chebyshev/Chebyshev.jl:16 [inlined]
 [2] Type at /home/macd/.julia/dev/ApproxFun/src/Spaces/Chebyshev/Chebyshev.jl:22 [inlined]
 [3] canonicalspace at /home/macd/.julia/dev/ApproxFun/src/Spaces/Ultraspherical/DirichletSpace.jl:30 [inlined]
 [4] plan_transform(::ChebyshevDirichlet{1,1,Segment{Float64},Float64}, ::Array{Float64,1}) at /home/macd/.julia/dev/ApproxFun/src/Fun/Space.jl:396
 [5] transform(::ChebyshevDirichlet{1,1,Segment{Float64},Float64}, ::Array{Float64,1}) at /home/macd/.julia/dev/ApproxFun/src/Fun/Space.jl:426
 [6] *(::ApproxFun.TransformPlan{Float64,ApproxFun.ContinuousSpace{Float64,Float64},false,Nothing}, ::Array{Float64,1}) at /home/macd/.julia/dev/ApproxFun/src/Spaces/Ultraspherical/ContinuousSpace.jl:79
...
dlfivefifty commented 5 years ago

High order zeros are a known problem with the approach to root finding, and so the code is quite brittle (hence why it used to work but doesn't anymore). It's not clear to me how to overcome this, though this particular error looks like it's not filtering out empty intervals correctly.