Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
41 changes: 23 additions & 18 deletions src/char-values.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,8 @@ export charλ,
charA,
charB
"""
charλ(q,ν,k)
charλ(q, ν, k::AbstractRange)
charλ(q, ν, k::Integer)

char value λ_(ν+k) for Mathieu's equation

Expand All @@ -13,11 +14,10 @@ where
q ∈ ℝ - parameter
ν ∈ [-1,1] - fractional part of the non-integer order
k ∈ ℤ⁺ - range of integer parts of the order

"""
function charλ(nu_::Real, q::Real; k::AbstractRange=1:1) # reduced = true
#nu = reduced ? rem(nu_+1,2)-1 : nu_;
nu = rem(nu_+1,2)-1;
function charλ(nu_::Real, q::Real, k::AbstractRange) # reduced = true
#nu = reduced ? rem(nu_+1,2)-1 : nu_
nu = rem(nu_+1,2)-1

# Set matrix size using formula from Shirts paper (1993), Eqs. (2.1)-(2.2).
nu0 = nu + maximum(k)
Expand All @@ -33,8 +33,11 @@ function charλ(nu_::Real, q::Real; k::AbstractRange=1:1) # reduced = true
return a
end

charλ(nu_::Real, q::Real, k::Integer) = charλ(nu_, q, k:k)[1]

"""
charA(q; k=0:4)
charA(q, k::AbstractRange)
charA(q, k::Integer)

char value A_k for Mathieu's equation

Expand All @@ -44,9 +47,8 @@ where

q ∈ ℝ - parameter
k ∈ ℤ⁺ - eigenvalue index

"""
function charA(q::Real; k::AbstractRange=0:0)
function charA(q::Real, k::AbstractRange)
all(x -> x >= 0, k) || throw(DomainError(k, "Indices must be non-negative integers."))

# Boolean indices of even and odd n values
Expand All @@ -55,22 +57,24 @@ function charA(q::Real; k::AbstractRange=0:0)

a = Array{Float64}(undef ,length(k))
k1 = k .+ 1
a[ie] = charλ(0.0, abs(q); k = k1)[ie]
a[ie] = charλ(0.0, abs(q), k1)[ie]
if q>=0
a[io] = charλ(one(q), q; k = k1)[io]
a[io] = charλ(one(q), q, k1)[io]
else
if 0 in k # maybe not the cleanest way to do it
a[io] = charλ(one(q), abs(q); k = k[2]:last(k))[io[2:end]]
a[io] = charλ(one(q), abs(q), k[2]:last(k))[io[2:end]]
else
a[io] = charλ(one(q), abs(q); k=k)[io]
a[io] = charλ(one(q), abs(q), k)[io]
end
end
return a
end

charA(q::Real, k::Integer) = charA(q, k:k)[1]

"""
charB(q,k)
charB(q, k::AbstractRange)
charB(q, k::Integer)

char value B_k for Mathieu's equation

Expand All @@ -80,20 +84,21 @@ where

q ∈ ℝ - parameter
k ∈ ℤ - eigenvalueindex

"""
function charB(q::Real; k::AbstractRange=1:1)
function charB(q::Real, k::AbstractRange)
all(x -> x > 0, k) || throw(DomainError(k, "Indices must be positive integers."))
# Boolean indices of even and odd n values
ie = map(iseven, k)
io = map(!, ie)

b = Array{Float64}(undef, length(k))
b[ie] = charλ(0.0,q,k=k)[ie]
b[ie] = charλ(0.0, q, k)[ie]
if q>=0
b[io] = charλ(1.0,q,k=k)[io]
b[io] = charλ(1.0, q, k)[io]
else
b[io] = charλ(1.0,abs(q),k = (k .+ 1))[io]
b[io] = charλ(1.0,abs(q), (k .+ 1))[io]
end
return b
end

charB(q::Real, k::Integer) = charB(q, k:k)[1]
26 changes: 17 additions & 9 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -11,52 +11,60 @@ function tapprox(a, b; atol=1e-15)
end

@testset "basic" begin
@test charλ(1, 0, k=1:10) == [1.0, 1.0, 9.0, 9.0, 25.0, 25.0, 49.0, 49.0, 81.0, 81.0]
@test maximum(charA(0,k=0:100) - [0:100;].^2) == 0
@test norm(charB(0,k=1:100) - [1:100;].^2) == 0
@test charλ(1, 0, 1:10) == [1.0, 1.0, 9.0, 9.0, 25.0, 25.0, 49.0, 49.0, 81.0, 81.0]
@test maximum(charA(0, 0:100) - [0:100;].^2) == 0
@test norm(charB(0, 1:100) - [1:100;].^2) == 0
end

@testset "interface" begin
q = 1.0
k = 2
nu = .5
@test charλ(nu, q, k) == charλ(nu, q, k:k)[1]
@test charA(q, k) == charA(q, k:k)[1]
@test charB(q, k) == charB(q, k:k)[1]
end

filename = "MathieuCharacteristicA-1.csv"
@testset "$filename" begin
test1 = readcsv(filename)
r = reduce(hcat,Vector{Float64}[charA(q,k=0:10) for q in [-10:.01:10;]])
r = reduce(hcat,Vector{Float64}[charA(q, 0:10) for q in [-10:.01:10;]])
@test tapprox(test1, r; atol=7.5e-13)
end

filename = "MathieuCharacteristicA-2.csv"
@testset "$filename" begin
test1 = readcsv(filename)
r = reduce(hcat,Vector{Float64}[charA(q,k=0:3) for q in [30:.01:50;]])
r = reduce(hcat,Vector{Float64}[charA(q, 0:3) for q in [30:.01:50;]])
@test tapprox(test1, r; atol=7.6e-13) # NOTE: was 7.5e-13
end

filename = "MathieuCharacteristicB-1.csv"
@testset "$filename" begin
test1 = readcsv(filename)
r = reduce(hcat,Vector{Float64}[charB(q,k=1:10) for q in [-10:.01:10;]])
r = reduce(hcat,Vector{Float64}[charB(q, 1:10) for q in [-10:.01:10;]])
@test tapprox(test1, r; atol=7.5e-13)
end

filename = "MathieuCharacteristicB-2.csv"
@testset "$filename" begin
test1 = readcsv(filename)
r = reduce(hcat,Vector{Float64}[charB(q,k=1:3) for q in [30:.01:50;]])
r = reduce(hcat,Vector{Float64}[charB(q, 1:3) for q in [30:.01:50;]])
@test tapprox(test1, r; atol=2.8e-11)
end

filename = "MathieuCharacteristicL-1.csv"
@testset "$filename" begin
test1 = readcsv(filename)[1:100,:]
test2 = Float64[charλ(q,ν,k=1:1)[1] for ν in [0:.01:0.99;], q in [-5:.01:5;]]
test2 = Float64[charλ(q,ν, 1:1)[1] for ν in [0:.01:0.99;], q in [-5:.01:5;]]
@test_broken tapprox(test1, test2, atol=7.5e-15)
# TODO: test ν > 1 (currently failing)
end

filename = "MathieuCharacteristicL-2.csv"
@testset "$filename" begin
test1 = readcsv(filename)[1:100,:]
test2 = Float64[charλ(q,ν,k=1:1)[1] for ν in [0:.01:0.99;], q in [30:.01:50;]]
test2 = Float64[charλ(q,ν, 1:1)[1] for ν in [0:.01:0.99;], q in [30:.01:50;]]
@test_broken tapprox(test1, test2, atol=4.5e-14)
# TODO: test ν > 1 (currently failing)
end