From 7740d0af45132b3576064769e397a309876c0883 Mon Sep 17 00:00:00 2001 From: Yuto Horikawa Date: Wed, 17 Feb 2021 08:07:13 +0900 Subject: [PATCH] Update documents (#96) * rename besselroots to approx_besselroots * add warning for compatibility of besselroots * use DocMeta.setdocmeta! instead of setup = * add test for besselroots * update benchmark.md to generate results automatically * fix typo in besselroots.jl * move some texts from README.md to docs/src/*.md * update ci.yml to trigger ci with release * update ci.yml for jldoctest * update ci.yml again for jldoctest * update autogeneration of benchmark * add more (a)repl block * bump version to v0.4.7 --- .github/workflows/ci.yml | 12 +++-- Project.toml | 2 +- README.md | 104 ------------------------------------ docs/Project.toml | 1 + docs/make.jl | 3 ++ docs/src/benchmark.md | 36 +++++-------- docs/src/besselroots.md | 6 +-- docs/src/gaussquadrature.md | 35 ++++++++++++ docs/src/index.md | 16 +++--- src/FastGaussQuadrature.jl | 1 + src/besselroots.jl | 18 ++++--- src/constants.jl | 4 +- src/gausschebyshev.jl | 8 +-- src/gausshermite.jl | 2 +- src/gaussjacobi.jl | 4 +- src/gausslaguerre.jl | 8 +-- src/gausslegendre.jl | 2 +- src/gausslobatto.jl | 4 +- src/gaussradau.jl | 4 +- test/test_besselroots.jl | 12 +++-- 20 files changed, 109 insertions(+), 173 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 314f5aa..23f3eed 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -1,7 +1,12 @@ name: CI on: - - push - - pull_request + pull_request: + branches: + - master + push: + branches: + - master + tags: '*' jobs: test: name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} @@ -56,8 +61,9 @@ jobs: Pkg.instantiate()' - run: | julia --project=docs -e ' - using Documenter: doctest + using Documenter using FastGaussQuadrature + DocMeta.setdocmeta!(FastGaussQuadrature, :DocTestSetup, :(using LinearAlgebra, SpecialFunctions, FastGaussQuadrature)) doctest(FastGaussQuadrature)' - run: julia --project=docs docs/make.jl env: diff --git a/Project.toml b/Project.toml index 1741923..58fbe68 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "FastGaussQuadrature" uuid = "442a2c76-b920-505d-bb47-c5924d526838" -version = "0.4.6" +version = "0.4.7" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/README.md b/README.md index 36d5d26..186f5f3 100644 --- a/README.md +++ b/README.md @@ -17,93 +17,6 @@ For a quirky account on the history of computing Gauss-Legendre quadrature, see * The fastest Julia code for Gauss quadrature nodes and weights (without tabulation). * Change the perception that Gauss quadrature rules are expensive to compute. -## Examples -Here we compute `100000` nodes and weights of the Gauss rules. -Try a million or ten million. - -```julia -julia> @time gausschebyshev( 100000 ); - 0.001788 seconds (4 allocations: 1.526 MiB) - -julia> @time gausslegendre( 100000 ); - 0.002976 seconds (10 allocations: 2.289 MiB) - -julia> @time gaussjacobi( 100000, .9, -.1 ); - 0.894373 seconds (3.59 k allocations: 1.255 GiB, 36.38% gc time) - -julia> @time gaussradau( 100000 ); - 0.684122 seconds (3.59 k allocations: 1.256 GiB, 21.71% gc time) - -julia> @time gausslobatto( 100000 ); - 0.748166 seconds (3.57 k allocations: 1.256 GiB, 27.78% gc time) - -julia> @time gausslaguerre( 100000 ); - 0.156867 seconds (7 allocations: 2.292 MiB) - -julia> @time gausshermite( 100000 ); - 0.175055 seconds (386 allocations: 67.916 MiB, 9.18% gc time) -``` - -The paper [[1]](http://epubs.siam.org/doi/abs/10.1137/140954969) computed a billion Gauss-Legendre nodes. -So here we will do a billion + 1. -```julia -julia> @time gausslegendre( 1000000001 ); - 24.441304 seconds (10 allocations: 22.352 GiB, 2.08% gc time) -``` -(The nodes near the endpoints coalesce in 16-digits of precision.) - -## The algorithm for Gauss-Chebyshev -There are four kinds of Gauss-Chebyshev quadrature rules, corresponding to four weight functions: - -1. 1st kind, weight function `w(x) = 1/sqrt(1-x^2)` - -2. 2nd kind, weight function `w(x) = sqrt(1-x^2)` - -3. 3rd kind, weight function `w(x) = sqrt((1+x)/(1-x))` - -4. 4th kind, weight function `w(x) = sqrt((1-x)/(1+x))` - -They are all have explicit simple formulas for the nodes and weights [[4]](https://books.google.co.jp/books?id=8FHf0P3to0UC). - -## The algorithm for Gauss-Legendre -Gauss quadrature for the weight function `w(x) = 1`. - -* For `n ≤ 5`: Use an analytic expression. -* For `n ≤ 60`: Use Newton's method to solve `Pₙ(x)=0`. Evaluate Legendre polynomials `Pₙ` and their derivatives `Pₙ'` by 3-term recurrence. Weights are related to `Pₙ'`. -* For `n > 60`: Use asymptotic expansions for the Legendre nodes and weights [[1]](http://epubs.siam.org/doi/abs/10.1137/140954969). - -## The algorithm for Gauss-Jacobi -Gauss quadrature for the weight functions `w(x) = (1-x)^a(1+x)^b`, `a,b > -1`. - -* For `n ≤ 100`: Use Newton's method to solve `Pₙ(x)=0`. Evaluate `Pₙ` and `Pₙ'` by three-term recurrence. -* For `n > 100`: Use Newton's method to solve `Pₙ(x)=0`. Evaluate `Pₙ` and `Pₙ'` by an asymptotic expansion (in the interior of `[-1,1]`) and the three-term recurrence `O(n^-2)` close to the endpoints. (This is a small modification to the algorithm described in [[3]](http://epubs.siam.org/doi/abs/10.1137/120889873).) -* For `max(a,b) > 5`: Use the Golub-Welsch algorithm requiring `O(n^2)` operations. - -## The algorithm for Gauss-Radau -Gauss quadrature for the weight function `w(x)=1`, except the endpoint `-1` is included as a quadrature node. - -The Gauss-Radau nodes and weights can be computed via the `(0,1)` Gauss-Jacobi nodes and weights[[3]](http://epubs.siam.org/doi/abs/10.1137/120889873). - -## The algorithm for Gauss-Lobatto -Gauss quadrature for the weight function `w(x)=1`, except the endpoints `-1` and `1` are included as nodes. - -The Gauss-Lobatto nodes and weights can be computed via the `(1,1)` Gauss-Jacobi nodes and weights[[3]](http://epubs.siam.org/doi/abs/10.1137/120889873). - -## The algorithm for Gauss-Laguerre -Gauss quadrature for the weight function `w(x) = exp(-x)` on `[0,Inf)` - -* For `n < 128`: Use the Golub-Welsch algorithm. -* For `method=GLR`: Use the Glaser-Lui-Rohklin algorithm. Evaluate Laguerre polynomials `Lₙ` and their derivatives `Lₙ'` by using Taylor series expansions near roots generated by solving the second-order differential equation that `Lₙ` satisfies, see [[2]](http://epubs.siam.org/doi/pdf/10.1137/06067016X). -* For `n ≥ 128`: Use a Newton procedure on Riemann-Hilbert asymptotics of Laguerre polynomials, see [5], based on [8]. There are some heuristics to decide which expression to use, it allows a general weight `w(x) = x^alpha exp(-q_m x^m)` and this is O(sqrt(n)) when allowed to stop when the weights are below the smallest positive floating point number. - -## The algorithm for Gauss-Hermite -Gauss quadrature for the weight function `w(x) = exp(-x^2)` on the real line. - -* For `n < 200`: Use Newton's method to solve `Hₙ(x)=0`. Evaluate Hermite polynomials `Hₙ` and their derivatives `Hₙ'` by three-term recurrence. -* For `n ≥ 200`: Use Newton's method to solve `Hₙ(x)=0`. Evaluate `Hₙ` and `Hₙ'` by a uniform asymptotic expansion, see [[7]](http://arxiv.org/abs/1410.5286). - -The paper [[7]](http://arxiv.org/abs/1410.5286) also derives an `O(n)` algorithm for generalized Gauss-Hermite nodes and weights associated to weight functions of the form `exp(-V(x))`, where `V(x)` is a real polynomial. - ## Example usage ```julia julia> @time nodes, weights = gausslegendre( 100000 ); @@ -114,20 +27,3 @@ julia> @time dot( weights, nodes.^2 ) 0.000184 seconds (7 allocations: 781.422 KiB) 0.6666666666666665 ``` - -## References: -[1] I. Bogaert, ["Iteration-free computation of Gauss-Legendre quadrature nodes and weights"](http://epubs.siam.org/doi/abs/10.1137/140954969), SIAM J. Sci. Comput., 36(3), A1008-A1026, 2014. - -[2] A. Glaser, X. Liu, and V. Rokhlin. ["A fast algorithm for the calculation of the roots of special functions."](http://epubs.siam.org/doi/pdf/10.1137/06067016X) SIAM J. Sci. Comput., 29 (2007), 1420-1438. - -[3] N. Hale and A. Townsend, ["Fast and accurate computation of Gauss-Legendre and Gauss-Jacobi quadrature nodes and weights"](http://epubs.siam.org/doi/abs/10.1137/120889873), SIAM J. Sci. Comput., 2012. - -[4] J. C. Mason and D. C. Handscomb, ["Chebyshev Polynomials"](https://books.google.co.jp/books?id=8FHf0P3to0UC), CRC Press, 2002. - -[5] P. Opsomer, (in preparation). - -[6] A. Townsend, [The race for high order Gauss-Legendre quadrature](http://pi.math.cornell.edu/~ajt/papers/QuadratureEssay.pdf), in SIAM News, March 2015. - -[7] A. Townsend, T. Trogdon, and S. Olver, ["Fast computation of Gauss quadrature nodes and weights on the whole real line"](http://arxiv.org/abs/1410.5286), to appear in IMA Numer. Anal., 2014. - -[8] M. Vanlessen, "Strong asymptotics of Laguerre-Type orthogonal polynomials and applications in Random Matrix Theory", Constr. Approx., 25:125-175, 2007. diff --git a/docs/Project.toml b/docs/Project.toml index 41f859e..a2aaa28 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,4 +1,5 @@ [deps] +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" diff --git a/docs/make.jl b/docs/make.jl index 1387079..08813ea 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,6 +1,9 @@ using Documenter using FastGaussQuadrature +# Setup for doctests in docstrings +DocMeta.setdocmeta!(FastGaussQuadrature, :DocTestSetup, :(using LinearAlgebra, SpecialFunctions, FastGaussQuadrature)) + makedocs(; modules = [FastGaussQuadrature], format = Documenter.HTML( diff --git a/docs/src/benchmark.md b/docs/src/benchmark.md index 279bed5..c533fcb 100644 --- a/docs/src/benchmark.md +++ b/docs/src/benchmark.md @@ -3,33 +3,23 @@ Here we compute `100000` nodes and weights of the Gauss rules. Try a million or ten million. -```julia -julia> @time gausschebyshev( 100000 ); - 0.001788 seconds (4 allocations: 1.526 MiB) - -julia> @time gausslegendre( 100000 ); - 0.002976 seconds (10 allocations: 2.289 MiB) - -julia> @time gaussjacobi( 100000, .9, -.1 ); - 0.894373 seconds (3.59 k allocations: 1.255 GiB, 36.38% gc time) - -julia> @time gaussradau( 100000 ); - 0.684122 seconds (3.59 k allocations: 1.256 GiB, 21.71% gc time) - -julia> @time gausslobatto( 100000 ); - 0.748166 seconds (3.57 k allocations: 1.256 GiB, 27.78% gc time) - -julia> @time gausslaguerre( 100000 ); - 0.156867 seconds (7 allocations: 2.292 MiB) - -julia> @time gausshermite( 100000 ); - 0.175055 seconds (386 allocations: 67.916 MiB, 9.18% gc time) +```@repl +using LinearAlgebra, BenchmarkTools, FastGaussQuadrature +@btime gausschebyshev(100000); +@btime gausslegendre(100000); +@btime gaussjacobi(100000, 0.9, -0.1); +@btime gaussradau(100000); +@btime gausslobatto(100000); +@btime gausslaguerre(100000); +@btime gausshermite(100000); ``` The paper [[1]](http://epubs.siam.org/doi/abs/10.1137/140954969) computed a billion Gauss-Legendre nodes. So here we will do a billion + 1. + ```julia -julia> @time gausslegendre( 1000000001 ); - 24.441304 seconds (10 allocations: 22.352 GiB, 2.08% gc time) +julia> @btime gausslegendre(1000000001); + 23.363 s (10 allocations: 22.35 GiB) ``` + (The nodes near the endpoints coalesce in 16-digits of precision.) diff --git a/docs/src/besselroots.md b/docs/src/besselroots.md index f327fc5..2852450 100644 --- a/docs/src/besselroots.md +++ b/docs/src/besselroots.md @@ -1,9 +1,9 @@ # Roots of Bessel function -Since [SpecialFunctions.jl](https://github.com/JuliaMath/SpecialFunctions.jl) doesn't have a method to calculate roots of [Bessel function](https://en.wikipedia.org/wiki/Bessel_function), we implemented `besselroots`. +Since [SpecialFunctions.jl](https://github.com/JuliaMath/SpecialFunctions.jl) doesn't have a method to calculate roots of [Bessel function](https://en.wikipedia.org/wiki/Bessel_function), we implemented `approx_besselroots`. ```@docs -besselroots(ν::Real, n::Integer) +approx_besselroots(ν::Real, n::Integer) ``` -This method `besselroots` is used to calculate `gaussjacobi` and `gausslaguerre`. +This method `approx_besselroots` is used to calculate `gaussjacobi` and `gausslaguerre`. diff --git a/docs/src/gaussquadrature.md b/docs/src/gaussquadrature.md index 1a491f4..cbadf70 100644 --- a/docs/src/gaussquadrature.md +++ b/docs/src/gaussquadrature.md @@ -2,6 +2,11 @@ ## Gauss-Legendre quadrature +Gauss quadrature for the weight function ``w(x) = 1``. + +* For ``n ≤ 5``: Use an analytic expression. +* For ``n ≤ 60``: Use Newton's method to solve ``P_n(x)=0``. Evaluate Legendre polynomials ``P_n`` and their derivatives ``P_n'`` by 3-term recurrence. Weights are related to ``P_n'``. +* For ``n > 60``: Use asymptotic expansions for the Legendre nodes and weights [[1]](http://epubs.siam.org/doi/abs/10.1137/140954969). ```@docs gausslegendre(n::Integer) @@ -9,6 +14,12 @@ gausslegendre(n::Integer) ## Gauss-Hermite quadrature +Gauss quadrature for the weight function ``w(x) = \exp(-x^2)`` on the real line. + +* For ``n < 200``: Use Newton's method to solve ``H_n(x)=0``. Evaluate Hermite polynomials ``H_n`` and their derivatives ``H_n'`` by three-term recurrence. +* For ``n ≥ 200``: Use Newton's method to solve ``H_n(x)=0``. Evaluate ``H_n`` and ``H_n'`` by a uniform asymptotic expansion, see [[7]](http://arxiv.org/abs/1410.5286). + +The paper [[7]](http://arxiv.org/abs/1410.5286) also derives an ``O(n)`` algorithm for generalized Gauss-Hermite nodes and weights associated to weight functions of the form ``\exp(-V(x))``, where ``V(x)`` is a real polynomial. ```@docs gausshermite(n::Integer) @@ -16,6 +27,11 @@ gausshermite(n::Integer) ## Gauss-Laguerre quadrature +Gauss quadrature for the weight function ``w(x) = \exp(-x)`` on ``[0,+\infty)`` + +* For ``n < 128``: Use the Golub-Welsch algorithm. +* For `method=GLR`: Use the Glaser-Lui-Rohklin algorithm. Evaluate Laguerre polynomials ``L_n`` and their derivatives ``L_n'`` by using Taylor series expansions near roots generated by solving the second-order differential equation that ``L_n`` satisfies, see [[2]](http://epubs.siam.org/doi/pdf/10.1137/06067016X). +* For ``n ≥ 128``: Use a Newton procedure on Riemann-Hilbert asymptotics of Laguerre polynomials, see [5], based on [8]. There are some heuristics to decide which expression to use, it allows a general weight ``w(x) = x^\alpha \exp(-q_m x^m)`` and this is ``O(\sqrt{n})`` when allowed to stop when the weights are below the smallest positive floating point number. ```@docs gausslaguerre(n::Integer) @@ -27,6 +43,14 @@ gausslaguerre(n::Integer, α::Real) ## Gauss-Chebyshev quadrature +There are four kinds of Gauss-Chebyshev quadrature rules, corresponding to four weight functions: + +* 1st kind, weight function ``w(x) = 1/\sqrt{1-x^2}`` +* 2nd kind, weight function ``w(x) = \sqrt{1-x^2}`` +* 3rd kind, weight function ``w(x) = \sqrt{(1+x)/(1-x)}`` +* 4th kind, weight function ``w(x) = \sqrt{(1-x)/(1+x)}`` + +They are all have explicit simple formulas for the nodes and weights [[4]](https://books.google.co.jp/books?id=8FHf0P3to0UC). ```@docs gausschebyshev(n::Integer, kind::Integer) @@ -34,6 +58,11 @@ gausschebyshev(n::Integer, kind::Integer) ## Gauss-Jacobi quadrature +Gauss quadrature for the weight functions ``w(x) = (1-x)^\alpha(1+x)^\beta``, ``\alpha,\beta > -1``. + +* For ``n ≤ 100``: Use Newton's method to solve ``P_n(x)=0``. Evaluate ``P_n`` and ``P_n'`` by three-term recurrence. +* For ``n > 100``: Use Newton's method to solve ``P_n(x)=0``. Evaluate ``P_n`` and ``P_n'`` by an asymptotic expansion (in the interior of ``[-1,1]``) and the three-term recurrence ``O(n^{-2})`` close to the endpoints. (This is a small modification to the algorithm described in [[3]](http://epubs.siam.org/doi/abs/10.1137/120889873).) +* For ``max(a,b) > 5``: Use the Golub-Welsch algorithm requiring ``O(n^2)`` operations. ```@docs gaussjacobi(n::Integer, α::Real, β::Real) @@ -41,6 +70,9 @@ gaussjacobi(n::Integer, α::Real, β::Real) ## Gauss-Radau quadrature +Gauss quadrature for the weight function ``w(x)=1``, except the endpoint ``-1`` is included as a quadrature node. + +The Gauss-Radau nodes and weights can be computed via the ``(0,1)`` Gauss-Jacobi nodes and weights[[3]](http://epubs.siam.org/doi/abs/10.1137/120889873). ```@docs gaussradau(n::Integer) @@ -48,6 +80,9 @@ gaussradau(n::Integer) ## Gauss-Lobatto quadrature +Gauss quadrature for the weight function ``w(x)=1``, except the endpoints ``-1`` and ``1`` are included as nodes. + +The Gauss-Lobatto nodes and weights can be computed via the ``(1,1)`` Gauss-Jacobi nodes and weights[[3]](http://epubs.siam.org/doi/abs/10.1137/120889873). ```@docs gausslobatto(n::Integer) diff --git a/docs/src/index.md b/docs/src/index.md index 8594530..d4cfdd2 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -20,15 +20,11 @@ To check an integral \int_{-1}^{1} x^4 dx = \frac{2}{5} ``` by numerically, try following code. -```julia -julia> using FastGaussQuadrature, LinearAlgebra -julia> x, w = gausslegendre(3); - -julia> f(x) = x^4; - -julia> I = dot(w, f.(x)); - -julia> I ≈ 2/5 -true +```@repl +using FastGaussQuadrature, LinearAlgebra +x, w = gausslegendre(3) +f(x) = x^4 +I = dot(w, f.(x)) +I ≈ 2/5 ``` diff --git a/src/FastGaussQuadrature.jl b/src/FastGaussQuadrature.jl index f37094f..ffec28e 100644 --- a/src/FastGaussQuadrature.jl +++ b/src/FastGaussQuadrature.jl @@ -11,6 +11,7 @@ export gausshermite export gaussjacobi export gausslobatto export gaussradau +export approx_besselroots export besselroots import SpecialFunctions: besselj, airyai, airyaiprime diff --git a/src/besselroots.jl b/src/besselroots.jl index c4f68c5..c33eba9 100644 --- a/src/besselroots.jl +++ b/src/besselroots.jl @@ -1,17 +1,18 @@ @doc raw""" - gausslegendre(ν::Real, n::Integer) -> Vector{Float64} + approx_besselroots(ν::Real, n::Integer) -> Vector{Float64} Return the first ``n`` roots of [Bessel function](https://en.wikipedia.org/wiki/Bessel_function). +Note that this function is only 12-digits of precision. ```math J_{\nu}(x) = \sum_{m=0}^{\infty}\frac{(-1)^j}{\Gamma(\nu+j+1)j!} \left(\frac{x}{2}\right)^{2j+\nu} ``` # Examples -```jldoctest; setup = :(using FastGaussQuadrature, SpecialFunctions) +```jldoctest julia> ν = 0.3; -julia> roots = besselroots(ν, 10); +julia> roots = approx_besselroots(ν, 10); julia> zeros = (x -> besselj(ν, x)).(roots); @@ -19,12 +20,12 @@ julia> all(zeros .< 1e-12) true ``` """ -function besselroots(ν::Real, n::Integer) +function approx_besselroots(ν::Real, n::Integer) # FIXME (related issue #22 and #80) - return besselroots(Float64(ν), n) + return approx_besselroots(Float64(ν), n) end -function besselroots(ν::Float64, n::Integer) +function approx_besselroots(ν::Float64, n::Integer) # DEVELOPERS NOTES: # ν = 0 --> Full Float64 precision for n ≤ 20 (Wolfram Alpha), and very # accurate approximations for n > 20 (McMahon's expansion) @@ -219,3 +220,8 @@ function besselJ1(m) end return Jk2 end + +function besselroots(ν::Real, n::Integer) + @warn "`besselroots` has been renamed to `approx_besselroots`, and `besselroots` will be removed in the next major release." + return approx_besselroots(ν, n) +end diff --git a/src/constants.jl b/src/constants.jl index bf1d014..c882c9e 100644 --- a/src/constants.jl +++ b/src/constants.jl @@ -71,8 +71,8 @@ const Piessens_C = @SMatrix [ Values of Bessel function ``J_1`` on first ten roots of Bessel function `J_0`. # Examples -```jldoctest; setup = :(using FastGaussQuadrature, SpecialFunctions) -julia> roots = besselroots(0,10); +```jldoctest +julia> roots = approx_besselroots(0,10); julia> (besselj1.(roots)).^2 ≈ FastGaussQuadrature.besselJ1_10 true diff --git a/src/gausschebyshev.jl b/src/gausschebyshev.jl index 76d5128..f860e3a 100644 --- a/src/gausschebyshev.jl +++ b/src/gausschebyshev.jl @@ -9,7 +9,7 @@ Return nodes and weights of [Gauss-Chebyshev quadrature](https://en.wikipedia.or ``` # Examples -```jldoctest; setup = :(using FastGaussQuadrature, LinearAlgebra) +```jldoctest julia> x, w = gausschebyshev(3); julia> f(x) = x^4; @@ -31,7 +31,7 @@ Return nodes and weights of [Gauss-Chebyshev quadrature](https://en.wikipedia.or ``` # Examples -```jldoctest; setup = :(using FastGaussQuadrature, LinearAlgebra) +```jldoctest julia> x, w = gausschebyshev(3, 2); julia> f(x) = x^4; @@ -53,7 +53,7 @@ Return nodes and weights of [Gauss-Chebyshev quadrature](https://en.wikipedia.or ``` # Examples -```jldoctest; setup = :(using FastGaussQuadrature, LinearAlgebra) +```jldoctest julia> x, w = gausschebyshev(3, 3); julia> f(x) = x^4; @@ -75,7 +75,7 @@ Return nodes and weights of [Gauss-Chebyshev quadrature](https://en.wikipedia.or ``` # Examples -```jldoctest; setup = :(using FastGaussQuadrature, LinearAlgebra) +```jldoctest julia> x, w = gausschebyshev(3, 4); julia> f(x) = x^4; diff --git a/src/gausshermite.jl b/src/gausshermite.jl index 8f60b5d..92fcd7c 100644 --- a/src/gausshermite.jl +++ b/src/gausshermite.jl @@ -8,7 +8,7 @@ Return nodes and weights of [Gauss-Hermite quadrature](https://en.wikipedia.org/ ``` # Examples -```jldoctest; setup = :(using FastGaussQuadrature, LinearAlgebra) +```jldoctest julia> x, w = gausshermite(3); julia> f(x) = x^4; diff --git a/src/gaussjacobi.jl b/src/gaussjacobi.jl index 3fb6e85..0aace56 100644 --- a/src/gaussjacobi.jl +++ b/src/gaussjacobi.jl @@ -8,7 +8,7 @@ Return nodes and weights of [Gauss-Jacobi quadrature](https://en.wikipedia.org/w ``` # Examples -```jldoctest; setup = :(using FastGaussQuadrature, LinearAlgebra) +```jldoctest julia> x, w = gaussjacobi(3, 1/3, -1/3); julia> f(x) = x^4; @@ -374,7 +374,7 @@ function boundary(n::Integer, α::Float64, β::Float64, npts::Integer) # Use Newton iterations to find the first few Bessel roots: smallK = min(30, npts) - jk = besselroots(α, smallK) + jk = approx_besselroots(α, smallK) # Approximate roots via asymptotic formula: (see Olver 1974) phik = jk/(n + .5*(α + β + 1)) diff --git a/src/gausslaguerre.jl b/src/gausslaguerre.jl index 09e32d1..241b48e 100644 --- a/src/gausslaguerre.jl +++ b/src/gausslaguerre.jl @@ -8,7 +8,7 @@ Return nodes and weights of [Gauss-Laguerre quadrature](https://en.wikipedia.org ``` # Examples -```jldoctest; setup = :(using FastGaussQuadrature, LinearAlgebra) +```jldoctest julia> x, w = gausslaguerre(3); julia> f(x) = x^4; @@ -34,7 +34,7 @@ Return nodes and weights of generalized [Gauss-Laguerre quadrature](https://en.w ``` # Examples -```jldoctest; setup = :(using FastGaussQuadrature, LinearAlgebra) +```jldoctest julia> x, w = gausslaguerre(3, 1.0); julia> f(x) = x^4; @@ -130,7 +130,7 @@ function gausslaguerre_asy(n::Integer, α; # The Bessel region # First compute the roots of the Bessel function of order α - jak_vector = besselroots(α, k_bessel) + jak_vector = approx_besselroots(α, k_bessel) bessel_wins = true k = 0 @@ -642,7 +642,7 @@ function gausslaguerre_rec(n, α; reduced = false) n_pre = min(n, 7) ν = 4n + 2α + 2 - x_pre = T.(besselroots(α, n_pre)).^2 / ν # this is a lower bound by [DLMF 18.16.10] + x_pre = T.(approx_besselroots(α, n_pre)).^2 / ν # this is a lower bound by [DLMF 18.16.10] noUnderflow = true # this flag turns false once the weights start to underflow for k in 1:n diff --git a/src/gausslegendre.jl b/src/gausslegendre.jl index 0fd1d05..103fd77 100644 --- a/src/gausslegendre.jl +++ b/src/gausslegendre.jl @@ -8,7 +8,7 @@ Return nodes and weights of [Gauss-Legendre quadrature](https://en.wikipedia.org ``` # Examples -```jldoctest; setup = :(using FastGaussQuadrature, LinearAlgebra) +```jldoctest julia> x, w = gausslegendre(3); julia> f(x) = x^4; diff --git a/src/gausslobatto.jl b/src/gausslobatto.jl index d454474..81cd525 100644 --- a/src/gausslobatto.jl +++ b/src/gausslobatto.jl @@ -8,7 +8,7 @@ Return nodes and weights of [Gauss-Lobatto quadrature](https://mathworld.wolfram ``` # Examples -```jldoctest; setup = :(using FastGaussQuadrature, LinearAlgebra) +```jldoctest julia> x, w = gausslobatto(4); julia> f(x) = x^4; @@ -21,7 +21,7 @@ true Note that the both ends of nodes are fixed at -1 and 1. -```jldoctest; setup = :(using FastGaussQuadrature, LinearAlgebra) +```jldoctest julia> x, w = gausslobatto(4); julia> x[1], x[end] diff --git a/src/gaussradau.jl b/src/gaussradau.jl index 42ab911..50e0bc3 100644 --- a/src/gaussradau.jl +++ b/src/gaussradau.jl @@ -8,7 +8,7 @@ Return nodes and weights of [Gauss-Radau quadrature](https://mathworld.wolfram.c ``` # Examples -```jldoctest; setup = :(using FastGaussQuadrature, LinearAlgebra) +```jldoctest julia> x, w = gaussradau(3); julia> f(x) = x^4; @@ -21,7 +21,7 @@ true Note that the first node is fixed at -1. -```jldoctest; setup = :(using FastGaussQuadrature, LinearAlgebra) +```jldoctest julia> x, w = gaussradau(3); julia> x[1] diff --git a/test/test_besselroots.jl b/test/test_besselroots.jl index 1251c96..f046371 100644 --- a/test/test_besselroots.jl +++ b/test/test_besselroots.jl @@ -2,14 +2,16 @@ import SpecialFunctions import SpecialFunctions: besselj @testset "Bessel Roots" begin - @test_throws DomainError besselroots(0.0, -1) - - tol = 1e-11 + @test_throws DomainError approx_besselroots(0.0, -1) + + msg = "`besselroots` has been renamed to `approx_besselroots`, and `besselroots` will be removed in the next major release." + @test_logs (:warn, msg) besselroots(0.2,3) - # Check if besselj(ν, besselroots(ν, n) ) is small: + # Check if besselj(ν, approx_besselroots(ν, n) ) is small: + tol = 1e-11 for ν = 0.:0.1:5. n = 10 - @test norm( besselj.(ν, besselroots(ν, n) ), Inf ) < tol + @test norm( besselj.(ν, approx_besselroots(ν, n) ), Inf ) < tol end end