@@ -9,68 +9,69 @@ x = [1.0, 2.0, 3.0]
9
9
p = [4.0 , 5.0 , 6.0 ]
10
10
11
11
prob = 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 )
13
13
14
14
# Parametric GTPSA map
15
15
desc = Descriptor (3 , 2 , 3 , 2 ) # 3 variables 3 parameters, both to 2nd order
16
16
dx = @vars (desc)
17
17
dp = @params (desc)
18
18
prob_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 )
20
20
21
21
@test sol. u[end ] ≈ scalar .(sol_GTPSA. u[end ]) # scalar gets 0th order part
22
22
23
23
# Compare Jacobian against ForwardDiff
24
24
J_FD = ForwardDiff. jacobian ([x... , p... ]) do t
25
25
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 )
27
27
sol. u[end ]
28
28
end
29
29
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 )
31
31
32
32
# Compare Hessians against ForwardDiff
33
33
for i in 1 : 3
34
34
Hi_FD = ForwardDiff. hessian ([x... , p... ]) do t
35
35
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 )
37
37
sol. u[end ][i]
38
38
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 )
40
40
end
41
41
42
-
43
42
# 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
48
49
49
50
prob = 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 )
51
52
52
53
desc = Descriptor (6 , 2 ) # 6 variables to 2nd order
53
- dx = @vars (desc) # identity map
54
+ dx = @vars (desc) # identity map
54
55
prob_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 )
56
57
57
58
@test sol. u[end ] ≈ scalar .(sol_GTPSA. u[end ]) # scalar gets 0th order part
58
59
59
60
# Compare Jacobian against ForwardDiff
60
61
J_FD = ForwardDiff. jacobian (zeros (6 )) do t
61
62
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 )
63
64
sol. u[end ]
64
65
end
65
66
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 )
67
68
68
69
# Compare Hessians against ForwardDiff
69
70
for i in 1 : 6
70
71
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 )
73
74
sol. u[end ][i]
74
75
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 )
76
77
end
0 commit comments