@@ -62,7 +62,7 @@ for i=2:Ntraj
6262end
6363=#
6464
65- using OrdinaryDiffEq, DiffEqCallbacks, Test
65+ using OrdinaryDiffEq, NonlinearSolve, DiffEqCallbacks, Test
6666
6767# Initial state
6868u0 = [0 , - 0.25 , 0.42081 , 0 ]
8484const E = Hhh (u0)
8585
8686function ghh (resid, u, p)
87- resid[1 ] = Hhh (u[1 ], u[2 ], u[3 ], u[4 ]) - E
88- resid[2 : 4 ] .= 0
87+ resid[1 ] = - Hhh (u[1 ], u[2 ], u[3 ], u[4 ]) + E
8988end
9089
9190# energy conserving callback:
9291# important to use save = false, I don't want rescaling points
93- cb = ManifoldProjection (ghh, abstol = 1e-13 , save = false , autodiff = AutoForwardDiff ())
92+ cb = ManifoldProjection (ghh, resid_prototype = ones ( 1 ), nlsolve = TrustRegion (), abstol = 1e-9 , save = false , autodiff = AutoForwardDiff ())
9493
9594# Callback for Poincare surface of section
9695function psos_callback (j, direction = + 1 , offset:: Real = 0 ,
@@ -113,7 +112,7 @@ totalcb = CallbackSet(poincarecb, cb)
113112prob = ODEProblem (hheom!, u0, (0.0 , 100.0 ), callback = totalcb)
114113
115114extra_kw = Dict (:save_start => false , :save_end => false )
116- DEFAULT_DIFFEQ_KWARGS = Dict {Symbol, Any} (:abstol => 1e-9 , :reltol => 1e-9 )
115+ DEFAULT_DIFFEQ_KWARGS = Dict {Symbol, Any} (:abstol => 1e-10 , :reltol => 1e-10 )
117116
118117sol = solve (prob, Vern9 (); extra_kw... , DEFAULT_DIFFEQ_KWARGS... , save_everystep = false )
119118
0 commit comments