Skip to content

Commit

Permalink
Faster gaussjacobi and gausslobatto, and refactoring (#85)
Browse files Browse the repository at this point in the history 
* fix typos and error handling

* update gitignore

* refactoring

* refactoring, and making type stable of gaussjacobi

* update comments, minor fixes

* minor fixes, update error handling

* add more tests

* fix runtests.jl to check all test even if some tests fails
  • Loading branch information
@hyrodium
hyrodium authored Jan 4, 2021
1 parent ba2da14 commit c52b3bc
Show file tree
Hide file tree
Showing 19 changed files with 317 additions and 241 deletions.
2 changes: 1 addition & 1 deletion .gitignore
View file
Original file line number Diff line number Diff line change
@@ -1,2 +1,2 @@

.DS_Store
/Manifest.toml
1 change: 0 additions & 1 deletion src/FastGaussQuadrature.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,6 @@ sum(A) = Base.sum(A)
flipdim(A, d) = reverse(A, dims=d)



export gausslegendre
export gausschebyshev
export gausslaguerre
Expand Down
41 changes: 22 additions & 19 deletions src/besselroots.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ function besselroots(nu::Float64, n::Integer)
#BESSELROOTS The first N roots of the function J_v(x)

# DEVELOPERS NOTES:
# V = 0 --> Full double precision for N <= 20 (Wolfram Alpha), and very
# V = 0 --> Full Float64 precision for N <= 20 (Wolfram Alpha), and very
# accurate approximations for N > 20 (McMahon's expansion)
# -1 <= V <= 5 : V ~= 0 --> 12 decimal figures for the 6 first zeros
# (Piessens's Chebyshev series approximations), and very accurate
Expand All @@ -14,7 +14,7 @@ function besselroots(nu::Float64, n::Integer)
# Later modified by A. Townsend to work in Julia

if n < 0
error("Input N must be a positive integer")
throw(DomainError(n, "Input N must be a non-negative integer"))
end

x = zeros(n)
Expand All @@ -25,7 +25,7 @@ function besselroots(nu::Float64, n::Integer)
for k in min(n,20)+1:n
x[k] = McMahon(nu, k)
end
elseif n>0 && nu >= -1 && nu <= 5
elseif n > 0 && nu >= -1 && nu <= 5
correctFirstFew = Piessens( nu )
for k in 1:min(n,6)
x[k] = correctFirstFew[k]
Expand All @@ -38,7 +38,7 @@ function besselroots(nu::Float64, n::Integer)
x[k] = McMahon(nu, k)
end
end
x
return x
end

function McMahon(nu::Float64, k::Integer)
Expand All @@ -58,8 +58,11 @@ function McMahon(nu::Float64, k::Integer)
b = 0.25 * (2nu+4k-1)*pi
# 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)
return x
end

# Roots of Bessel funcion ``J_0`` in Float64.
# https://mathworld.wolfram.com/BesselFunctionZeros.html
const J0_roots =
[ 2.4048255576957728
5.5200781102863106
Expand Down Expand Up @@ -100,23 +103,23 @@ const Piessens_C = [
0.000001164419 0.000001903082 0.000000051084 0.000000003677 0.000000000448 0.000000000077
-0.000000485189 -0.000000781030 -0.000000015501 -0.000000000896 -0.000000000092 -0.000000000014
0.000000203309 0.000000322648 0.000000004736 0.000000000220 0.000000000019 0.000000000002
-0.000000085602 -0.000000134047 -0.000000001456 -0.000000000054 -0.000000000004 0.
0.000000036192 0.000000055969 0.000000000450 0.000000000013 0. 0.
-0.000000015357 -0.000000023472 -0.000000000140 -0.000000000003 0. 0.
0.000000006537 0.000000009882 0.000000000043 0.000000000001 0. 0.
-0.000000002791 -0.000000004175 -0.000000000014 0. 0. 0.
0.000000001194 0.000000001770 0.000000000004 0. 0. 0.
-0.000000000512 -0.000000000752 0. 0. 0. 0.
0.000000000220 0.000000000321 0. 0. 0. 0.
-0.000000000095 -0.000000000137 0. 0. 0. 0.
0.000000000041 0.000000000059 0. 0. 0. 0.
-0.000000000018 -0.000000000025 0. 0. 0. 0.
0.000000000008 0.000000000011 0. 0. 0. 0.
-0.000000000003 -0.000000000005 0. 0. 0. 0.
0.000000000001 0.000000000002 0. 0. 0. 0.]
-0.000000085602 -0.000000134047 -0.000000001456 -0.000000000054 -0.000000000004 0
0.000000036192 0.000000055969 0.000000000450 0.000000000013 0 0
-0.000000015357 -0.000000023472 -0.000000000140 -0.000000000003 0 0
0.000000006537 0.000000009882 0.000000000043 0.000000000001 0 0
-0.000000002791 -0.000000004175 -0.000000000014 0 0 0
0.000000001194 0.000000001770 0.000000000004 0 0 0
-0.000000000512 -0.000000000752 0 0 0 0
0.000000000220 0.000000000321 0 0 0 0
-0.000000000095 -0.000000000137 0 0 0 0
0.000000000041 0.000000000059 0 0 0 0
-0.000000000018 -0.000000000025 0 0 0 0
0.000000000008 0.000000000011 0 0 0 0
-0.000000000003 -0.000000000005 0 0 0 0
0.000000000001 0.000000000002 0 0 0 0]


function Piessens( nu::Float64 )
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:
C = Piessens_C
Expand Down
16 changes: 10 additions & 6 deletions src/gausschebyshev.jl
View file
Original file line number Diff line number Diff line change
@@ -1,22 +1,26 @@
function gausschebyshev(n::Integer, kind::Integer=1)
# GAUSS-CHEBYSHEV NODES AND WEIGHTS.

if n < 0
throw(DomainError(n, "Input n must be a non-negative integer"))
end

# Use known explicit formulas. Complexity O(n).
if kind == 1
# Gauss-ChebyshevT quadrature, i.e., w(x) = 1/sqrt(1-x^2)
[cos((2 * k - 1) * π / (2 * n)) for k = n:-1:1], fill(π / n, n)
return ([cos((2 * k - 1) * π / (2 * n)) for k = n:-1:1], fill(π / n, n))
elseif kind == 2
# Gauss-ChebyshevU quadrature, i.e., w(x) = sqrt(1-x^2)
([cos(k * π / (n + 1)) for k = n:-1:1],
[π/(n + 1) * sin(k / (n + 1) * π)^2 for k = n:-1:1])
return ([cos(k * π / (n + 1)) for k = n:-1:1],
[π/(n + 1) * sin(k / (n + 1) * π)^2 for k = n:-1:1])
elseif kind == 3
# Gauss-ChebyshevV quadrature, i.e., w(x) = sqrt((1+x)/(1-x))
([cos((k - .5) * π / (n + .5)) for k = n:-1:1],
return ([cos((k - .5) * π / (n + .5)) for k = n:-1:1],
[2π / (n + .5) * cos((k - .5) * π / (2 * (n + .5)))^2 for k = n:-1:1])
elseif kind == 4
# Gauss-ChebyshevW quadrature, i.e., w(x) = sqrt((1-x)/(1+x))
([cos(k * π / (n + .5)) for k = n:-1:1],
[2π / (n + .5) * sin(k * π / (2 * (n + .5)))^2 for k = n:-1:1])
return ([cos(k * π / (n + .5)) for k = n:-1:1],
[2π / (n + .5) * sin(k * π / (2 * (n + .5)))^2 for k = n:-1:1])
else
throw(ArgumentError("Chebyshev kind should be 1, 2, 3, or 4"))
end
Expand Down
72 changes: 36 additions & 36 deletions src/gausshermite.jl
View file
Original file line number Diff line number Diff line change
@@ -1,45 +1,44 @@
function gausshermite( n::Integer )
function gausshermite(n::Integer)
x,w = unweightedgausshermite(n)
w .*= exp.(-x.^2)
x, w
end
function unweightedgausshermite( n::Integer )
function unweightedgausshermite(n::Integer)
# GAUSSHERMITE(n) COMPUTE THE GAUSS-HERMITE NODES AND WEIGHTS IN O(n) time.
if n < 0
x = (Float64[],Float64[])
return x
throw(DomainError(n, "Input n must be a non-negative integer"))
elseif n == 0
x = (Float64[],Float64[])
return x
return Float64[],Float64[]
elseif n == 1
x = ([0.0],[sqrt(pi)])
return x
return [0.0],[sqrt(pi)]
elseif n <= 20
# GW algorithm
x = hermpts_gw( n )
x = hermpts_gw(n)
elseif n <= 200
# REC algorithm
x = hermpts_rec( n )
x = hermpts_rec(n)
else
# ASY algorithm
x = hermpts_asy( n )
x = hermpts_asy(n)
end

if mod(n,2) == 1 # fold out
# fold out
if mod(n,2) == 1
w = [flipdim(x[2][:],1); x[2][2:end]]
x = [-flipdim(x[1],1) ; x[1][2:end]]
else
w = [flipdim(x[2][:],1); x[2][:]]
x = [-flipdim(x[1],1) ; x[1]]
end
w .*= sqrt(π)/sum(exp.(-x.^2).*w)
(x, w)

return x, w
end

function hermpts_asy( n::Integer )
function hermpts_asy(n::Integer)
# Compute Hermite nodes and weights using asymptotic formula

x0 = HermiteInitialGuesses( n ) # get initial guesses
x0 = HermiteInitialGuesses(n) # get initial guesses
t0 = x0./sqrt(2n+1)
theta0 = acos.(t0) # convert to theta-variable
val = x0;
Expand All @@ -56,13 +55,13 @@ function hermpts_asy( n::Integer )
w = x.*val[1] .+ sqrt(2).*val[2]
w .= 1 ./ w.^2; # quadrature weights

(x, w)
return x, w
end

function hermpts_rec( n::Integer )
function hermpts_rec(n::Integer)
# Compute Hermite nodes and weights using recurrence relation.

x0 = HermiteInitialGuesses( n )
x0 = HermiteInitialGuesses(n)
x0 .*= sqrt(2)
val = x0
for _ = 1:10
Expand All @@ -77,12 +76,12 @@ function hermpts_rec( n::Integer )
x0 ./= sqrt(2)
w = 1 ./ last.(val).^2 # quadrature weights

x = (x0, w)
return x0, w
end

function hermpoly_rec( n::Integer, x0)
function hermpoly_rec(n::Integer, x0)
# HERMPOLY_rec evaluation of scaled Hermite poly using recurrence
n < 0 && throw(ArgumentError("n = $n must be positive"))
n < 0 && throw(DomainError(n, "Input n must be a non-negative integer"))
# evaluate:
w = exp(-x0^2 / (4*n))
wc = 0 # 0 times we've applied wc
Expand All @@ -103,15 +102,14 @@ function hermpoly_rec( n::Integer, x0)
Hold *= w
end

# return (value, derivative):
val = (H, -x0*H + sqrt(n)*Hold)
return H, -x0*H + sqrt(n)*Hold
end

function hermpoly_rec( r::Base.OneTo, x0)
function hermpoly_rec(r::Base.OneTo, x0)
isempty(r) && return [1.0]
n = maximum(r)
# HERMPOLY_rec evaluation of scaled Hermite poly using recurrence
n < 0 && throw(ArgumentError("n = $n must be positive"))
n < 0 && throw(DomainError(n, "Input n must be a non-negative integer"))
n == 0 && return [exp(-x0^2 / 4)]
p = max(1,floor(Int,x0^2/100))
w = exp(-x0^2 / (4*p))
Expand All @@ -134,10 +132,10 @@ function hermpoly_rec( r::Base.OneTo, x0)
end
ret .*= w^(p-wc)

ret
return ret
end

hermpoly_rec( r::AbstractRange, x0) = hermpoly_rec(Base.OneTo(maximum(r)), x0)[r.+1]
hermpoly_rec(r::AbstractRange, x0) = hermpoly_rec(Base.OneTo(maximum(r)), x0)[r.+1]


function hermpoly_asy_airy(n::Integer, theta)
Expand Down Expand Up @@ -167,11 +165,11 @@ function hermpoly_asy_airy(n::Integer, theta)
u2 = @. (-9*cosT^4 + 249*cosT^2 + 145)/1152
u3 = @. (-4042*cosT^9+18189*cosT^7-28287*cosT^5-151995*cosT^3-259290*cosT)/414720

#first term
# first term
A0 = 1
val = A0*Airy0

#second term
# second term
B0 = @. -(a0*phi^6 * u1+a1*u0)/chi^2
val .+= B0.*Airy1./musq.^(4/3)

Expand Down Expand Up @@ -210,17 +208,17 @@ function hermpoly_asy_airy(n::Integer, theta)
C1 = @. -(phi^18 * v3 + b1*phi^12 * v2 + b2*phi^6 * v1 + b3*v0)/chi^4
dval = dval + C1.*Airy0/musq.^(2/3+2)

#fourth term
# fourth term
D1 = @. (a0*phi^12 * v2 + a1*phi^6 * v1 + a2*v0)/chi^3
dval = dval + D1.*Airy1/musq.^2

dval = C.*dval

val = (val, dval)
return val, dval
end

let T(t) = @. t^(2/3)*(1+5/48*t^(-2)-5/36*t^(-4)+(77125/82944)*t^(-6) -108056875/6967296*t^(-8)+162375596875/334430208*t^(-10))
global function HermiteInitialGuesses( n::Integer )
global function HermiteInitialGuesses(n::Integer)
#HERMITEINTITIALGUESSES(N), Initial guesses for Hermite zeros.
#
# [1] L. Gatteschi, Asymptotics and bounds for the zeros of Laguerre
Expand All @@ -247,7 +245,9 @@ let T(t) = @. t^(2/3)*(1+5/48*t^(-2)-5/36*t^(-4)+(77125/82944)*t^(-6) -108056875

airyrts = -T(3/8*pi*(4*(1:m) .- 1))

airyrts_exact = [-2.338107410459762 # Exact Airy roots.
# Roots of Airy function in Float64
# https://mathworld.wolfram.com/AiryFunctionZeros.html
airyrts_exact = [-2.338107410459762
-4.087949444130970
-5.520559828095555
-6.786708090071765
Expand All @@ -262,7 +262,7 @@ let T(t) = @. t^(2/3)*(1+5/48*t^(-2)-5/36*t^(-4)+(77125/82944)*t^(-6) -108056875
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) )
x_init_airy = flipdim(x_init,1)

# Tricomi initial guesses. Equation (2.1) in [1]. Originally in [2].
# These initial guesses are good near x = 0 . Note: zeros of besselj(+/-.5,x)
Expand Down Expand Up @@ -299,10 +299,10 @@ let T(t) = @. t^(2/3)*(1+5/48*t^(-2)-5/36*t^(-4)+(77125/82944)*t^(-6) -108056875
end


function hermpts_gw( n::Integer )
function hermpts_gw(n::Integer)
# Golub--Welsch algorithm. Used here for n<=20.

beta = sqrt.(0.5 .* (1:n-1)) # 3-term recurrence coeffs
beta = sqrt.((1:n-1)/2) # 3-term recurrence coeffs
T = SymTridiagonal(zeros(n), beta) # Jacobi matrix
(D, V) = eigen(T) # Eigenvalue decomposition
indx = sortperm(D) # Hermite points
Expand Down
Loading

0 comments on commit c52b3bc

Please sign in to comment.