diff --git a/test/SecondPoincareInvariants/test_SecondPoincareInvariants.jl b/test/SecondPoincareInvariants/test_SecondPoincareInvariants.jl index 2f8c7fe..1e6ca27 100644 --- a/test/SecondPoincareInvariants/test_SecondPoincareInvariants.jl +++ b/test/SecondPoincareInvariants/test_SecondPoincareInvariants.jl @@ -30,8 +30,8 @@ end @safetestset "Consistency Between Implementations" begin using PoincareInvariants - using PoincareInvariants.SecondPoincareInvariants.Chebyshev: ChebyshevPlan - using PoincareInvariants.SecondPoincareInvariants.FiniteDifferences: FiniteDiffPlan + using PoincareInvariants.SecondPoincareInvariants: ChebyshevPlan + using PoincareInvariants.SecondPoincareInvariants: FiniteDiffPlan D = 6 Ω(z, t, p) = [ @@ -52,15 +52,15 @@ end x + y ) - Idefault = let pinv = PI2{Float64}(Ω, D, 10_000) + Idefault = let pinv = PI2{Float64}(Ω, D, 1_000) compute!(pinv, getpoints(f, pinv), 0, nothing) end - Icheb = let pinv = PI2{Float64}(Ω, D, 10_000, ChebyshevPlan) + Icheb = let pinv = PI2{Float64}(Ω, D, 1_000, ChebyshevPlan) compute!(pinv, getpoints(f, pinv), 0, nothing) end - Ifindiff = let pinv = PI2{Float64}(Ω, D, (100, 100), FiniteDiffPlan) + Ifindiff = let pinv = PI2{Float64}(Ω, D, (31, 33), FiniteDiffPlan) compute!(pinv, getpoints(f, pinv), 0, nothing) end @@ -68,75 +68,89 @@ end @test Icheb ≈ Ifindiff rtol=1e-3 end -@safetestset "Free Particle" begin - using PoincareInvariants - - function free_particle!(points, t) - mid = length(points) ÷ 2 - for i in 1:mid - points[i] .+= points[mid+i] .* t - end - end - - @testset "Square" begin - D = 8 - N = 1_000 - Ω(z, t, p) = CanonicalSymplecticMatrix(D) - pinv = SecondPoincareInvariant{Float64}(Ω, D, N) - - phasepoints = getpoints(pinv) do x, y - (x, 0, 0, 0, y, 0, 0, 0) - end +@safetestset "Linear Symplectic Maps" begin + @safetestset "In 2D" begin + using PoincareInvariants + using PoincareInvariants.SecondPoincareInvariants: FiniteDiffPlan - @test abs(1 - compute!(pinv, phasepoints, 0, nothing)) / eps() < 10 + Ω(z, t, p) = CanonicalSymplecticMatrix(2) + pinv1 = PI2{Float64}(Ω, 2, 1_000) + pinv2 = PI2{Float64}(Ω, 2, 1_000, FiniteDiffPlan) - free_particle!(phasepoints, 10) - @test abs(1 - compute!(pinv, phasepoints, 0, nothing)) / eps() < 20 + A = [1 5; + 0 1] + B = [2 0; + 0 0.5] + f(r, θ) = A * B * [r*cospi(θ), r*sinpi(θ)] + # half circle with radius one - free_particle!(phasepoints, 100) - @test abs(1 - compute!(pinv, phasepoints, 0, nothing)) / eps() < 100 + @test abs(π/2 - compute!(pinv1, getpoints(f, pinv1), 0, nothing)) / eps() < 10 + @test abs(π/2 - compute!(pinv2, getpoints(f, pinv2), 0, nothing)) / π/2 < 1e-3 end - @testset "Quarter Circle" begin - D = 4 - N = 2_000 - Ω(z, t, p) = CanonicalSymplecticMatrix(D) - pinv = SecondPoincareInvariant{Float64}(Ω, D, N) + @safetestset "In 8D" begin + using PoincareInvariants + using PoincareInvariants.SecondPoincareInvariants: FiniteDiffPlan + + Ω(z, t, p) = CanonicalSymplecticMatrix(8) + pinv1 = PI2{Float64}(Ω, 8, 1_000) + pinv2 = PI2{Float64}(Ω, 8, 1_000, FiniteDiffPlan) + + A = CanonicalSymplecticMatrix(8) + + B = [1 0 0 0 1 5 8 1; + 0 1 0 0 5 2 6 9; + 0 0 1 0 8 6 3 7; + 0 0 0 1 1 9 7 4; + 0 0 0 0 1 0 0 0; + 0 0 0 0 0 1 0 0; + 0 0 0 0 0 0 1 0; + 0 0 0 0 0 0 0 1] + + C = [1.0 0 2.0 0 0 0 0 0; + 0 2.0 0 1.0 0 0 0 0; + 0 1.0 0 -2.0 0 0 0 0; + 2.0 0 -1.0 0 0 0 0 0; + 0 0 0 0 0.2 0 0.4 0; + 0 0 0 0 0 0.4 0 0.2; + 0 0 0 0 0.2 0 0 -0.4; + 0 0 0 0 0.4 0 -0.2 0] + + f(r, θ) = A * B * C * [r*sinpi(θ), 0, 0, 0, r*cospi(θ), 0, 0, 0] + # half circle with radius one and negative orientation + + @test abs(-π/2 - compute!(pinv1, getpoints(f, pinv1), 0, nothing)) / eps() < 10 + @test abs(-π/2 - compute!(pinv2, getpoints(f, pinv2), 0, nothing)) / π/2 < 1e-3 + end +end - phasepoints = getpoints(pinv) do r, θ - (0, r * cos(θ * π/2), 0, r * sin(θ * π/2)) - end +@safetestset "Henon Map" begin + using PoincareInvariants + using PoincareInvariants.SecondPoincareInvariants: FiniteDiffPlan - @test abs(π/4 - compute!(pinv, phasepoints, 0, nothing)) / eps() < 10 + Ω(z, t, p) = CanonicalSymplecticMatrix(4) + pinv1 = PI2{Float64}(Ω, 4, 1_000) + pinv2 = PI2{Float64}(Ω, 4, 1_000, FiniteDiffPlan) - free_particle!(phasepoints, 10) - @test abs(π/4 - compute!(pinv, phasepoints, 0, nothing)) / eps() < 20 + init(x, y) = (x, -x, y, -y) + henon((q, p), a) = (p + a - q^2, -q) - free_particle!(phasepoints, 100) - @test abs(π/4 - compute!(pinv, phasepoints, 0, nothing)) / eps() < 200 + # henon plus symplectic intermixing of q1/p1 with q2/p2 + function f((q1, q2, p1, p2)) + q1, p1 = henon((q1, p1), 2) + q2, p2 = henon((q2, p2), 3) + return (q1 - p2, q2 - p1, p1, p2) end - @testset "Half Sphere in 4D" begin - D = 4 - N = 1500 - Ω(z, t, p) = CanonicalSymplecticMatrix(D) - pinv = SecondPoincareInvariant{Float64}(Ω, D, N) - - phasepoints = getpoints(pinv) do θ, ϕ - sinθ, cosθ = sincospi(θ) - sinϕ, cosϕ = sincospi(ϕ) - return (sinθ * cosϕ, sinθ * sinϕ, cosθ, -1) - end + f1(x, y) = init(x, y) |> f + @test abs(2 - compute!(pinv1, getpoints(f1, pinv1), 0, nothing)) / eps() < 10 + @test abs(2 - compute!(pinv2, getpoints(f1, pinv2), 0, nothing)) / eps() < 500 - # invariant is -π - # TODO: do the calculation by hand and double check it + f2(x, y) = init(x, y) |> f |> f + @test abs(2 - compute!(pinv1, getpoints(f2, pinv1), 0, nothing)) / eps() < 50 + @test abs(2 - compute!(pinv2, getpoints(f2, pinv2), 0, nothing)) < 0.01 - @test abs(-π - compute!(pinv, phasepoints, 0, nothing)) / eps() < 15 - - free_particle!(phasepoints, 10) - @test abs(-π - compute!(pinv, phasepoints, 0, nothing)) / eps() < 15 - - free_particle!(phasepoints, 100) - @test abs(-π - compute!(pinv, phasepoints, 0, nothing)) / eps() < 150 - end + f3(x, y) = init(x, y) |> f |> f |> f + @test abs(2 - compute!(pinv1, getpoints(f3, pinv1), 0, nothing)) / eps() < 50 + @test abs(2 - compute!(pinv2, getpoints(f3, pinv2), 0, nothing)) < 0.5 end