fix + test trajectory (#80)

Author: Datseris

src/continuous.jl

@@ -8,8 +8,6 @@ export CDS_KWARGS
 using SimpleDiffEq: SimpleATsit5
 const DEFAULT_SOLVER = SimpleATsit5()
 const DEFAULT_DIFFEQ_KWARGS = (abstol = 1e-6, reltol = 1e-6)
-
-
 const CDS_KWARGS = (alg = DEFAULT_SOLVER, DEFAULT_DIFFEQ_KWARGS...)
 
 _get_solver(a) = haskey(a, :alg) ? a[:alg] : DEFAULT_SOLVER
@@ -132,7 +130,7 @@ function trajectory(ds::ContinuousDynamicalSystem, T, u = ds.u0;
     t0 = ds.t0
     tvec = (t0+Ttr):dt:(T+t0+Ttr)
     sol = Vector{SVector{dimension(ds), stateeltype(ds)}}(undef, length(tvec))
-    integ = integrator(ds, u; dt = dt, tfinal = tvec[end]+2dt, diffeq...)
+    integ = integrator(ds, u; dt = dt, diffeq...)
     step!(integ, Ttr)
     for (i, t) in enumerate(tvec)
         while t > integ.t

test/continuous_systems.jl

@@ -6,17 +6,16 @@ println("\nTesting continuous system evolution...")
 
 @testset "Lorenz System" begin
 
-  lo11 = Systems.lorenz() #with Jac
-  # lo33 = Systems.lorenz(big.([0.0, 10.0, 0.0]))
-
-  @testset "trajectory" begin
+    lo11 = Systems.lorenz() #with Jac
+    # lo33 = Systems.lorenz(big.([0.0, 10.0, 0.0]))
     # trajectory pure:
     ts1 = trajectory(lo11, 1.0, dt = 0.01)
     # ts3 = trajectory(lo33, 1.0)
     # @test eltype(ts3[1]) == BigFloat
     @test size(ts1) == (101, 3)
     # @test ts1[1, :] ≈ ts3[end,:]
     @test ts1[1, :] == SVector{3}(lo11.u0)
+    @test !any(x -> abs(x) > 1e3, ts1[end])
     # trajectory with diff_eq_kwargs and dt:
     ts1 = trajectory(lo11, 1.0; dt=0.1,
     diff_eq_kwargs=Dict(:abstol=>1e-8, :reltol=>1e-8))
@@ -25,32 +24,29 @@ println("\nTesting continuous system evolution...")
     @test size(ts1) == (11, 3)
     @test ts1[1, :] == SVector{3}(lo11.u0)
     # @test ts1[end, :] ≈ ts3[end,:]
-  end
 
 end
 
 @testset "Lorenz96" begin
     u = ones(5)
     lo11 = Systems.lorenz96(5, u)
     # lo33 = Systems.lorenz96(5, big.(ones(5)))
-
-    @testset "trajectory" begin
-      # trajectory pure:
-      ts1 = trajectory(lo11, 2.0, dt = 0.01)
-      # ts3 = trajectory(lo33, 2.0)
-      # @test eltype(ts3[1]) == BigFloat
-      @test size(ts1) == (201, 5)
-      # @test ts1[end, :] ≈ ts3[end,:]
-      @test ts1[1, :] == SVector{5}(u)
-      # trajectory with diff_eq_kwargs and dt:
-      ts1 = trajectory(lo11, 2.0; dt=0.1,
-      diff_eq_kwargs=Dict(:abstol=>1e-9, :reltol=>1e-9))
-      # ts3 = trajectory(lo33, 2.0; dt=0.1,
-      # diff_eq_kwargs=Dict(:abstol=>1e-9, :reltol=>1e-9))
-      @test size(ts1) == (21, 5)
-      # @test ts1[end, :] ≈ ts3[end,:]
-      @test ts1[1, :] == SVector{5}(u)
-    end
+    # trajectory pure:
+    ts1 = trajectory(lo11, 2.0, dt = 0.01)
+    # ts3 = trajectory(lo33, 2.0)
+    # @test eltype(ts3[1]) == BigFloat
+    @test size(ts1) == (201, 5)
+    # @test ts1[end, :] ≈ ts3[end,:]
+    @test ts1[1, :] == SVector{5}(u)
+    # trajectory with diff_eq_kwargs and dt:
+    ts1 = trajectory(lo11, 2.0; dt=0.1,
+    diff_eq_kwargs=Dict(:abstol=>1e-9, :reltol=>1e-9))
+    # ts3 = trajectory(lo33, 2.0; dt=0.1,
+    # diff_eq_kwargs=Dict(:abstol=>1e-9, :reltol=>1e-9))
+    @test size(ts1) == (21, 5)
+    # @test ts1[end, :] ≈ ts3[end,:]
+    @test ts1[1, :] == SVector{5}(u)
+    @test !any(x -> abs(x) > 1e3, ts1[end])
 end
 
 @testset "Gissinger Columns" begin
@@ -69,4 +65,14 @@ end
     ds = Systems.qbh()
     tr = trajectory(ds, 20.0)
     @test size(tr) == (2001, 4)
+    @test !any(x -> abs(x) > 1e3, tr[end])
   end
+
+@testset "Finite trajectory HH" begin
+    ds = Systems.henonheiles()
+    u0 = [0, -0.25, 0.42081, 0]
+    tr = trajectory(ds, 100.0, u0; dt = 0.05)
+    @test !any(x -> abs(x) > 1e3, tr[end])
+    tr = trajectory(ds, 100.0, u0; dt = 0.05, alg = Vern9())
+    @test !any(x -> abs(x) > 1e3, tr[end])
+end