Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Faster gaussjacobi and gausslobatto, and refactoring #85

Merged
merged 8 commits into from
Jan 4, 2021
Merged

Faster gaussjacobi and gausslobatto, and refactoring #85

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
02c25d8
fix typos and error handling
hyrodium Dec 31, 2020
6c15641
update gitignore
hyrodium Dec 31, 2020
177d77a
refactoring
hyrodium Dec 31, 2020
9c32545
refactoring, and making type stable of gaussjacobi
hyrodium Dec 31, 2020
f566b88
update comments, minor fixes
hyrodium Dec 31, 2020
5a806e9
minor fixes, update error handling
hyrodium Dec 31, 2020
78eea8e
add more tests
hyrodium Dec 31, 2020
1a975de
fix runtests.jl to check all test even if some tests fails
hyrodium Jan 3, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion .gitignore
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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