Fixes

Author: ParamThakkar123

benchmarks/DynamicalODE/Henon-Heiles_energy_conservation_benchmark.jmd

@@ -19,20 +19,22 @@ In order to separate these two, we will use `iip` for the in-place names and `oo
 using OrdinaryDiffEq, Plots, DiffEqCallbacks
 using SciMLBenchmarks
 using TaylorIntegration, LinearAlgebra, StaticArrays
+
+# Set plotting formats
 gr(fmt=:png)
 default(fmt=:png)
 
 T(p) = 1//2 * norm(p)^2
-V(q) = 1//2 * (q[1]^2 + q[2]^2 + 2q[1]^2 * q[2]- 2//3 * q[2]^3)
-H(p,q, params) = T(p) + V(q)
+V(q) = 1//2 * (q[1]^2 + q[2]^2 + 2 * q[1]^2 * q[2] - 2//3 * q[2]^3)
+H(p, q, params) = T(p) + V(q)
 
-function iip_dq(dq,p,q,params,t)
+function iip_dq(dq, p, q, params, t)
     dq[1] = p[1]
     dq[2] = p[2]
 end
 
-function iip_dp(dp,p,q,params,t)
-    dp[1] = -q[1] * (1 + 2q[2])
+function iip_dp(dp, p, q, params, t)
+    dp[1] = -q[1] * (1 + 2 * q[2])
     dp[2] = -q[2] - (q[1]^2 - q[2]^2)
 end
 
@@ -44,7 +46,7 @@ function oop_dq(p, q, params, t)
 end
 
 function oop_dp(p, q, params, t)
-    dp1 = -q[1] * (1 + 2q[2])
+    dp1 = -q[1] * (1 + 2 * q[2])
     dp2 = -q[2] - (q[1]^2 - q[2]^2)
     @SVector [dp1, dp2]
 end
@@ -56,19 +58,33 @@ function hamilton(du, u, p, t)
     dq, q = @views u[3:4], du[3:4]
     dp, p = @views u[1:2], du[1:2]
 
-    dp[1] = -q[1] * (1 + 2q[2])
+    dp[1] = -q[1] * (1 + 2 * q[2])
     dp[2] = -q[2] - (q[1]^2 - q[2]^2)
     dq .= p
 
     return nothing
 end
 
+function manifold_jacobian!(J, u)
+    p, q = u[1:2], u[3:4]
+    
+    J[1, 1] = p[1]
+    J[1, 2] = p[2]
+    
+    J[1, 3] = q[1] + 4 * q[1] * q[2]
+    J[1, 4] = q[2] + 2 * q[1]^2 - 2 * q[2]^2
+    
+    J[2:4, :] .= 0
+    
+    nothing
+end
+
 function g(resid, u, p)
     resid[1] = H([u[1], u[2]], [u[3], u[4]], nothing) - E
     resid[2:4] .= 0
 end
 
-const cb = ManifoldProjection(g, nlopts=Dict(:ftol => 1e-13), autodiff = AutoForwardDiff())
+const cb = ManifoldProjection(g, nlopts=Dict(:ftol => 1e-13), manifold_jacobian=manifold_jacobian)
 
 const E = H(iip_p0, iip_q0, nothing)
 ```

benchmarks/DynamicalODE/Project.toml

@@ -22,4 +22,4 @@ Plots = "2.0.0, 1"
 SciMLBenchmarks = "0.1.3"
 StaticArrays = "1.9.7"
 Statistics = "1.11.1"
-TaylorIntegration = "0.16.1"
+TaylorIntegration = "0.16.1"

benchmarks/DynamicalODE/Quadrupole_boson_Hamiltonian_energy_conservation_benchmark.jmd

@@ -29,6 +29,23 @@ function iip_dq(dq, p, q, params, t)
     dq[2] = A * p[2]
 end
 
+function manifold_jacobian!(J, u)
+    p₁, p₂ = u[1], u[2]
+    q₁, q₂ = u[3], u[4]
+    
+    q₁_sq = q₁^2
+    q₂_sq = q₂^2
+    q_sum_sq = q₁_sq + q₂_sq
+    
+    J[1,1] = A * p₁
+    J[1,2] = A * p₂
+    J[1,3] = A * q₁ + (B/√2) * (3q₂_sq - 3q₁_sq) + D * q₁ * q_sum_sq
+    J[1,4] = A * q₂ + (B/√2) * 6q₁q₂ + D * q₂ * q_sum_sq
+    J[2:4,:] .= 0
+    
+    nothing
+end
+
 function iip_dp(dp, p, q, params, t)
     dp[1] = -A * q[1] - 3 * B / sqrt(2) * (q[2]^2 - q[1]^2) - D * q[1] * (q[1]^2 + q[2]^2)
     dp[2] = -q[2] * (A + 3 * sqrt(2) * B * q[1] + D * (q[1]^2 + q[2]^2))
@@ -67,7 +84,7 @@ function g(resid, u, p)
 end
 
 const E = H(iip_p0, iip_q0, nothing)
-const cb = ManifoldProjection(g, nlopts=Dict(:ftol=>1e-13), autodiff = AutoForwardDiff());
+const cb = ManifoldProjection(g, nlopts=Dict(:ftol=>1e-13), manifold_jacobian! = manifold_jacobian!);
 ```
 
 For the comparison we will use the following function