@@ -9,68 +9,69 @@ x = [1.0, 2.0, 3.0]
99p = [4.0 , 5.0 , 6.0 ]
1010
1111prob = ODEProblem (f!, x, (0.0 , 1.0 ), p)
12- sol = solve (prob, Tsit5 (), reltol= 1e-16 , abstol= 1e-16 )
12+ sol = solve (prob, Tsit5 (), reltol = 1e-16 , abstol = 1e-16 )
1313
1414# Parametric GTPSA map
1515desc = Descriptor (3 , 2 , 3 , 2 ) # 3 variables 3 parameters, both to 2nd order
1616dx = @vars (desc)
1717dp = @params (desc)
1818prob_GTPSA = ODEProblem (f!, x .+ dx, (0.0 , 1.0 ), p .+ dp)
19- sol_GTPSA = solve (prob_GTPSA, Tsit5 (), reltol= 1e-16 , abstol= 1e-16 )
19+ sol_GTPSA = solve (prob_GTPSA, Tsit5 (), reltol = 1e-16 , abstol = 1e-16 )
2020
2121@test sol. u[end ] ≈ scalar .(sol_GTPSA. u[end ]) # scalar gets 0th order part
2222
2323# Compare Jacobian against ForwardDiff
2424J_FD = ForwardDiff. jacobian ([x... , p... ]) do t
2525 prob = ODEProblem (f!, t[1 : 3 ], (0.0 , 1.0 ), t[4 : 6 ])
26- sol = solve (prob, Tsit5 (), reltol= 1e-16 , abstol= 1e-16 )
26+ sol = solve (prob, Tsit5 (), reltol = 1e-16 , abstol = 1e-16 )
2727 sol. u[end ]
2828end
2929
30- @test J_FD ≈ GTPSA. jacobian (sol_GTPSA. u[end ], include_params= true )
30+ @test J_FD ≈ GTPSA. jacobian (sol_GTPSA. u[end ], include_params = true )
3131
3232# Compare Hessians against ForwardDiff
3333for i in 1 : 3
3434 Hi_FD = ForwardDiff. hessian ([x... , p... ]) do t
3535 prob = ODEProblem (f!, t[1 : 3 ], (0.0 , 1.0 ), t[4 : 6 ])
36- sol = solve (prob, Tsit5 (), reltol= 1e-16 , abstol= 1e-16 )
36+ sol = solve (prob, Tsit5 (), reltol = 1e-16 , abstol = 1e-16 )
3737 sol. u[end ][i]
3838 end
39- @test Hi_FD ≈ GTPSA. hessian (sol_GTPSA. u[end ][i], include_params= true )
39+ @test Hi_FD ≈ GTPSA. hessian (sol_GTPSA. u[end ][i], include_params = true )
4040end
4141
42-
4342# ODEProblem 2 =======================
44- pdot! (dq, p, q, params, t) = dq .= [0.0 , 0.0 , 0.0 ]
45- qdot! (dp, p, q, params, t) = dp .= [p[1 ] / sqrt ((1 + p[3 ])^ 2 - p[1 ]^ 2 - p[2 ]^ 2 ),
46- p[2 ] / sqrt ((1 + p[3 ])^ 2 - p[1 ]^ 2 - p[2 ]^ 2 ),
47- p[3 ] / sqrt (1 + p[3 ]^ 2 ) - (p[3 ] + 1 )/ sqrt ((1 + p[3 ])^ 2 - p[1 ]^ 2 - p[2 ]^ 2 )]
43+ pdot! (dq, p, q, params, t) = dq .= [0.0 , 0.0 , 0.0 ]
44+ function qdot! (dp, p, q, params, t)
45+ dp .= [p[1 ] / sqrt ((1 + p[3 ])^ 2 - p[1 ]^ 2 - p[2 ]^ 2 ),
46+ p[2 ] / sqrt ((1 + p[3 ])^ 2 - p[1 ]^ 2 - p[2 ]^ 2 ),
47+ p[3 ] / sqrt (1 + p[3 ]^ 2 ) - (p[3 ] + 1 ) / sqrt ((1 + p[3 ])^ 2 - p[1 ]^ 2 - p[2 ]^ 2 )]
48+ end
4849
4950prob = DynamicalODEProblem (pdot!, qdot!, [0.0 , 0.0 , 0.0 ], [0.0 , 0.0 , 0.0 ], (0.0 , 25.0 ))
50- sol = solve (prob, Yoshida6 (), dt = 1.0 , reltol= 1e-16 , abstol= 1e-16 )
51+ sol = solve (prob, Yoshida6 (), dt = 1.0 , reltol = 1e-16 , abstol = 1e-16 )
5152
5253desc = Descriptor (6 , 2 ) # 6 variables to 2nd order
53- dx = @vars (desc) # identity map
54+ dx = @vars (desc) # identity map
5455prob_GTPSA = DynamicalODEProblem (pdot!, qdot!, dx[1 : 3 ], dx[4 : 6 ], (0.0 , 25.0 ))
55- sol_GTPSA = solve (prob_GTPSA, Yoshida6 (), dt = 1.0 , reltol= 1e-16 , abstol= 1e-16 )
56+ sol_GTPSA = solve (prob_GTPSA, Yoshida6 (), dt = 1.0 , reltol = 1e-16 , abstol = 1e-16 )
5657
5758@test sol. u[end ] ≈ scalar .(sol_GTPSA. u[end ]) # scalar gets 0th order part
5859
5960# Compare Jacobian against ForwardDiff
6061J_FD = ForwardDiff. jacobian (zeros (6 )) do t
6162 prob = DynamicalODEProblem (pdot!, qdot!, t[1 : 3 ], t[4 : 6 ], (0.0 , 25.0 ))
62- sol = solve (prob, Yoshida6 (), dt = 1.0 , reltol= 1e-16 , abstol= 1e-16 )
63+ sol = solve (prob, Yoshida6 (), dt = 1.0 , reltol = 1e-16 , abstol = 1e-16 )
6364 sol. u[end ]
6465end
6566
66- @test J_FD ≈ GTPSA. jacobian (sol_GTPSA. u[end ], include_params= true )
67+ @test J_FD ≈ GTPSA. jacobian (sol_GTPSA. u[end ], include_params = true )
6768
6869# Compare Hessians against ForwardDiff
6970for i in 1 : 6
7071 Hi_FD = ForwardDiff. hessian (zeros (6 )) do t
71- prob = DynamicalODEProblem (pdot!, qdot!, t[1 : 3 ], t[4 : 6 ], (0.0 , 25.0 ))
72- sol = solve (prob, Yoshida6 (), dt = 1.0 , reltol= 1e-16 , abstol= 1e-16 )
72+ prob = DynamicalODEProblem (pdot!, qdot!, t[1 : 3 ], t[4 : 6 ], (0.0 , 25.0 ))
73+ sol = solve (prob, Yoshida6 (), dt = 1.0 , reltol = 1e-16 , abstol = 1e-16 )
7374 sol. u[end ][i]
7475 end
75- @test Hi_FD ≈ GTPSA. hessian (sol_GTPSA. u[end ][i], include_params= true )
76+ @test Hi_FD ≈ GTPSA. hessian (sol_GTPSA. u[end ][i], include_params = true )
7677end
0 commit comments