From f849563c0befc20b8d6b7016fa97d55d78d95ce7 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Thu, 2 Feb 2017 22:06:05 +1100 Subject: [PATCH 1/4] Jacobi tests pass --- REQUIRE | 2 +- src/besselroots.jl | 100 ++++++++++++++++++------------------ src/gaussjacobi.jl | 76 +++++++++++++-------------- src/gausslegendre.jl | 12 ++--- test/test_gausschebyshev.jl | 14 ++--- test/test_gaussjacobi.jl | 20 ++++---- test/test_gausslegendre.jl | 22 ++++---- test/test_gausslobatto.jl | 22 ++++---- test/test_gaussradau.jl | 14 ++--- 9 files changed, 140 insertions(+), 142 deletions(-) diff --git a/REQUIRE b/REQUIRE index 98f71e9..84feefe 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,2 +1,2 @@ -julia 0.4 +julia 0.5 Compat 0.8.8 diff --git a/src/besselroots.jl b/src/besselroots.jl index 0cc012c..c2faa1c 100644 --- a/src/besselroots.jl +++ b/src/besselroots.jl @@ -1,4 +1,4 @@ -function besselroots(nu::Float64, n::Int64) +function besselroots(nu::Float64, n::Int64) #BESSELROOTS The first N roots of the function J_v(x) # DEVELOPERS NOTES: @@ -10,9 +10,9 @@ function besselroots(nu::Float64, n::Int64) # V > 5 --> moderately accurate for the 6 first zeros and good # approximations for the others (McMahon's expansion) -# This code was originally written by L. L. Peixoto in MATLAB. +# This code was originally written by L. L. Peixoto in MATLAB. # Later modified by A. Townsend to work in Julia - + if ( n == 0 ) x = [] elseif( n > 0 && nu == 0 ) @@ -20,7 +20,7 @@ elseif( n > 0 && nu == 0 ) correctFirstFew = Tabulate( ) x[1:min(n,20)] = correctFirstFew[1:min(n,20)] x -elseif ( n>0 && nu >= -1 && nu <= 5 ) +elseif ( n>0 && nu >= -1 && nu <= 5 ) x = McMahon( nu, n ) correctFirstFew = Piessens( nu ) x[1:min(n,6)] = correctFirstFew[1:min(n,6)] @@ -33,55 +33,55 @@ end function McMahon( nu::Float64, n::Int64 ) -# McMahon's expansion. This expansion gives very accurate approximation -# for the sth zero (s >= 7) in the whole region V >=- 1, and moderate -# approximation in other cases. -b = Array(Float64,n) -mu = 4*nu^2 -a1 = 1. / 8. -a3 = (7.*mu-31.) / 384. -a5 = 4.*(3779.+mu*(-982.+83.*mu)) / 61440. # Evaluate via Horner's method. -a7 = 6.*(-6277237.+mu*(1585743.+mu*(-153855.+6949.*mu))) / 20643840.; -a9 = 144.*(2092163573.+mu*(-512062548.+mu*(48010494.+mu*(-2479316.+70197.*mu)))) / 11890851840. -a11 = 720.*(-8249725736393.+mu*(1982611456181.+mu*(-179289628602.+mu*(8903961290. + - mu*(-287149133.+5592657.*mu))))) / 10463949619200. -a13 = 576.*(423748443625564327. + mu*(-100847472093088506.+mu*(8929489333108377. + - mu*(-426353946885548.+mu*(13172003634537.+mu*(-291245357370. + mu*4148944183.)))))) / 13059009124761600. -for k=1:n b[k] = .25*(2*nu+4*k-1)*pi end # beta -# Evaluate using Horner's scheme: -x = b - (mu-1)*( ((((((a13./b.^2 + a11)./b.^2 + a9)./b.^2 + a7)./b.^2 + a5)./b.^2 + a3)./b.^2 + a1)./b) + # McMahon's expansion. This expansion gives very accurate approximation + # for the sth zero (s >= 7) in the whole region V >=- 1, and moderate + # approximation in other cases. + b = Array{Float64}(n) + mu = 4*nu^2 + a1 = 1. / 8. + a3 = (7.*mu-31.) / 384. + a5 = 4.*(3779.+mu*(-982.+83.*mu)) / 61440. # Evaluate via Horner's method. + a7 = 6.*(-6277237.+mu*(1585743.+mu*(-153855.+6949.*mu))) / 20643840.; + a9 = 144.*(2092163573.+mu*(-512062548.+mu*(48010494.+mu*(-2479316.+70197.*mu)))) / 11890851840. + a11 = 720.*(-8249725736393.+mu*(1982611456181.+mu*(-179289628602.+mu*(8903961290. + + mu*(-287149133.+5592657.*mu))))) / 10463949619200. + a13 = 576.*(423748443625564327. + mu*(-100847472093088506.+mu*(8929489333108377. + + mu*(-426353946885548.+mu*(13172003634537.+mu*(-291245357370. + mu*4148944183.)))))) / 13059009124761600. + for k=1:n b[k] = .25*(2*nu+4*k-1)*pi end # beta + # Evaluate using Horner's scheme: + x = b - (mu-1)*( ((((((a13./b.^2 + a11)./b.^2 + a9)./b.^2 + a7)./b.^2 + a5)./b.^2 + a3)./b.^2 + a1)./b) end function Tabulate( ) -# First 20 roots of J0(x) are precomputed (using Wolfram Alpha): -y = Array(Float64,20) -y = [ 2.4048255576957728 - 5.5200781102863106 - 8.6537279129110122 - 11.791534439014281 - 14.930917708487785 - 18.071063967910922 - 21.211636629879258 - 24.352471530749302 - 27.493479132040254 - 30.634606468431975 - 33.775820213573568 - 36.917098353664044 - 40.058425764628239 - 43.199791713176730 - 46.341188371661814 - 49.482609897397817 - 52.624051841114996 - 55.765510755019979 - 58.906983926080942 - 62.048469190227170 ] + # First 20 roots of J0(x) are precomputed (using Wolfram Alpha): + y = Array{Float64}(20) + y = [ 2.4048255576957728 + 5.5200781102863106 + 8.6537279129110122 + 11.791534439014281 + 14.930917708487785 + 18.071063967910922 + 21.211636629879258 + 24.352471530749302 + 27.493479132040254 + 30.634606468431975 + 33.775820213573568 + 36.917098353664044 + 40.058425764628239 + 43.199791713176730 + 46.341188371661814 + 49.482609897397817 + 52.624051841114996 + 55.765510755019979 + 58.906983926080942 + 62.048469190227170 ] end -function Piessens( nu::Float64 ) -# Piessens's Chebyshev series approximations (1984). Calculates the 6 first -# zeros to at least 12 decimal figures in region -1 <= V <= 5: - y = Array(Float64,6); +function Piessens( nu::Float64 ) + # Piessens's Chebyshev series approximations (1984). Calculates the 6 first + # zeros to at least 12 decimal figures in region -1 <= V <= 5: + y = Array{Float64}(6); C = [ 2.883975316228 8.263194332307 11.493871452173 14.689036505931 17.866882871378 21.034784308088 @@ -114,8 +114,8 @@ function Piessens( nu::Float64 ) 0.000000000008 0.000000000011 0. 0. 0. 0. -0.000000000003 -0.000000000005 0. 0. 0. 0. 0.000000000001 0.000000000002 0. 0. 0. 0.] - - T = Array(Float64,size(C,1)) + + T = Array{Float64}(size(C,1)) pt = (nu-2)/3 T[1], T[2] = 1., pt for k = 2:size(C,1)-1 @@ -124,4 +124,4 @@ function Piessens( nu::Float64 ) y[1:6] = T'*C; y[1] *= sqrt(nu+1); # Scale the first root. y -end \ No newline at end of file +end diff --git a/src/gaussjacobi.jl b/src/gaussjacobi.jl index 042dceb..dd737c8 100644 --- a/src/gaussjacobi.jl +++ b/src/gaussjacobi.jl @@ -38,8 +38,8 @@ function JacobiRec(n::Int, a::Float64, b::Float64) x11, x12 = HalfRec(n, a, b, 1) x21, x22 = HalfRec(n, b, a, 0) - x = Array(Float64, n) - w = Array(Float64, n) + x = Array{Float64}(n) + w = Array{Float64}(n) m1 = length(x11) m2 = length(x21) sum_w = 0.0 @@ -76,15 +76,15 @@ function HalfRec(n::Int, a::Float64, b::Float64, flag) a1 = .25 - a^2 b1 = .25 - b^2 c1² = c1^2 - x = Array(Float64, m) + x = Array{Float64}(m) @inbounds for i in 1:m C = muladd(2.0, r[i], a - .5) * (π * c1) C_2 = 0.5 * C x[i] = cos(muladd(c1², muladd(-b1, tan(C_2), a1 * cot(C_2)), C)) end - P1 = Array(Float64, m) - P2 = Array(Float64, m) + P1 = Array{Float64}(m) + P2 = Array{Float64}(m) # Loop until convergence: for _ in 1:10 innerJacobiRec!(n, x, a, b, P1, P2) @@ -133,8 +133,8 @@ end function innerJacobiRec(n, x, a, b) # EVALUATE JACOBI POLYNOMIALS AND ITS DERIVATIVE USING THREE-TERM RECURRENCE. N = length(x) - P = Array(Float64, N) - PP = Array(Float64, N) + P = Array{Float64}(N) + PP = Array{Float64}(N) innerJacobiRec!(n, x, a, b, P, PP) P, PP end @@ -182,8 +182,8 @@ function asy1(n::Int, a::Float64, b::Float64, nbdy) # Algorithm for computing nodes and weights in the interior. # Approximate roots via asymptotic formula: (Gatteschi and Pittaluga, 1985) - K = (2*(n:-1:1)+a-.5)*pi/(2*n+a+b+1) - tt = K + 1/(2*n+a+b+1)^2*((.25-a^2)*cot(.5*K)-(.25-b^2)*tan(.5*K)) + K = (2*(n:-1:1).+a.-0.5).*pi./(2*n+a+b+1) + tt = K .+ (1/(2*n+a+b+1)^2).*((0.25-a^2).*cot.(0.5.*K).-(0.25-b^2).*tan.(0.5*K)) # First half (x > 0): t = tt[tt .<= pi/2]' @@ -204,7 +204,7 @@ function asy1(n::Int, a::Float64, b::Float64, nbdy) t += vals[1]./vals[2] # Newton update. # Store: - x = cos(t) + x = cos.(t) w = 1./vals[2].^2 # Second half (x < 0): @@ -226,7 +226,7 @@ function asy1(n::Int, a::Float64, b::Float64, nbdy) t += vals[1]./vals[2] # Newton update. # Store: - [-cos(vec(t));vec(x)],[1./vec(vals[2]).^2;vec(w)] + [.-(cos.(vec(t)));vec(x)],[1./vec(vals[2]).^2;vec(w)] end function feval_asy1(n::Int, a::Float64, b::Float64, t, idx) @@ -236,23 +236,23 @@ function feval_asy1(n::Int, a::Float64, b::Float64, t, idx) M = 20 # Some often used vectors/matrices: - onesT = ones(1,length(t)); onesM = ones(M); MM = (collect(0:M-1)')'; + onesT = ones(length(t))'; onesM = ones(M); MM = collect(0:M-1); # The sine and cosine terms: - alpha = (.5*(2*n+a+b+1+MM))*onesT .* (onesM*t) - .5*(a+.5)*pi - cosA = cos(alpha); sinA = sin(alpha) + alpha = (0.5*(2*n+a+b+1+MM))*onesT .* (onesM*t) - 0.5*(a+0.5)*pi + cosA = cos.(alpha); sinA = sin.(alpha) - sinT = onesM*sin(t) - cosT = onesM*cos(t) - cosA2 = cosA.*cosT + sinA.*sinT - sinA2 = sinA.*cosT - cosA.*sinT + sinT = onesM*sin.(t) + cosT = onesM*cos.(t) + cosA2 = cosA.*cosT .+ sinA.*sinT + sinA2 = sinA.*cosT .- cosA.*sinT one = ones(t) - sinT = vcat( one , cumprod(onesM[2:end]*(.5*csc(.5*t)))) - cosT = .5*sec(.5*t) + sinT = vcat( one , cumprod(onesM[2:end].*(0.5.*csc.(0.5.*t)))) + cosT = 0.5.*sec.(0.5.*t) j = transpose(0:M-2) - vec = (.5+a+j).*(.5-a+j)./(j+1)./(2*n+a+b+j+2) + vec = (0.5.+a.+j).*(0.5.-a.+j)./(j.+1)./(2*n.+a.+b.+j.+2) P1 = [1 cumprod(vec,2)] P1[3:4:end] = -P1[3:4:end] P1[4:4:end] = -P1[4:4:end] @@ -281,25 +281,25 @@ function feval_asy1(n::Int, a::Float64, b::Float64, t, idx) SC = sinT for m = 1:M l = 1:2:m - phi = PHI[m:m, l] - dS1 = phi * SC[l, :] .* cosA[m:m, :] - phi2 = PHI2[m:m, l] - dS12 = phi2*SC[l, :] .* cosA2[m:m, :] + phi = PHI[m, l]' + dS1 = phi * SC[l, :] .* cosA[m, :]' + phi2 = PHI2[m, l]' + dS12 = phi2*SC[l, :] .* cosA2[m, :]' l = 2:2:m - phi = PHI[m:m, l] - dS2 = phi * SC[l, :] .* sinA[m:m, :] - phi2 = PHI2[m:m, l] - dS22 = phi2 * SC[l, :] .* sinA2[m:m, :] + phi = PHI[m, l]' + dS2 = phi * SC[l, :] .* sinA[m, :]' + phi2 = PHI2[m, l]' + dS22 = phi2 * SC[l, :] .* sinA2[m, :]' if m - 1 > 10 && norm(dS1[idx] + dS2[idx], Inf) < eps(Float64) / 100 break end S = S + dS1 + dS2 S2 = S2 + dS12 + dS22 - SC[1:m,:] = broadcast(.*,SC[1:m,:],cosT) + SC[1:m,:] = SC[1:m,:].*cosT end # Constant out the front: - dsa = .5*(a^2)/n; dsb = .5*(b^2)/n; dsab = .25*(a+b)^2/n + dsa = 0.5*(a^2)/n; dsb = 0.5*(b^2)/n; dsab = 0.25*(a+b)^2/n ds = dsa + dsb - dsab; s = ds; j = 1 dsold = ds # to fix a = -b bug. while ( (abs(ds/s) + dsold) > eps(Float64)/10 ) @@ -324,9 +324,9 @@ function feval_asy1(n::Int, a::Float64, b::Float64, t, idx) S2 = C2*S2 # Use relation for derivative: - ders = (n*(a-b-(2*n+a+b)*cos(t)).*vals + 2*(n+a)*(n+b)*S2)/(2*n+a+b)./sin(t) - t = abs( t ) - denom = 1./real(sin(t/2).^(a+.5).*cos(t/2).^(b+.5)) + ders = (n.*(a.-b.-(2*n+a+b).*cos.(t)).*vals .+ 2.*(n+a).*(n+b).*S2)./(2*n+a+b)./sin.(t) + t .= abs.( t ) + denom = 1./real.(sin.(t/2).^(a+0.5).*cos.(t./2).^(b+0.5)) vals = vals.*denom ders = ders.*denom @@ -342,7 +342,7 @@ function boundary(n::Int, a::Float64, b::Float64, npts) # Approximate roots via asymptotic formula: (see Olver 1974) phik = jk/(n + .5*(a + b + 1)) - x = cos( phik + ((a^2-.25)*(1-phik.*cot(phik))./(8*phik) - .25*(a^2-b^2)*tan(.5*phik))/(n + .5*(a + b + 1))^2 ) + x = cos.( phik .+ ((a^2-0.25).*(1.-phik.*cot.(phik))./(8*phik) .- 0.25.*(a^2-b^2).*tan.(0.5.*phik))./(n + 0.5*(a + b + 1))^2 ) dx = 1.0; counter = 0; # Newton iteration: @@ -375,11 +375,11 @@ function JacobiGW( n::Int64, a::Float64, b::Float64 ) (b^2 - a^2)./((abi - 2).*abi) (b^2 - a^2)./((2*n - 2+ab).*(2*n+ab))] bb = [2*sqrt( (1 + a)*(1 + b)/(ab + 3))/(ab + 2) ; - 2*sqrt(ii.*(ii + a).*(ii + b).*(ii + ab)./(abi.^2 - 1))./abi] + 2.*sqrt.(ii.*(ii .+ a).*(ii .+ b).*(ii .+ ab)./(abi.^2 .- 1))./abi] TT = diagm(bb,1) + diagm(aa) + diagm(bb,-1) # Jacobi matrix. x, V = eig( TT ) # Eigenvalue decomposition. # Quadrature weights: - w = V[1,:].^2*( 2^(ab+1)*gamma(a+1)*gamma(b+1)/gamma(2+ab) ); - w = w/sum(w); + w = V[1,:].^2.*( 2^(ab+1)*gamma(a+1)*gamma(b+1)/gamma(2+ab) ); + w .= w./sum(w); x, vec(w) end diff --git a/src/gausslegendre.jl b/src/gausslegendre.jl index 490550f..f7fc140 100644 --- a/src/gausslegendre.jl +++ b/src/gausslegendre.jl @@ -57,7 +57,7 @@ function legpts_nodes(n, a) # ASYMPTOTIC EXPANSION FOR THE GAUSS-LEGENDRE NODES. vn = 1 / (n + 0.5) m = length(a) - nodes = cot(a) + nodes = cot.(a) vn² = vn * vn vn⁴ = vn² * vn² @inbounds if n <= 255 @@ -109,7 +109,7 @@ function legpts_weights(n, m, a) # ASYMPTOTIC EXPANSION FOR THE GAUSS-LEGENDRE WEIGHTS. vn = 1 / (n + 0.5) vn² = vn^2 - weights = Array(Float64, m) + weights = Array{Float64}(m) if n <= 850000 @inbounds for i in 1:m weights[i] = cot(a[i]) @@ -209,8 +209,8 @@ end function innerRec(n, x) # EVALUATE LEGENDRE AND ITS DERIVATIVE USING THREE-TERM RECURRENCE RELATION. N = size(x, 1) - myPm1 = Array(Float64, N) - myPPm1 = Array(Float64, N) + myPm1 = Array{Float64}(N) + myPPm1 = Array{Float64}(N) @inbounds for j = 1:N xj = x[j] Pm2 = 1.0 @@ -242,7 +242,7 @@ const besselZeros_20 = [2.4048255576957728, 5.5200781102863106, function besselZeroRoots(m) # BESSEL0ROOTS ROOTS OF BESSELJ(0,x). USE ASYMPTOTICS. # Use McMahon's expansion for the remainder (NIST, 10.21.19): - jk = Array(Float64, m) + jk = Array{Float64}(m) p = (1071187749376 / 315, 0.0, -401743168 / 105, 0.0, 120928 / 15, 0.0, -124 / 3, 0.0, 1.0, 0.0) # First 20 are precomputed: @@ -281,7 +281,7 @@ function besselJ1(m) # BESSELJ1 EVALUATE BESSELJ(1,x)^2 AT ROOTS OF BESSELJ(0,x). # USE ASYMPTOTICS. Use Taylor series of (NIST, 10.17.3) and McMahon's # expansion (NIST, 10.21.19): - Jk2 = Array(Float64, m) + Jk2 = Array{Float64}(m) c = (-171497088497 / 15206400, 461797 / 1152, -172913 / 8064, 151 / 80, -7 / 24, 0.0, 2.0) # First 10 are precomputed: diff --git a/test/test_gausschebyshev.jl b/test/test_gausschebyshev.jl index ab9c788..79c4239 100644 --- a/test/test_gausschebyshev.jl +++ b/test/test_gausschebyshev.jl @@ -1,16 +1,16 @@ # Test for gausschebyshev(). -n = 10 +n = 10 # Test x.^2: x, w = gausschebyshev(n,1) -@test_approx_eq dot(x.^2,w) pi/2 +@test dot(x.^2,w) ≈ pi/2 x, w = gausschebyshev(n,2) -@test_approx_eq dot(x.^2,w) pi/8 +@test dot(x.^2,w) ≈ pi/8 x, w = gausschebyshev(n,3) -@test_approx_eq dot(x.^2,w) pi/2 +@test dot(x.^2,w) ≈ pi/2 x, w = gausschebyshev(n,4) -@test_approx_eq dot(x.^2,w) pi/2 +@test dot(x.^2,w) ≈ pi/2 # Test x^3: x, w = gausschebyshev(n,1) @@ -18,6 +18,6 @@ x, w = gausschebyshev(n,1) x, w = gausschebyshev(n,2) @test abs(dot(x.^3,w))<1e-15 x, w = gausschebyshev(n,3) -@test_approx_eq dot(x.^3,w) 3*pi/8 +@test dot(x.^3,w) ≈ 3*pi/8 x, w = gausschebyshev(n,4) -@test_approx_eq dot(x.^3,w) -3*pi/8 \ No newline at end of file +@test dot(x.^3,w) ≈ -3*pi/8 diff --git a/test/test_gaussjacobi.jl b/test/test_gaussjacobi.jl index d0d33b5..5d503b8 100644 --- a/test/test_gaussjacobi.jl +++ b/test/test_gaussjacobi.jl @@ -9,7 +9,7 @@ x, w = gaussjacobi(n, a, b) @test abs(x[37] - 0.912883347814032) < tol @test abs(w[37] - 0.046661910947553) < tol -@test_approx_eq dot(w,exp(x)) 2.7568520361985516 +@test dot(w,exp.(x)) ≈ 2.7568520361985516 # Test a larger n (using ASY) a = -.7 @@ -20,7 +20,7 @@ x, w = gaussjacobi(n, a, b) @test abs(x[37] + 0.893103435898983) < tol @test abs(w[37] - 1.962534523788093e-04) < tol -@test_approx_eq dot(w,exp(x)) 16.722957039404044 +@test dot(w,exp.(x)) ≈ 16.722957039404044 # Test n = 1: a = 1.0 @@ -29,7 +29,7 @@ x, w = gaussjacobi(1, a, b) @test abs(x[1] - (b - a) / (a + b + 2)) < tol @test abs(w[1] - 2^(a + b + 1) * beta(a + 1, b + 1)) < tol -@test_approx_eq dot(w,ones(x)) 1.3333333333333333 +@test dot(w,ones(x)) ≈ 1.3333333333333333 x, w = gaussjacobi(1013, .9, -.1) @test abs(x[2] + 0.999986012231899) < tol @@ -37,7 +37,7 @@ x, w = gaussjacobi(1013, .9, -.1) @test abs(w[2] - 9.314674169892358e-05) < tol @test abs(w[13] - 4.654651764553262e-04) < tol -@test_approx_eq dot(w,exp(x)) 1.6915068974063106 +@test dot(w,exp.(x)) ≈ 1.6915068974063106 x, w = gaussjacobi(10013, .9, -.1) @@ -46,21 +46,21 @@ x, w = gaussjacobi(10013, .9, -.1) @test abs(w[2] - 1.509654630405615e-06) < tol @test abs(w[13] - 7.548275262993863e-06) < tol -@test_approx_eq dot(w,exp(x)) 1.6915068974063106 +@test dot(w,exp.(x)) ≈ 1.6915068974063106 # tests bug where row vectors were returned @test isa(x,Vector{Float64}) @test isa(w,Vector{Float64}) -# test last alpha and beta parameters: +# test last alpha and beta parameters: x, w = gaussjacobi(100, 19., 21.) -@test abs(x[87] - 0.832211446176040) < tol +@test abs(x[87] - 0.832211446176040) < tol @test abs(w[50] - 0.064530500882703) < tol -# test for small alpha and beta: +# test for small alpha and beta: x, w = gaussjacobi(10000, .1, .2) -@test abs(x[1] - -0.999999963363548) < tol +@test abs(x[1] - -0.999999963363548) < tol @test abs(w[500] - 2.183393039546711e-05) < tol @@ -69,5 +69,3 @@ x,w=gaussjacobi(20,11.,0.) # tests bug where row vectors were returned @test isa(x,Vector{Float64}) @test isa(w,Vector{Float64}) - - diff --git a/test/test_gausslegendre.jl b/test/test_gausslegendre.jl index 96ce9a2..fc2d77d 100644 --- a/test/test_gausslegendre.jl +++ b/test/test_gausslegendre.jl @@ -3,18 +3,18 @@ for n=2:10 # check all special cases x, w = gausslegendre(n) - @test_approx_eq dot( w,(x.^2)) 2/3 + @test dot( w,(x.^2)) ≈ 2/3 end tol = 1e-14 n = 42 x, w = gausslegendre(n) @test length(x) == n && length(w) == n -@test_approx_eq_eps x[37] 0.910959724904127 tol -@test_approx_eq_eps w[37] 0.030479240699603 tol +@test x[37] ≈ 0.910959724904127 atol=tol +@test w[37] ≈ 0.030479240699603 atol=tol -@test_approx_eq dot( w,(x.^2)) 2/3 -@test_approx_eq dot( w,exp(x)) exp(1)-exp(-1) +@test dot( w,(x.^2)) ≈ 2/3 +@test dot( w,exp.(x)) ≈ exp(1)-exp(-1) # Test a larger n (using ASY) n = 251 @@ -23,8 +23,8 @@ x, w = gausslegendre(n) @test abs(x[37] + 0.896467746955729) < tol @test abs(w[37] - 0.005535005742012) < tol -@test_approx_eq dot( w,(x.^2)) 2/3 -@test_approx_eq dot( w,exp(x)) exp(1)-exp(-1) +@test dot( w,(x.^2)) ≈ 2/3 +@test dot( w,exp.(x)) ≈ exp(1)-exp(-1) x, w = gausslegendre(1013) @test norm(x[2] - -0.999985167586110, Inf) < tol @@ -32,8 +32,8 @@ x, w = gausslegendre(1013) @test norm(w[2] - 1.681691163200592e-05, Inf) < tol @test norm(w[13] - 1.224755309137936e-04, Inf) < tol -@test_approx_eq dot( w,(x.^2)) 2/3 -@test_approx_eq dot( w,exp(x)) exp(1)-exp(-1) +@test dot( w,(x.^2)) ≈ 2/3 +@test dot( w,exp.(x)) ≈ exp(1)-exp(-1) x, w = gausslegendre(10013) @test norm(x[2] - -0.999999848054223, Inf) < tol @@ -41,5 +41,5 @@ x, w = gausslegendre(10013) @test norm(w[2] - 1.722757320118474e-07, Inf) < tol @test norm(w[13] - 1.254980540032470e-06, Inf) < tol -@test_approx_eq dot( w,(x.^2)) 2/3 -@test_approx_eq dot( w,exp(x)) exp(1)-exp(-1) +@test dot( w,(x.^2)) ≈ 2/3 +@test dot( w,exp.(x)) ≈ exp(1)-exp(-1) diff --git a/test/test_gausslobatto.jl b/test/test_gausslobatto.jl index 037e914..995ebf5 100644 --- a/test/test_gausslobatto.jl +++ b/test/test_gausslobatto.jl @@ -1,19 +1,19 @@ # Test for GaussLobatto() n = 2 x,w = gausslobatto(n) -@test_approx_eq x[1] -1. -@test_approx_eq x[2] 1. -@test_approx_eq w[1] 1. -@test_approx_eq w[2] 1. +@test x[1] ≈ -1. +@test x[2] ≈ 1. +@test w[1] ≈ 1. +@test w[2] ≈ 1. n = 3 x,w = gausslobatto(n) -@test_approx_eq x[1] -1. +@test x[1] ≈ -1. @test abs(x[2])<1e-15 -@test_approx_eq x[3] 1. -@test_approx_eq w[1] 1./3 -@test_approx_eq w[2] 4./3 -@test_approx_eq w[3] 1./3 +@test x[3] ≈ 1. +@test w[1] ≈ 1./3 +@test w[2] ≈ 4./3 +@test w[3] ≈ 1./3 tol = 1e-14 n = 42 @@ -23,5 +23,5 @@ x,w = gausslobatto(n) @test abs(w[37] - 0.029306411216166) < tol @test ( x[1] == -1 && x[n] == 1 ) -@test_approx_eq dot( w,(x.^2)) 2/3 -@test_approx_eq dot( w,exp(x)) exp(1)-exp(-1) +@test dot( w,(x.^2)) ≈ 2/3 +@test dot( w,exp(x)) ≈ exp(1)-exp(-1) diff --git a/test/test_gaussradau.jl b/test/test_gaussradau.jl index 629c1aa..98099ea 100644 --- a/test/test_gaussradau.jl +++ b/test/test_gaussradau.jl @@ -2,15 +2,15 @@ n = 1 x, w = gaussradau(n) -@test_approx_eq x[1] -1. -@test_approx_eq w[1] 2. +@test x[1] ≈ -1. +@test w[1] ≈ 2. n = 2 x, w = gaussradau(n) -@test_approx_eq x[1] -1. -@test_approx_eq x[2] 1./3. -@test_approx_eq w[1] .5 -@test_approx_eq w[2] 1.5 +@test x[1] ≈ -1. +@test x[2] ≈ 1./3. +@test w[1] ≈ .5 +@test w[2] ≈ 1.5 tol = 1e-14 n = 42 @@ -21,4 +21,4 @@ x, w = gaussradau(n) @test abs(w[37] - 0.031190846817016) < tol @test x[1] == -1 -@test_approx_eq dot( w,exp(x)) exp(1)-exp(-1) \ No newline at end of file +@test dot( w,exp(x)) ≈ exp(1)-exp(-1) \ No newline at end of file From 6cec4a1d9867295312b113ba9f5b233c3b94091b Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Thu, 2 Feb 2017 22:08:35 +1100 Subject: [PATCH 2/4] three tests pass --- test/test_besselroots.jl | 8 ++++---- test/test_gausslobatto.jl | 2 +- test/test_gaussradau.jl | 6 +++--- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/test/test_besselroots.jl b/test/test_besselroots.jl index 1277b6a..aad8bce 100644 --- a/test/test_besselroots.jl +++ b/test/test_besselroots.jl @@ -1,9 +1,9 @@ tol = 1e-11 -# Check if besselj(nu, besselroots(nu, n) ) is small: +# Check if besselj(nu, besselroots(nu, n) ) is small: -for nu = 0.:0.1:5. +for nu = 0.:0.1:5. n = 10 - @test norm( besselj(nu, besselroots(nu, n) ), Inf ) < tol + @test norm( besselj.(nu, besselroots(nu, n) ), Inf ) < tol -end \ No newline at end of file +end diff --git a/test/test_gausslobatto.jl b/test/test_gausslobatto.jl index 995ebf5..fe539e7 100644 --- a/test/test_gausslobatto.jl +++ b/test/test_gausslobatto.jl @@ -24,4 +24,4 @@ x,w = gausslobatto(n) @test ( x[1] == -1 && x[n] == 1 ) @test dot( w,(x.^2)) ≈ 2/3 -@test dot( w,exp(x)) ≈ exp(1)-exp(-1) +@test dot( w,exp.(x)) ≈ exp(1)-exp(-1) diff --git a/test/test_gaussradau.jl b/test/test_gaussradau.jl index 98099ea..cc4cfa0 100644 --- a/test/test_gaussradau.jl +++ b/test/test_gaussradau.jl @@ -1,11 +1,11 @@ # Test for GaussRadau() -n = 1 +n = 1 x, w = gaussradau(n) @test x[1] ≈ -1. @test w[1] ≈ 2. -n = 2 +n = 2 x, w = gaussradau(n) @test x[1] ≈ -1. @test x[2] ≈ 1./3. @@ -21,4 +21,4 @@ x, w = gaussradau(n) @test abs(w[37] - 0.031190846817016) < tol @test x[1] == -1 -@test dot( w,exp(x)) ≈ exp(1)-exp(-1) \ No newline at end of file +@test dot( w,exp.(x)) ≈ exp(1)-exp(-1) From d4fb98c868ce9fc8cc7ace6c4b49ef2c77f0c773 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Fri, 3 Feb 2017 16:17:15 +1100 Subject: [PATCH 3/4] 0.5 and 0.6 pass tests without errors! --- .travis.yml | 1 - src/gausshermite.jl | 58 +++++++++++++++++++------------------- src/gaussjacobi.jl | 2 +- src/gausslaguerre.jl | 10 +++---- test/test_gausshermite.jl | 10 +++++-- test/test_gausslaguerre.jl | 14 +++++++-- test/test_gausslegendre.jl | 6 ++-- 7 files changed, 59 insertions(+), 42 deletions(-) diff --git a/.travis.yml b/.travis.yml index ea2973f..7354c42 100644 --- a/.travis.yml +++ b/.travis.yml @@ -3,7 +3,6 @@ os: - linux - osx julia: - - 0.4 - 0.5 - nightly notifications: diff --git a/src/gausshermite.jl b/src/gausshermite.jl index 1608cf0..8172526 100644 --- a/src/gausshermite.jl +++ b/src/gausshermite.jl @@ -39,20 +39,20 @@ function hermpts_asy( n::Int ) x0 = HermiteInitialGuesses( n ) # get initial guesses t0 = x0./sqrt(2n+1) -theta0 = acos(t0) # convert to theta-variable +theta0 = acos.(t0) # convert to theta-variable val = x0; for k = 1:20 val = hermpoly_asy_airy(n, theta0); - dt = -val[1]./(sqrt(2)*sqrt(2n+1)*val[2].*sin(theta0)) - theta0 = theta0 - dt; # Newton update + dt = -val[1]./(sqrt(2).*sqrt(2n+1).*val[2].*sin.(theta0)) + theta0 .-= dt; # Newton update if norm(dt,Inf) < sqrt(eps(Float64))/10 break end end -t0 = cos(theta0) +t0 = cos.(theta0) x = sqrt(2n+1)*t0 #back to x-variable -ders = x.*val[1] + sqrt(2)*val[2] -w = (exp(-x.^2)./ders.^2)'; # quadrature weights +ders = x.*val[1] .+ sqrt(2).*val[2] +w = exp.(-x.^2)./ders.^2; # quadrature weights x = (x, w) end @@ -66,14 +66,14 @@ val = x0 for kk = 1:10 val = hermpoly_rec(n, x0) dx = val[1]./val[2] - dx[ isnan( dx ) ] = 0 + dx[ isnan.( dx ) ] = 0 x0 = x0 - dx if norm(dx, Inf)=20 - +# Error if n < 20 because initial guesses are based on asymptotic expansions: +@assert n>=20 + # Gatteschi formula involving airy roots [1]. # These initial guess are good near x = sqrt(n+1/2); if mod(n,2) == 1 @@ -209,9 +209,9 @@ airyrts_exact = [-2.338107410459762 # Exact Airy roots. -12.828776752865757] airyrts[1:10] = airyrts_exact # correct first 10. -x_init = sqrt(abs(nu + 2^(2/3)*airyrts*nu^(1/3) + 1/5*2^(4/3)*airyrts.^2*nu^(-1/3) + - (11/35-a^2-12/175*airyrts.^3)/nu + (16/1575*airyrts+92/7875*airyrts.^4)*2^(2/3)*nu^(-5/3) - - (15152/3031875*airyrts.^5+1088/121275*airyrts.^2)*2^(1/3)*nu^(-7/3))) +x_init = sqrt.(abs.(nu .+ (2^(2/3)).*airyrts.*nu^(1/3) .+ (1/5*2^(4/3)).*airyrts.^2.*nu^(-1/3) .+ + (11/35-a^2-12/175).*airyrts.^3./nu .+ ((16/1575).*airyrts.+(92/7875).*airyrts.^4).*2^(2/3).*nu^(-5/3) .- + ((15152/3031875).*airyrts.^5.+(1088/121275).*airyrts.^2).*2^(1/3).*nu^(-7/3))) x_init_airy = real( flipdim(x_init,1) ) # Tricomi initial guesses. Equation (2.1) in [1]. Originally in [2]. @@ -223,14 +223,14 @@ nu = (4*m+2*a+2) rhs = (4*m-4*collect(1:m)+3)./nu*pi for k = 1:7 - val = Tnk0 - sin(Tnk0) - rhs - dval = 1 - cos(Tnk0) + val = Tnk0 .- sin.(Tnk0) .- rhs + dval = 1 .- cos.(Tnk0) dTnk0 = val./dval - Tnk0 = Tnk0 - dTnk0 + Tnk0 = Tnk0 .- dTnk0 end -tnk = cos(Tnk0/2).^2 -x_init_sin = sqrt(nu*tnk - (5./(4*(1-tnk).^2) - 1./(1-tnk)-1+3*a^2)/3/nu) +tnk = cos.(Tnk0./2).^2 +x_init_sin = sqrt.(nu*tnk .- (5./(4.*(1-tnk).^2) .- 1./(1.-tnk).-1.+3*a^2)./3./nu) # Patch together p = 0.4985+eps(Float64) @@ -251,7 +251,7 @@ end function hermpts_gw( n::Int ) # Golub--Welsch algorithm. Used here for n<=20. - beta = sqrt(.5*(1:n-1)) # 3-term recurrence coeffs + beta = sqrt.(0.5.*(1:n-1)) # 3-term recurrence coeffs T = diagm(beta, 1) + diagm(beta, -1) # Jacobi matrix (D, V) = eig(T) # Eigenvalue decomposition indx = sortperm(D) # Hermite points @@ -263,4 +263,4 @@ function hermpts_gw( n::Int ) x = x[ii] w = w[ii] return (x,w) -end \ No newline at end of file +end diff --git a/src/gaussjacobi.jl b/src/gaussjacobi.jl index dd737c8..f97c6c3 100644 --- a/src/gaussjacobi.jl +++ b/src/gaussjacobi.jl @@ -226,7 +226,7 @@ function asy1(n::Int, a::Float64, b::Float64, nbdy) t += vals[1]./vals[2] # Newton update. # Store: - [.-(cos.(vec(t)));vec(x)],[1./vec(vals[2]).^2;vec(w)] + [(-).(cos.(vec(t)));vec(x)],[1./vec(vals[2]).^2;vec(w)] end function feval_asy1(n::Int, a::Float64, b::Float64, t, idx) diff --git a/src/gausslaguerre.jl b/src/gausslaguerre.jl index d1990c5..29f949f 100644 --- a/src/gausslaguerre.jl +++ b/src/gausslaguerre.jl @@ -37,7 +37,7 @@ function laguerreGW( n::Int64, alpha::Float64 ) # Calculate Gauss-Laguerre nodes and weights based on Golub-Welsch alph = 2*(1:n)-1+alpha # 3-term recurrence coeffs - beta = sqrt( (1:n-1).*(alpha + (1:n-1) ) ) + beta = sqrt.( (1:n-1).*(alpha .+ (1:n-1) ) ) T = diagm(beta,1) + diagm(alph) + diagm(beta,-1) # Jacobi matrix x, V = eig( T ) # eigenvalue decomposition w = gamma(alpha+1)*V[1,:].^2 # Quadrature weights @@ -766,7 +766,7 @@ function asyAiry(np, y, alpha, T) fn = (np*3im*( sqrt(z)*sqrt(1 - z) - acos(sqrt(z) ) ) )^(2/3) d = z - 1 if T == 1 - return real( 4*sqrt(pi)/z^(1/4)/(d + 0im)^(1/4)*z^(-alpha/2)*(cos( (alpha + 1)/2*acos(2*z - 1) )*fn^(1/4)*airy(0,fn) + -1im*sin( (alpha + 1)/2*acos(2*z - 1) )*ifelse(angle(z-1) <= 0, -one(z), one(z) )*fn^(-1/4)*airy(1,fn) ) ) + return real( 4*sqrt(pi)/z^(1/4)/(d + 0im)^(1/4)*z^(-alpha/2)*(cos( (alpha + 1)/2*acos(2*z - 1) )*fn^(1/4)*airyai(fn) + -1im*sin( (alpha + 1)/2*acos(2*z - 1) )*ifelse(angle(z-1) <= 0, -one(z), one(z) )*fn^(-1/4)*airyaiprime(fn) ) ) end R1 = 0.0 R2 = 0.0 @@ -1092,7 +1092,7 @@ function asyAiry(np, y, alpha, T) R1 = R1 + (1/48*alpha^4 + 7/240*alpha^2 - 1/105)/np^1 + 1 R2 = R2 + (1/48*alpha^4 + 1/12*alpha^3 + 7/240*alpha^2 + 1/60*alpha + 1/140)/np^1 end - p = real( 4*sqrt(pi)/z^(1/4)/(d + 0im)^(1/4)*z^(-alpha/2)*( (R1*cos( (alpha + 1)/2*acos(2*z - 1) ) -cos( (alpha - 1)/2*acos(2*z - 1) )*R2)*fn^(1/4)*airy(0,fn) + 1im*(-sin( (alpha + 1)/2*acos(2*z - 1) )*R1 +sin( (alpha - 1)/2*acos(2*z - 1) )*R2)*ifelse(angle(z-1) <= 0, -one(z), one(z) )*fn^(-1/4)*airy(1,fn) ) ) + p = real( 4*sqrt(pi)/z^(1/4)/(d + 0im)^(1/4)*z^(-alpha/2)*( (R1*cos( (alpha + 1)/2*acos(2*z - 1) ) -cos( (alpha - 1)/2*acos(2*z - 1) )*R2)*fn^(1/4)*airyai(fn) + 1im*(-sin( (alpha + 1)/2*acos(2*z - 1) )*R1 +sin( (alpha - 1)/2*acos(2*z - 1) )*R2)*ifelse(angle(z-1) <= 0, -one(z), one(z) )*fn^(-1/4)*airyaiprime(fn) ) ) end @@ -1274,7 +1274,7 @@ end function asyAirygen(np, z, alpha, T::Int64, qm, m::Int64, UQ, fn, useQ::Bool, xin=NaN+NaN*1im) d = z - 1.0 +0im if T == 1 - return real( 4*sqrt(pi)/(z+0im)^(1/4+alpha/2)/d^(1/4)*(cos( (alpha + 1)/2*acos(2*z - 1+0im) )*fn^(1/4)*airy(0,fn) -1im*sin( (alpha + 1)/2*acos(2*z - 1+0im) )*ifelse(angle(z-1) <= 0, -one(z), one(z) )*fn^(-1/4)*airy(1,fn) ) ) + return real( 4*sqrt(pi)/(z+0im)^(1/4+alpha/2)/d^(1/4)*(cos( (alpha + 1)/2*acos(2*z - 1+0im) )*fn^(1/4)*airyai(fn) -1im*sin( (alpha + 1)/2*acos(2*z - 1+0im) )*ifelse(angle(z-1) <= 0, -one(z), one(z) )*fn^(-1/4)*airyaiprime(fn) ) ) end R = (1+0im)*[1 0] if useQ @@ -1300,7 +1300,7 @@ function asyAirygen(np, z, alpha, T::Int64, qm, m::Int64, UQ, fn, useQ::Bool, xi end end end - p = real( 4*sqrt(pi)/(z+0im)^(1/4+alpha/2)/d^(1/4)*( (R[1]*cos( (alpha + 1)/2*acos(2*z - 1+0im) ) -cos( (alpha - 1)/2*acos(2*z - 1+0im) )*R[2]*1im*4^alpha)*fn^(1/4)*airy(0,fn) + 1im*(-sin( (alpha + 1)/2*acos(2*z - 1+0im) )*R[1] +sin( (alpha - 1)/2*acos(2*z - 1+0im) )*R[2]*1im*4^alpha)*ifelse(angle(z-1) <= 0, -one(z), one(z) )*fn^(-1/4)*airy(1,fn) ) ) + p = real( 4*sqrt(pi)/(z+0im)^(1/4+alpha/2)/d^(1/4)*( (R[1]*cos( (alpha + 1)/2*acos(2*z - 1+0im) ) -cos( (alpha - 1)/2*acos(2*z - 1+0im) )*R[2]*1im*4^alpha)*fn^(1/4)*airyai(fn) + 1im*(-sin( (alpha + 1)/2*acos(2*z - 1+0im) )*R[1] +sin( (alpha - 1)/2*acos(2*z - 1+0im) )*R[2]*1im*4^alpha)*ifelse(angle(z-1) <= 0, -one(z), one(z) )*fn^(-1/4)*airyaiprime(fn) ) ) end # Additional short functions diff --git a/test/test_gausshermite.jl b/test/test_gausshermite.jl index 29258a5..d42de1e 100644 --- a/test/test_gausshermite.jl +++ b/test/test_gausshermite.jl @@ -2,14 +2,18 @@ tol = 1e-14 -n = 18; +n = 18; x,w = gausshermite( n ) +@test isa(x,Vector{Float64}) +@test isa(w,Vector{Float64}) @test (length(x) == n && length(w) == n) @test (dot(w,x) < tol && abs(dot(w,x.^2) - sqrt(pi)/2) < tol) # Test a small n: n = 42 x,w = gausshermite( n ) +@test isa(x,Vector{Float64}) +@test isa(w,Vector{Float64}) @test (length(x) == n && length(w) == n) @test (dot(w,x) < tol && abs(dot(w,x.^2) - sqrt(pi)/2) < tol) @test abs(x[37] - 5.660357581283058) < tol @@ -18,7 +22,9 @@ x,w = gausshermite( n ) # Test a larger n: n = 251 x,w = gausshermite( n ) +@test isa(x,Vector{Float64}) +@test isa(w,Vector{Float64}) @test (length(x) == n && length(w) == n) @test (dot(w,x) < tol && abs(dot(w,x.^2) - sqrt(pi)/2) < 300*tol) @test abs(x[37] - -13.292221459334638) < tol -@test abs(w[123] - 0.117419270715955) < 2*tol \ No newline at end of file +@test abs(w[123] - 0.117419270715955) < 2*tol diff --git a/test/test_gausslaguerre.jl b/test/test_gausslaguerre.jl index 181cc2c..024aece 100644 --- a/test/test_gausslaguerre.jl +++ b/test/test_gausslaguerre.jl @@ -14,6 +14,14 @@ tol = 3.e-13 tolGen = 2.e-9 # Possibly set srand(...) to get the same results for integration of a random polynomial. ns = [42; 251; 5000]; + +# the following were computed using quadgk +quadgk_exa = [3.2176248255905825e14 4.298744794628054e18 844.0328465154408 + 3.0483026539043077 902.5428056960185 1.6962837924111536 + 1.120730668666291 29.98453680681862 1.0146978179749486] + + +srand(0) for wei = 1:3 if wei == 1; alpha = 0.0;qm = 1.0; m = 1; elseif wei == 2; alpha = pi+0.1; qm = 1.0; m = 1; @@ -37,7 +45,9 @@ for wei = 1:3 end fcti(x) = fct(x).*exp(-qm*x.^m).*x^(alpha) # Test integration of a random polynomial which should be integrated exactly: - exa = quadgk(fcti, 0, Inf, reltol = tol)[1] + exa = quadgk_exa[ni,wei] + + for methi = 1:ifelse(ni == 1, 1, 5) if (ni == 1) # Golub-Welsch is the accurate algorithm for low n but does not allow m or qm != 1. @@ -137,7 +147,7 @@ for alpha = [0.0; 4.15] @test abs(-2395.51952258921326919097391744358588135/exp(0.01*4*np/2)/factorp -aRHe)/abs(aRHe) < tolEx end - for n = ceil(Int64, linspace(111, 5000, 5) ) + for n = ceil.([Int64], linspace(111, 5000, 5) ) for z = linspace(1e-4, 0.99, 10) aRH = polyAsyRH(n, 4*n*z, alpha, T) @test abs(polyAsyRHgen(n, 4*n*z, alpha, T, 1.0, 1, UQ) -aRH)/abs(aRH) < 2*tolEx diff --git a/test/test_gausslegendre.jl b/test/test_gausslegendre.jl index fc2d77d..c61bbf2 100644 --- a/test/test_gausslegendre.jl +++ b/test/test_gausslegendre.jl @@ -1,3 +1,5 @@ + + # Test for gausslegendre(). @@ -10,8 +12,8 @@ tol = 1e-14 n = 42 x, w = gausslegendre(n) @test length(x) == n && length(w) == n -@test x[37] ≈ 0.910959724904127 atol=tol -@test w[37] ≈ 0.030479240699603 atol=tol +@test ≈(x[37],0.910959724904127;atol=tol) +@test ≈(w[37],0.030479240699603;atol=tol) @test dot( w,(x.^2)) ≈ 2/3 @test dot( w,exp.(x)) ≈ exp(1)-exp(-1) From 51834b5d456ece953f60e8a4c27bab799dc8e1e0 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Wed, 15 Feb 2017 11:06:46 +1100 Subject: [PATCH 4/4] update for latest 0.6 --- REQUIRE | 1 + src/FastGaussQuadrature.jl | 4 +++- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/REQUIRE b/REQUIRE index 84feefe..2a2c164 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,2 +1,3 @@ julia 0.5 Compat 0.8.8 +SpecialFunctions 0.1 diff --git a/src/FastGaussQuadrature.jl b/src/FastGaussQuadrature.jl index ad251f2..52b17df 100644 --- a/src/FastGaussQuadrature.jl +++ b/src/FastGaussQuadrature.jl @@ -1,6 +1,6 @@ module FastGaussQuadrature -using Compat +using Compat, SpecialFunctions export gausslegendre export gausschebyshev @@ -11,6 +11,8 @@ export gausslobatto export gaussradau export besselroots +import SpecialFunctions: besselj, airyai, airyaiprime + include("gausslegendre.jl") include("gausschebyshev.jl") include("gausslaguerre.jl")