|
| 1 | +using OrdinaryDiffEq, ForwardDiff, GTPSA, Test |
| 2 | + |
| 3 | +f!(du, u, p, t) = du .= p .* u |
| 4 | + |
| 5 | +# Initial variables and parameters |
| 6 | +x = [1.0, 2.0, 3.0] |
| 7 | +p = [4.0, 5.0, 6.0] |
| 8 | + |
| 9 | +prob = ODEProblem(f!, x, (0.0, 1.0), p) |
| 10 | +sol = solve(prob, Tsit5(), reltol=1e-16, abstol=1e-16) |
| 11 | + |
| 12 | +# Parametric GTPSA map |
| 13 | +desc = Descriptor(3, 2, 3, 2) # 3 variables 3 parameters, both to 2nd order |
| 14 | +dx = vars(desc) |
| 15 | +dp = params(desc) |
| 16 | +prob_GTPSA = ODEProblem(f!, x .+ dx, (0.0, 1.0), p .+ dp) |
| 17 | +sol_GTPSA = solve(prob_GTPSA, Tsit5(), reltol=1e-16, abstol=1e-16) |
| 18 | + |
| 19 | +@test sol.u[end] ≈ scalar.(sol_GTPSA.u[end]) # scalar gets 0th order part |
| 20 | + |
| 21 | +# Compare Jacobian against ForwardDiff |
| 22 | +J_FD = ForwardDiff.jacobian([x..., p...]) do t |
| 23 | + prob = ODEProblem(f!, t[1:3], (0.0, 1.0), t[4:6]) |
| 24 | + sol = solve(prob, Tsit5(), reltol=1e-16, abstol=1e-16) |
| 25 | + sol.u[end] |
| 26 | +end |
| 27 | + |
| 28 | +@test J_FD ≈ GTPSA.jacobian(sol_GTPSA.u[end], include_params=true) |
| 29 | + |
| 30 | +# Compare Hessians against ForwardDiff |
| 31 | +for i in 1:3 |
| 32 | + Hi_FD = ForwardDiff.hessian([x..., p...]) do t |
| 33 | + prob = ODEProblem(f!, t[1:3], (0.0, 1.0), t[4:6]) |
| 34 | + sol = solve(prob, Tsit5(), reltol=1e-16, abstol=1e-16) |
| 35 | + sol.u[end][i] |
| 36 | + end |
| 37 | + @test Hi_FD ≈ GTPSA.hessian(sol_GTPSA.u[end][i], include_params=true) |
| 38 | +end |
| 39 | + |
0 commit comments