Updates

Author: ParamThakkar123

benchmarks/DynamicalODE/Henon-Heiles_energy_conservation_benchmark.jmd

@@ -63,12 +63,27 @@ function hamilton(du, u, p, t)
     return nothing
 end
 
+function custom_manifold_jacobian(u, params)
+    p1, p2, q1, q2 = u
+
+    J = zeros(4, 4)
+
+    J[1, 1] = 1.0
+    J[2, 2] = 1.0 
+    J[3, 3] = 1 + 2 * q_2
+    J[3, 4] = 2 * q1    
+    J[4, 3] = 2 * q1      
+    J[4, 4] = 1 - 4 * q2
+
+    return J
+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))
+const cb = ManifoldProjection(g, nlopts=Dict(:ftol => 1e-13), manifold_jacobian = custom_manifold_jacobian)
 const E = H(iip_p0, iip_q0, nothing)
 ```
 

benchmarks/DynamicalODE/Quadrupole_boson_Hamiltonian_energy_conservation_benchmark.jmd

@@ -61,13 +61,27 @@ function hamilton(z, params, t)
     )
 end
 
+function custom_manifold_jacobian(u, params)
+    p0, p2, q0, q2 = u
+
+    J = zeros(4, 4)
+
+    J[1, 1] = A * p0 
+    J[2, 2] = A * p2
+
+    J[3, 3] = A * q0 - sqrt(2) * B * (3 * q2^2 - q0^2) + D * q0 * (q0^2 + q2^2)  # ∂H/∂q0
+    J[4, 4] = A * q2 + 6 * B * q0 * q2 + D * q2 * (q0^2 + q2^2)                   # ∂H/∂q2
+
+    return J
+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 E = H(iip_p0, iip_q0, nothing)
-const cb = ManifoldProjection(g, nlopts=Dict(:ftol=>1e-13));
+const cb = ManifoldProjection(g, nlopts=Dict(:ftol=>1e-13), manifold_jacobian = custom_manifold_jacobian);
 ```
 
 For the comparison we will use the following function