Skip to content

Commit

Permalink
Deprecate gausschebyshev and add gausschebyshevt, ..., `gausscheb…
Browse files Browse the repository at this point in the history
…yshevw` (#123)

* deprecate gausschebyshev and add gausschebyshev1 to gausschebyshev4

* fix docs

* fix error handling for gausschebyshev

* update docstrings

* rename gausschebyshev* function

* rename `gausschebyshevT` to `gausschebyshev` etc.

* bump version to v0.5.2
  • Loading branch information
hyrodium authored Oct 2, 2023
1 parent 7023121 commit de3eb8f
Show file tree
Hide file tree
Showing 5 changed files with 84 additions and 44 deletions.
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
name = "FastGaussQuadrature"
uuid = "442a2c76-b920-505d-bb47-c5924d526838"
version = "0.5.1"
version = "0.5.2"

[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Expand Down
5 changes: 4 additions & 1 deletion docs/src/gaussquadrature.md
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,10 @@ There are four kinds of Gauss-Chebyshev quadrature rules, corresponding to four
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)
gausschebyshevt(n::Integer)
gausschebyshevu(n::Integer)
gausschebyshevv(n::Integer)
gausschebyshevw(n::Integer)
```


Expand Down
4 changes: 4 additions & 0 deletions src/FastGaussQuadrature.jl
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,10 @@ using StaticArrays

export gausslegendre
export gausschebyshev
export gausschebyshevt
export gausschebyshevu
export gausschebyshevv
export gausschebyshevw
export gausslaguerre
export gausshermite
export gaussjacobi
Expand Down
77 changes: 51 additions & 26 deletions src/gausschebyshev.jl
Original file line number Diff line number Diff line change
@@ -1,6 +1,5 @@
@doc raw"""
gausschebyshev(n::Integer) -> x, w # nodes, weights
gausschebyshev(n::Integer, 1) -> x, w # nodes, weights
gausschebyshevt(n::Integer) -> x, w # nodes, weights
Return nodes `x` and weights `w` of [Gauss-Chebyshev quadrature](https://en.wikipedia.org/wiki/Chebyshev%E2%80%93Gauss_quadrature) of the 1st kind.
Expand All @@ -10,7 +9,7 @@ Return nodes `x` and weights `w` of [Gauss-Chebyshev quadrature](https://en.wiki
# Examples
```jldoctest
julia> x, w = gausschebyshev(3);
julia> x, w = gausschebyshevt(3);
julia> f(x) = x^4;
Expand All @@ -19,10 +18,16 @@ julia> I = dot(w, f.(x));
julia> I ≈ 3π/8
true
```
"""
function gausschebyshevt(n::Integer)
if n < 0
throw(DomainError(n, "Input n must be a non-negative integer"))
end
return [cos((2 * k - 1) * π / (2 * n)) for k = n:-1:1], fill/ n, n)
end

---
gausschebyshev(n::Integer, 2) -> x, w # nodes, weights
@doc raw"""
gausschebyshevu(n::Integer) -> x, w # nodes, weights
Return nodes `x` and weights `w` of [Gauss-Chebyshev quadrature](https://en.wikipedia.org/wiki/Chebyshev%E2%80%93Gauss_quadrature) of the 2nd kind.
Expand All @@ -32,7 +37,7 @@ Return nodes `x` and weights `w` of [Gauss-Chebyshev quadrature](https://en.wiki
# Examples
```jldoctest
julia> x, w = gausschebyshev(3, 2);
julia> x, w = gausschebyshevu(3);
julia> f(x) = x^4;
Expand All @@ -41,10 +46,16 @@ julia> I = dot(w, f.(x));
julia> I ≈ π/16
true
```
"""
function gausschebyshevu(n::Integer)
if n < 0
throw(DomainError(n, "Input n must be a non-negative integer"))
end
return [cos(k * π / (n + 1)) for k = n:-1:1], [π/(n + 1) * sin(k / (n + 1) * π)^2 for k = n:-1:1]
end

---
gausschebyshev(n::Integer, 3) -> x, w # nodes, weights
@doc raw"""
gausschebyshevv(n::Integer) -> x, w # nodes, weights
Return nodes `x` and weights `w` of [Gauss-Chebyshev quadrature](https://en.wikipedia.org/wiki/Chebyshev%E2%80%93Gauss_quadrature) of the 3rd kind.
Expand All @@ -54,7 +65,7 @@ Return nodes `x` and weights `w` of [Gauss-Chebyshev quadrature](https://en.wiki
# Examples
```jldoctest
julia> x, w = gausschebyshev(3, 3);
julia> x, w = gausschebyshevv(3);
julia> f(x) = x^4;
Expand All @@ -63,10 +74,16 @@ julia> I = dot(w, f.(x));
julia> I ≈ 3π/8
true
```
"""
function gausschebyshevv(n::Integer)
if n < 0
throw(DomainError(n, "Input n must be a non-negative integer"))
end
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]
end

---
gausschebyshev(n::Integer, 4) -> x, w # nodes, weights
@doc raw"""
gausschebyshevw(n::Integer) -> x, w # nodes, weights
Return nodes `x` and weights `w` of [Gauss-Chebyshev quadrature](https://en.wikipedia.org/wiki/Chebyshev%E2%80%93Gauss_quadrature) of the 4th kind.
Expand All @@ -76,7 +93,7 @@ Return nodes `x` and weights `w` of [Gauss-Chebyshev quadrature](https://en.wiki
# Examples
```jldoctest
julia> x, w = gausschebyshev(3, 4);
julia> x, w = gausschebyshevw(3);
julia> f(x) = x^4;
Expand All @@ -86,6 +103,13 @@ julia> I ≈ 3π/8
true
```
"""
function gausschebyshevw(n::Integer)
if n < 0
throw(DomainError(n, "Input n must be a non-negative integer"))
end
return [cos(k * π / (n + .5)) for k = n:-1:1], [2π / (n + .5) * sin(k * π / (2 * (n + .5)))^2 for k = n:-1:1]
end

function gausschebyshev(n::Integer, kind::Integer=1)
# GAUSS-CHEBYSHEV NODES AND WEIGHTS.

Expand All @@ -95,20 +119,21 @@ function gausschebyshev(n::Integer, kind::Integer=1)

# Use known explicit formulas. Complexity O(n).
if kind == 1
# Gauss-ChebyshevT quadrature, i.e., w(x) = 1/sqrt(1-x^2)
return ([cos((2 * k - 1) * π / (2 * n)) for k = n:-1:1], fill/ n, n))
# Gauss-Chebyshevt quadrature, i.e., w(x) = 1/sqrt(1-x^2)
Base.depwarn("`gausschebyshev(n, 1)` is deprecated and will be removed in the next breaking release. Please use `gausschebyshevt(n)` instead.", :gausschebyshev)
return gausschebyshevt(n)
elseif kind == 2
# Gauss-ChebyshevU quadrature, i.e., w(x) = sqrt(1-x^2)
return ([cos(k * π / (n + 1)) for k = n:-1:1],
/(n + 1) * sin(k / (n + 1) * π)^2 for k = n:-1:1])
# Gauss-Chebyshevu quadrature, i.e., w(x) = sqrt(1-x^2)
Base.depwarn("`gausschebyshev(n, 2)` is deprecated and will be removed in the next breaking release. Please use `gausschebyshevu(n)` instead.", :gausschebyshev)
return gausschebyshevu(n)
elseif kind == 3
# Gauss-ChebyshevV quadrature, i.e., w(x) = sqrt((1+x)/(1-x))
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])
# Gauss-Chebyshevv quadrature, i.e., w(x) = sqrt((1+x)/(1-x))
Base.depwarn("`gausschebyshev(n, 3)` is deprecated and will be removed in the next breaking release. Please use `gausschebyshevv(n)` instead.", :gausschebyshev)
return gausschebyshevv(n)
elseif kind == 4
# Gauss-ChebyshevW quadrature, i.e., w(x) = sqrt((1-x)/(1+x))
return ([cos(k * π / (n + .5)) for k = n:-1:1],
[2π / (n + .5) * sin(k * π / (2 * (n + .5)))^2 for k = n:-1:1])
# Gauss-Chebyshevw quadrature, i.e., w(x) = sqrt((1-x)/(1+x))
Base.depwarn("`gausschebyshev(n, 4)` is deprecated and will be removed in the next breaking release. Please use `gausschebyshevw(n)` instead.", :gausschebyshev)
return gausschebyshevw(n)
else
throw(ArgumentError("Chebyshev kind should be 1, 2, 3, or 4"))
end
Expand Down
40 changes: 24 additions & 16 deletions test/test_gausschebyshev.jl
Original file line number Diff line number Diff line change
@@ -1,34 +1,42 @@
@testset "Gauss–Chebyshev" begin

@testset "Check error" begin
@test_throws DomainError gausschebyshev(-1,1)
@test_throws DomainError gausschebyshev(-1,2)
@test_throws DomainError gausschebyshev(-1,3)
@test_throws DomainError gausschebyshev(-1,4)
@test_throws ArgumentError gausschebyshev(0,5)
@test_throws DomainError gausschebyshevt(-1)
@test_throws DomainError gausschebyshevu(-1)
@test_throws DomainError gausschebyshevv(-1)
@test_throws DomainError gausschebyshevw(-1)
end

n = 10

@testset "x.^2" begin
x, w = gausschebyshev(n,1)
x, w = gausschebyshevt(n)
@test dot(x.^2,w) π/2
x, w = gausschebyshev(n,2)
x, w = gausschebyshevu(n)
@test dot(x.^2,w) π/8
x, w = gausschebyshev(n,3)
x, w = gausschebyshevv(n)
@test dot(x.^2,w) π/2
x, w = gausschebyshev(n,4)
x, w = gausschebyshevw(n)
@test dot(x.^2,w) π/2
end

@testset "x^3" begin
x, w = gausschebyshev(n,1)
@test abs(dot(x.^3,w))<1e-15
x, w = gausschebyshev(n,2)
@test abs(dot(x.^3,w))<1e-15
x, w = gausschebyshev(n,3)
x, w = gausschebyshevt(n)
@test abs(dot(x.^3,w)) < 1e-15
x, w = gausschebyshevu(n)
@test abs(dot(x.^3,w)) < 1e-15
x, w = gausschebyshevv(n)
@test dot(x.^3,w) 3*π/8
x, w = gausschebyshev(n,4)
x, w = gausschebyshevw(n)
@test dot(x.^3,w) -3*π/8
end
end

@testset "deprecated" begin
n = 42
@test gausschebyshevt(n) == gausschebyshev(n, 1)
@test gausschebyshevu(n) == gausschebyshev(n, 2)
@test gausschebyshevv(n) == gausschebyshev(n, 3)
@test gausschebyshevw(n) == gausschebyshev(n, 4)
@test_throws ArgumentError gausschebyshev(0,5)
end
end

0 comments on commit de3eb8f

Please sign in to comment.