diff --git a/doc/equations/harmonic_oscillator.png b/doc/equations/harmonic_oscillator.png
new file mode 100644
index 0000000000..427a3ea508
Binary files /dev/null and b/doc/equations/harmonic_oscillator.png differ
diff --git a/doc/equations/harmonic_oscillator.svg b/doc/equations/harmonic_oscillator.svg
new file mode 100644
index 0000000000..7320acb737
--- /dev/null
+++ b/doc/equations/harmonic_oscillator.svg
@@ -0,0 +1,32 @@
+
+
+
+
+
+
+
+
diff --git a/doc/html/quadrature.html b/doc/html/quadrature.html
index 9041b488ae..dadad954c4 100644
--- a/doc/html/quadrature.html
+++ b/doc/html/quadrature.html
@@ -51,6 +51,7 @@
Fourier Integrals Naive Monte Carlo Integration Symplectic Integration Wavelet Transforms Numerical Differentiation Automatic Differentiation
+ namespace boost{ namespace math{ namespace quadrature {
+
+ template
+ std::pair, std::vector > integrate_hamiltonian(const ReturnType p0,
+ const ReturnType q0,
+ const RealType dt,
+ const unsigned steps,
+ Func dHdp,
+ Func dHdq,
+ std::string method = "Y6")
+
+ template
+ std::pair, std::vector > integrate_hamiltonian(const ReturnType p0,
+ const ReturnType q0,
+ const RealType dt,
+ const unsigned steps,
+ Func dHdp,
+ Func dHdq,
+ std::string method,
+ const ``__Policy``& pol)
+
+ }}} // namespaces
+
+[heading Description]
+
+The functional `integrate_hamiltonian` calculates the phase space trajectory for a given Hamiltonian.
+The trajectories are calculated using symplectic integration which preserves the energy of
+a system. Even higher order traditional numerical integration algorithms will gain or lose
+energy at long times. The functional implements the methods in [@https://www.sciencedirect.com/science/article/abs/pii/0375960190900923 Yoshida's]
+landmark paper and methods by [@https://www.sciencedirect.com/science/article/pii/S0377042701004927?fr=RR-2&ref=pdf_download&rr=a1b51ec059e1fbec Blanes and Moan].
+We assume that the Hamiltonian is separable so that it can be written as `H = T(p) + V(q)`.
+
+We now give an example for a simple harmonic oscillator with the Hamiltonian
+[/ $H = \frac{p^2}{2m} + \frac{1}{2}kx^2$ ]
+[equation harmonic_oscillator]
+
+For simplicity we will assume `k = m = 1`. Then the partial derivatives of the Hamiltonian
+with respect to `p` and `q` are
+
+ double dHdp(double p)
+ {
+ return p;
+ }
+
+ double dHdq(double q)
+ {
+ return q;
+ }
+
+Note that the functional can be readily generalized to multiple coordinates by changing the
+signature of `dHdp` to
+
+ std::vector dHdp(std::vector p)
+ {
+ // calculate the partial derivatives with respect to each p_i
+ }
+
+We then define the timestep size and number of steps of the algorithm to go from `t=0` to `t=100`
+
+ RealType dt = 0.05;
+ RealType t_end = 100;
+ unsigned int steps = t_end / dt;
+
+Lastly, we define the initial conditions so that the oscillator starts from rest at `x=1`
+
+ RealType q0 = 1;
+ RealType p0 = 0;
+
+We now evolve the system using the following
+
+ std::vector p;
+ std::vector q;
+
+ std::tie(p, q) = boost::math::quadrature::integrate_hamiltonian(p0, q0, dt, steps, dHdp, dHdq, "Y6");
+
+The ouput vectors `q, p` are the position and momentum of the system at each time. In
+higher dimensions the output will be a vector of vectors.
+
+The last argument that we pass to `integrate_hamiltonian`, "Y6" here, sets the integration
+method to use. The string "Y6" stands for Yoshida's 6th order integrator. Yoshida's 2nd and
+4th order methods are also available by passing the string "Y2", or "Y4" respectively.
+The 4th order and 6th order symplectic Runge-Kutta-Nystrom methods from Table 3 of Blanes and Moan
+are available using the strings "SRKNB6" and "SRKNB11".
+
+[optional_policy]
+
+References:
+
+Yoshida, Haruo. [`Construction of higher order symplectic integrators], Physics Letters A, 150.5-7 (1990): 262-268.
+
+Blanes, Sergio, and Per Christian Moan. [`Practical symplectic partitioned runge–kutta and runge–kutta–nyström methods`], Journal of Computational and Applied Mathematics 142.2 (2002): 313-330.
+
+[endsect] [/section:symplectic Symplectic Integration]
+
diff --git a/include/boost/math/quadrature/symplectic.hpp b/include/boost/math/quadrature/symplectic.hpp
new file mode 100644
index 0000000000..3f82d69ec8
--- /dev/null
+++ b/include/boost/math/quadrature/symplectic.hpp
@@ -0,0 +1,340 @@
+// Copyright Jacob Hass, 2026
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_QUADRATURE_SYMPLECTIC_HPP
+#define BOOST_MATH_QUADRATURE_SYMPLECTIC_HPP
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+namespace boost{ namespace math { namespace quadrature { namespace detail {
+
+template
+using void_t = void;
+
+template
+struct has_plus : std::false_type {};
+
+template
+struct has_plus() + std::declval())> > : std::true_type {};
+
+template
+typename std::enable_if::value, T>::type
+add(T x, U y)
+{
+ return x + y;
+}
+
+template
+typename std::enable_if::value, T>::type
+add(T vec1, U vec2)
+{
+ for (unsigned i=0; i < vec1.size(); i++)
+ {
+ vec1[i] = vec1[i] + vec2[i];
+ }
+ return vec1;
+}
+
+template
+using void_t = void;
+
+template
+struct has_mult : std::false_type {};
+
+template
+struct has_mult() * std::declval())> > : std::true_type {};
+
+template
+typename std::enable_if::value, T>::type
+mult_prefactor(T x, U prefactor)
+{
+ return x * prefactor;
+}
+
+template
+typename std::enable_if::value, T>::type
+mult_prefactor(T vec1, U prefactor)
+{
+ for (unsigned i=0; i < vec1.size(); i++)
+ {
+ vec1[i] = vec1[i] * prefactor;
+ }
+ return vec1;
+}
+
+template
+std::pair second_order_yoshida(const RandomAccessContainer p0,
+ const RandomAccessContainer q0,
+ const RealType dt,
+ Func dHdp,
+ Func dHdq)
+{
+ RandomAccessContainer p = p0;
+ RandomAccessContainer q = q0;
+
+ // Half step in q
+ RandomAccessContainer dHdp_val = dHdp(p);
+ RandomAccessContainer dq = mult_prefactor(dHdp_val, 0.5 * dt);
+ q = add(q, dq);
+
+ // Full step in p
+ RandomAccessContainer dHdq_val = dHdq(q);
+ RandomAccessContainer dp = mult_prefactor(dHdq_val, -dt);
+ p = add(p, dp);
+
+ // Half step in q
+ dHdp_val = dHdp(p);
+ dq = mult_prefactor(dHdp_val, 0.5 * dt);
+ q = add(q, dq);
+
+ return std::make_pair(p, q);
+}
+
+template
+std::pair fourth_order_yoshida(const RandomAccessContainer p0,
+ const RandomAccessContainer q0,
+ const RealType dt,
+ Func dHdp,
+ Func dHdq)
+{
+ BOOST_MATH_STD_USING
+
+ RandomAccessContainer p = p0;
+ RandomAccessContainer q = q0;
+
+ // RealType x0 = -(std::pow(2, 1/3) / (2 - std::pow(2, 1/3)));
+ RealType x1 = 1. / (2. - std::cbrt(2.));
+ RealType x0 = 1. - 2. * x1;
+
+ std::vector weights = { x1, x0, x1 };
+
+ for (unsigned i=0; i < weights.size(); i++)
+ {
+ std::tie(p, q) = second_order_yoshida(p, q, weights[i] * dt, dHdp, dHdq);
+ }
+
+ return std::make_pair(p, q);
+}
+
+template
+std::pair sixth_order_yoshida(const RandomAccessContainer p0,
+ const RandomAccessContainer q0,
+ const RealType dt,
+ Func dHdp,
+ Func dHdq)
+{
+ RandomAccessContainer p = p0;
+ RandomAccessContainer q = q0;
+
+ // Choosing "System A" solution
+ // The following Mathematica command can calculate these coefficients to arbitrary precision
+ // FindRoot[{w0+2(w1+w2+w3)==1,
+ // w0^3 + 2(w1^3 + w2^3+w3^3)==0,
+ // w0^5 + 2(w1^5 + w2^5+w3^5)==0,
+ // 1/6(w0^2w1^3-w0^4*w1) - 1/6(w0^3w1^2-w0*w1^4)+1/6((w0+2w1)^2w2^3-(w0+2w1)(w0^3+2w1^3)w2) - 1/6((w0^3+2w1^3)w2^2-(w0+2w1)w2^4) +1/6((w0+2w1+2w2)^2w3^3-(w0+2w1+2w2)(w0^3+2w1^3+2w2^3)w3)-1/6((w0^3+2w1^3+2w2^3)w3^2-(w0+2w1+2w2)w3^4)==0}, {{w0,1.3151863206839063}, {w1,-1.17767998417887}, {w2,0.235573213359357}, {w3,0.784513610477560}}, WorkingPrecision->64]
+ RealType w1 = static_cast(-1.17767998417887100694641568096431573463926925263459848447536851379674155618156L);
+ RealType w2 = static_cast(0.23557321335935813368479318297853460168646808210340111900349313095621471215223L);
+ RealType w3 = static_cast(0.78451361047755726381949763386634987577682441745149338456794779895125997479548L);
+ // w0 = 1.31518632068391121888424972823886251435195350615940796180785516777853373846773
+ RealType w0 = 1. - 2. * (w1 + w2 + w3);
+ std::vector weights = { w3, w2, w1, w0, w1, w2, w3};
+
+ for (unsigned i=0; i < weights.size(); i++)
+ {
+ std::tie(p, q) = second_order_yoshida(p, q, weights[i] * dt, dHdp, dHdq);
+ }
+
+ return std::make_pair(p, q);
+}
+
+template
+std::pair SRKN_b_order_6(const RandomAccessContainer p0,
+ const RandomAccessContainer q0,
+ const RealType dt,
+ Func dHdp,
+ Func dHdq)
+{ // This method implements SRKN_b^6 in Table 3 here
+ // https://www.sciencedirect.com/science/article/pii/S0377042701004927
+
+ RandomAccessContainer p = p0;
+ RandomAccessContainer q = q0;
+
+ RealType b1 = static_cast(0.0829844064174052);
+ RealType b2 = static_cast(0.396309801498368);
+ RealType b3 = static_cast(-0.0390563049223486);
+ RealType b4 = 1. - 2. * (b1 + b2 + b3);
+
+ RealType a1 = static_cast(0.245298957184271);
+ RealType a2 = static_cast(0.604872665711080);
+ RealType a3 = 0.5 - (a1 + a2);
+
+ std::vector b_weights = {b1, b2, b3, b4, b3, b2};
+ std::vector a_weights = {a1, a2, a3, a3, a2, a1};
+
+ RealType a, b;
+ for (unsigned int i=0; i < b_weights.size(); i++)
+ {
+ b = b_weights[i];
+ a = a_weights[i];
+
+ RandomAccessContainer dHdp_val = dHdp(p);
+ RandomAccessContainer dq = mult_prefactor(dHdp_val, dt * b);
+ q = add(q, dq);
+
+ RandomAccessContainer dHdq_val = dHdq(q);
+ RandomAccessContainer dp = mult_prefactor(dHdq_val, -a * dt);
+ p = add(p, dp);
+ }
+ // Need to do one more step in q
+ RandomAccessContainer dHdp_val = dHdp(p);
+ RandomAccessContainer dq = mult_prefactor(dHdp_val, dt * b1);
+ q = add(q, dq);
+
+ return std::make_pair(p, q);
+}
+
+template
+std::pair SRKN_b_order_11(const RandomAccessContainer p0,
+ const RandomAccessContainer q0,
+ const RealType dt,
+ Func dHdp,
+ Func dHdq)
+{ // This method implements SRKN_b^11 in Table 3 here
+ // https://www.sciencedirect.com/science/article/pii/S0377042701004927
+
+ RandomAccessContainer p = p0;
+ RandomAccessContainer q = q0;
+
+ RealType b1 = static_cast(0.0414649985182624);
+ RealType b2 = static_cast(0.198128671918067);
+ RealType b3 = static_cast(-0.0400061921041533);
+ RealType b4 = static_cast(0.0752539843015807);
+ RealType b5 = static_cast(-0.0115113874206879);
+ RealType b6 = 0.5 - (b1 + b2 + b3 + b4 + b5);
+
+ RealType a1 = static_cast(0.123229775946271);
+ RealType a2 = static_cast(0.290553797799558);
+ RealType a3 = static_cast(-0.127049212625417);
+ RealType a4 = static_cast(-0.246331761062075);
+ RealType a5 = static_cast(0.357208872795928);
+ RealType a6 = 1. - 2. * (a1 + a2 + a3 + a4 + a5);
+
+ std::vector b_weights = {b1, b2, b3, b4, b5, b6, b6, b5, b4, b3, b2};
+ std::vector a_weights = {a1, a2, a3, a4, a5, a6, a5, a4, a3, a2, a1};
+
+ RealType a, b;
+ for (unsigned int i=0; i < b_weights.size(); i++)
+ {
+ b = b_weights[i];
+ a = a_weights[i];
+
+ RandomAccessContainer dHdp_val = dHdp(p);
+ RandomAccessContainer dq = mult_prefactor(dHdp_val, dt * b);
+ q = add(q, dq);
+
+ RandomAccessContainer dHdq_val = dHdq(q);
+ RandomAccessContainer dp = mult_prefactor(dHdq_val, -a * dt);
+ p = add(p, dp);
+ }
+ // Need to do one more step in q
+ RandomAccessContainer dHdp_val = dHdp(p);
+ RandomAccessContainer dq = mult_prefactor(dHdp_val, dt * b1);
+ q = add(q, dq);
+
+ return std::make_pair(p, q);
+}
+
+template
+std::pair, std::vector > integrate_hamiltonian_imp(const RandomAccessContainer p0,
+ const RandomAccessContainer q0,
+ const RealType dt,
+ const unsigned steps,
+ Func dHdp,
+ Func dHdq,
+ std::string method,
+ const Policy& pol)
+{
+ // Not sure how to make this function string nicer
+ static const char* function = "boost::math::quadrature::integrate_hamiltonian(p0, q0, %1%, steps, dHdp, dHdq)";
+
+ if ((dt <= 0) || !(boost::math::isfinite)(dt))
+ {
+ boost::math::policies::raise_domain_error(function, "Time step must be positive and finite but got: dt = %1%.\n", dt, pol);
+ }
+
+ // Check if method is available
+ std::vector available_methods = {"Y6", "Y4", "Y2"};
+
+ typedef std::pair (*stepperType)(RandomAccessContainer, RandomAccessContainer, RealType, Func, Func);
+
+ std::map m{{"Y6", sixth_order_yoshida},
+ {"Y4", fourth_order_yoshida},
+ {"Y2", second_order_yoshida},
+ {"SRKNB6", SRKN_b_order_6},
+ {"SRKNB11", SRKN_b_order_11}};
+ stepperType stepper = m.at(method);
+
+ std::vector p(steps);
+ std::vector q(steps);
+ p[0] = p0;
+ q[0] = q0;
+
+ RandomAccessContainer p_current = p0;
+ RandomAccessContainer q_current = q0;
+ for (unsigned i=1; i < steps; i++)
+ {
+ std::tie(p_current, q_current) = stepper(p_current, q_current, dt, dHdp, dHdq);
+ p[i] = p_current;
+ q[i] = q_current;
+ }
+ return std::make_pair(p, q);
+}
+} // namespace detail
+
+template
+std::pair, std::vector > integrate_hamiltonian(const RandomAccessContainer p0,
+ const RandomAccessContainer q0,
+ const RealType dt,
+ const unsigned steps,
+ Func dHdp,
+ Func dHdq,
+ std::string method,
+ const Policy& pol)
+{
+ return detail::integrate_hamiltonian_imp(p0, q0, dt, steps, dHdp, dHdq, method, pol);
+}
+
+template
+std::pair, std::vector > integrate_hamiltonian(const RandomAccessContainer p0,
+ const RandomAccessContainer q0,
+ const RealType dt,
+ const unsigned steps,
+ Func dHdp,
+ Func dHdq,
+ std::string method)
+{
+ return integrate_hamiltonian(p0, q0, dt, steps, dHdp, dHdq, method, boost::math::policies::policy<>());
+}
+
+template
+std::pair, std::vector > integrate_hamiltonian(const RandomAccessContainer p0,
+ const RandomAccessContainer q0,
+ const RealType dt,
+ const unsigned steps,
+ Func dHdp,
+ Func dHdq)
+{
+ return integrate_hamiltonian(p0, q0, dt, steps, dHdp, dHdq, "Y6", boost::math::policies::policy<>());
+}
+}}}
+#endif
diff --git a/test/Jamfile.v2 b/test/Jamfile.v2
index 18d5d9ae01..b19a34ff27 100644
--- a/test/Jamfile.v2
+++ b/test/Jamfile.v2
@@ -1343,6 +1343,7 @@ test-suite quadrature :
[ run test_trapezoidal.cpp /boost/test//boost_unit_test_framework : : :
release [ requires cxx11_lambdas cxx11_auto_declarations cxx11_decltype cxx11_unified_initialization_syntax cxx11_variadic_templates ]
$(float128_type) ]
+ [ run test_symplectic.cpp ]
;
test-suite autodiff :
diff --git a/test/test_symplectic.cpp b/test/test_symplectic.cpp
new file mode 100644
index 0000000000..8be7cf102b
--- /dev/null
+++ b/test/test_symplectic.cpp
@@ -0,0 +1,228 @@
+/*
+ * Copyright Jacob Hass, 2026
+ * Use, modification and distribution are subject to the
+ * Boost Software License, Version 1.0. (See accompanying file
+ * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+#define BOOST_TEST_MODULE symplectic_quadrature
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#if __has_include()
+# include
+#endif
+
+using boost::math::quadrature::integrate_hamiltonian;
+
+// Equations of motion for simple harmonic oscillator
+template
+Real oscillator_dHdp(Real p) { return p; }
+template
+Real oscillator_dHdq(Real q) { return q; }
+
+// Equations of motion for simple pendulum
+template
+std::vector pendulum_vector_dHdp(std::vector p) { return p; }
+template
+std::vector pendulum_vector_dHdq(std::vector q)
+{
+ BOOST_MATH_STD_USING
+ std::vector q_return(q.size());
+ for (unsigned i=0; i
+std::vector hh_dHdp(std::vector p)
+{
+ return p;
+}
+
+template
+std::vector hh_dHdq(std::vector q)
+{
+ BOOST_MATH_STD_USING
+ std::vector dHdq(q.size());
+ dHdq[0] = q[0] + 2.0 * q[0] * q[1]; // x + 2*xy
+ dHdq[1] = q[1] + pow(q[0], 2.0) - pow(q[1], 2.0); // y + x^2 - y^2
+ return dHdq;
+}
+
+// The energy of the system is H = 1/2 (px^2 + py^2) + 1/2(x^2 + y^2) + x^2 y - y^3/3
+template
+Real hh_energy(std::vector p, std::vector q)
+{
+ BOOST_MATH_STD_USING
+ Real x = q[0];
+ Real y = q[1];
+ Real px = p[0];
+ Real py = p[1];
+ return 0.5 * (pow(px, 2.) + pow(py, 2.)) + 0.5 * (pow(x, 2.) + pow(y, 2.)) + pow(x, 2.) * y - pow(y, 3.) / 3.;
+}
+
+template
+void test_invalid_parameters()
+{
+ RealType q0 = 1;
+ RealType p0 = 0;
+ // Negative timestep
+ BOOST_CHECK_THROW(boost::math::quadrature::integrate_hamiltonian(q0, p0, -0.1, 10, oscillator_dHdp, oscillator_dHdq), std::domain_error);
+
+ // Method not in {'Y6', 'Y4', 'Y2'}
+ BOOST_CHECK_THROW(boost::math::quadrature::integrate_hamiltonian(q0, p0, 0.1, 10, oscillator_dHdp, oscillator_dHdq, "InvalidMethod"), std::out_of_range);
+}
+
+/* Test if SHO energy fluctuations are below a given tolerance*/
+template
+void test_harmonic_oscillator(const RealType tol, const std::string method)
+{
+ BOOST_MATH_STD_USING
+
+ RealType dt = 0.05;
+ RealType t_end = 100;
+ unsigned int steps = t_end / dt;
+
+ RealType q0 = 1;
+ RealType p0 = 0;
+
+ std::vector p;
+ std::vector q;
+
+ std::tie(p, q) = boost::math::quadrature::integrate_hamiltonian(p0, q0, dt, steps, oscillator_dHdp, oscillator_dHdq, method);
+
+ RealType p_val;
+ RealType q_val;
+ std::vector abs_energy_error(p.size());
+ for (unsigned i=0; i < p.size(); i++)
+ {
+ p_val = p[i];
+ q_val = q[i];
+
+ abs_energy_error[i] = std::abs(std::pow(p_val, 2) + std::pow(q_val, 2) - 1);
+ }
+
+ RealType max_error = *std::max_element(std::begin(abs_energy_error), std::end(abs_energy_error));
+ BOOST_CHECK_LE(max_error, tol);
+}
+
+/* Test if SHO energy fluctuations are below a given tolerance*/
+template
+void test_pendulum(const RealType tol, const std::string method)
+{
+ BOOST_MATH_STD_USING
+
+ RealType dt = 0.05;
+ RealType t_end = 100;
+ unsigned int steps = t_end / dt;
+
+ std::vector q0 = {boost::math::constants::pi() * 0.5};
+ std::vector p0 = {0.};
+
+ std::vector > p;
+ std::vector > q;
+
+ std::tie(p, q) = boost::math::quadrature::integrate_hamiltonian(p0, q0, dt, steps, pendulum_vector_dHdp, pendulum_vector_dHdq, method);
+
+ RealType p_val;
+ RealType q_val;
+ std::vector abs_energy_error(p.size());
+ for (unsigned i=0; i < p.size(); i++)
+ {
+ p_val = p[i][0];
+ q_val = q[i][0];
+ // The energy of the system is H = p^2 / 2 + (1-cos(q)). With our initial condition
+ // this yields H = 1. Thus, we just want p^2 / 2 - cos(q) = 0
+ abs_energy_error[i] = std::abs(std::pow(p_val, 2.) / 2. - std::cos(q_val));
+ }
+
+ RealType max_error = *std::max_element(std::begin(abs_energy_error), std::end(abs_energy_error));
+ BOOST_CHECK_LE(max_error, tol);
+}
+
+/* Test if SHO energy fluctuations are below a given tolerance*/
+template
+void test_hh_model(const RealType tol, const std::string method)
+{
+ BOOST_MATH_STD_USING
+
+ RealType dt = 0.005;
+ unsigned int steps = 20000;
+
+ std::vector q0 = {0.5, 0};
+ std::vector p0 = {0, 0.25};
+
+ RealType total_energy = hh_energy(p0, q0);
+
+ std::vector > p;
+ std::vector > q;
+
+ std::tie(p, q) = boost::math::quadrature::integrate_hamiltonian(p0, q0, dt, steps, hh_dHdp, hh_dHdq, method);
+
+ std::vector p_val;
+ std::vector q_val;
+ std::vector abs_energy_error(p.size());
+ RealType sum = 0;
+ for (unsigned i=0; i < p.size(); i++)
+ {
+ p_val = p[i];
+ q_val = q[i];
+
+ abs_energy_error[i] = abs(hh_energy(p_val, q_val) - total_energy);
+ sum += abs_energy_error[i];
+ }
+
+ RealType max_error = *std::max_element(std::begin(abs_energy_error), std::end(abs_energy_error));
+ BOOST_CHECK_LE(max_error, tol);
+}
+
+BOOST_AUTO_TEST_CASE(symplectic_quadrature)
+{
+ test_invalid_parameters();
+
+ // Test doubles
+ // Simple Harmonic Oscillator Tests
+ test_harmonic_oscillator(1e-10, "Y6");
+ test_harmonic_oscillator(7e-4, "Y4");
+ test_harmonic_oscillator(7e-4, "Y2");
+ test_harmonic_oscillator(1e-11, "SRKNB6");
+ test_harmonic_oscillator(2e-14, "SRKNB11");
+
+ // Pendulum Tests
+ test_pendulum(1e-10, "Y6");
+ test_pendulum(5e-4, "Y4");
+ test_pendulum(5e-4, "Y2");
+ test_pendulum(1e-8, "SRKNB6");
+ test_pendulum(1e-10, "SRKNB11");
+
+ // Henon Heiles Model
+ test_hh_model(1e-14, "SRKNB11");
+ test_hh_model(1e-14, "Y6");
+ test_hh_model(5e-11, "Y4");
+ test_hh_model(5e-6, "Y2");
+ test_hh_model(1e-13, "SRKNB6");
+
+ // Test floats
+ test_harmonic_oscillator(6e-6, "Y6");
+ test_pendulum(5e-6, "Y6");
+ test_hh_model(5e-6, "Y6");
+
+ // Test long doubles
+ test_harmonic_oscillator(1e-10, "Y6");
+ test_pendulum(1e-10, "Y6");
+ test_hh_model(1e-14, "Y6");
+
+ // Test multiprecision
+ test_hh_model(1e-16, "Y6");
+}