gwater / struve.jl Goto Github PK

3.0 3.0 4.0 168 KB

Struve functions for Julia

License: MIT License

Julia 100.00%

struve.jl's Introduction

Struve.jl

This package provides methods to compute the Struve functions H, K, L, and M.

The default methods currently use a mixture of power series, large argument expansions, and integral representations of the functions which are computed numerically using QuadGK.jl and SpecialFunctions.jl.

Struve.jl exports four functions:

struveh(ν, z)
struvek(ν, z)
sruvel(ν, z)
struvem(ν, z)

which compute the Struve functions of the first and second kind (struveh and struvek), and the modified Struve functions of the first and second kind (struvel and struvem).

It also implements fast approximations for H₀ and H₁ on the real axis (with absolute error below 2×10⁻³). For fast, high accuracy approximations ApproxFun.jl may be used.

Struve.H0_fast(x)
Struve.H1_fast(x)

Please note: Implementations are verified against function's power series computed in higher precision as well as explicit forms. In general, we try to provide maximum relative errors throughout the entire range better than 1e-13 in double precision. There may be some regions around cutoffs that are less accurate. Bug reports are very welcome as it is difficult to test over all ranges of order and argument. Rigorous implementations of Struve functions are difficult to find but comparisons to Mathematica or Keisan would be helpful. The power series can also be called directly (e.g., Struve.struveh_power_series(big"1.5", big"90.0")) using higher precision which will be a good comparison for the lower precisions calculations.

struve.jl's People

Contributors

Stargazers

Watchers

Forkers

yoshinobu-ishizaki juliatagbot standardgalactic heltonmc

struve.jl's Issues

Struve and related functions?

I'm wondering what your thoughts are on including the 'related functions' to Struve functions here?

On NIST chapter 11 (https://dlmf.nist.gov/11) it is Struve and related functions. These might include

Struve functions
Lommel
Anger
Weber
Anger-Weber
Kelvin (compute with imaginary bessel functions)
Scorer (might be best with Airy functions)

I do not have implementations ready for these yet to contribute, but would like to define the scope of this package. I would really like to build up the special function libraries in Julia and would think this would be a nice companion to Bessels.jl that I can link to from there!

verify integrals

implement test cases based on https://dlmf.nist.gov/11.7.ii

fix 0.7 deprecations

verify zeros

using tables in https://doi.org/10.1016/S0377-0427(01)00580-5

Bug for x=0

When trying

Struve.H(0,z) I get the error

integrand produced NaN in the interval (0.9999999999999929, 1.0)
Stacktrace:
  [1] evalrule(f::QuadGK.var"#14#23"{Struve.var"#3#4"{Float64, Int64, Float64}, Float64, Float64}, a::Float64, b::Float64, x::Vector{Float64}, w::Vector{Float64}, gw::Vector{Float64}, nrm::typeof(norm))
    @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/evalrule.jl:37
  [2] adapt(f::QuadGK.var"#14#23"{Struve.var"#3#4"{Float64, Int64, Float64}, Float64, Float64}, segs::Vector{QuadGK.Segment{Float64, Float64, Float64}}, I::Float64, E::Float64, numevals::Int64, x::Vector{Float64}, w::Vector{Float64}, gw::Vector{Float64}, n::Int64, atol::Float64, rtol::Float64, maxevals::Int64, nrm::typeof(norm))
    @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:48
  [3] do_quadgk(f::QuadGK.var"#14#23"{Struve.var"#3#4"{Float64, Int64, Float64}, Float64, Float64}, s::Tuple{Float64, Float64}, n::Int64, atol::Nothing, rtol::Nothing, maxevals::Int64, nrm::typeof(norm), segbuf::Nothing)
    @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:37
  [4] #46
    @ ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:219 [inlined]
  [5] handle_infinities(workfunc::QuadGK.var"#46#47"{Nothing, Nothing, Int64, Int64, typeof(norm), Nothing}, f::Struve.var"#3#4"{Float64, Int64, Float64}, s::Tuple{Float64, Float64})
    @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:110
  [6] #quadgk#45
    @ ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:218 [inlined]
  [7] quadgk
    @ ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:216 [inlined]
  [8] integrate
    @ ~/.julia/packages/Struve/BgIsA/src/integral.jl:5 [inlined]
  [9] _K_integral(ν::Int64, z::Float64)
    @ Struve ~/.julia/packages/Struve/BgIsA/src/integral.jl:13
 [10] _H_integral(ν::Int64, z::Float64)
    @ Struve ~/.julia/packages/Struve/BgIsA/src/integral.jl:17
 [11] Laplace2D(ξ::Vector{Float64}, R::Int64)
    @ Main ~/Documents/Julia_develop/InterpFourierChebyshev/src/KernelsFourier.jl:15
 [12] top-level scope
    @ ~/Documents/Julia_develop/InterpFourierChebyshev/src/KernelsFourier.jl:30

verify inequilities

implement test cases, like https://dlmf.nist.gov/11.4.E13

Thoughts on updating this package?

Hi there @gwater !! Thanks for contributing this package. I came across this while thinking about bringing together native Julia implementations of special functions and this caught my interest due to the overlap with other functions.

I was wondering how you felt about updating this package or if you are still interested in developing some of these functions? I read through the code base and noticed there are a couple of places where it might be possible to make some optimizations. In general, I would be open to contributing some PRs for this. Would that be ok?

I don't want to change too much if you are happy with this so let me know.

Some of the changes I have in mind would be to - (1) update and optimize all the power series forms so that they can be run to a desired tolerance (and be much faster), (2) not rely on QuadGK or any integration tools that require allocation, (3) use some of the relations to bessel function (Bessels.jl) to give some of the identities.

Let me know your general thoughts on the future of this package and where you would want to see it go!!
-Michael

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!