Interoperability with SymPy#
Starting with version 0.10, symbolic expressions in heyoka.py can be converted to/from SymPy expressions. This feature can be very useful, as it unlocks the possibility of using the wide arsenal of SymPy algorithms on heyoka.py’s expressions. Possible use cases include:
the automatic simplification of heyoka.py’s expressions,
the application of symbolic algorithms not available in heyoka.py (e.g., symbolic integration),
LaTeX pretty-printing.
In this tutorial, we will show how SymPy interoperability can be used for expression simplification purposes, which ultimately leads to a measurable performance increase in the numerical integration of a system of ODEs.
For this example, we will consider the relativistic dynamics of Mercury’s orbit around the Sun, which, in the first post-Newtonian approximation, can be described by the following Hamiltonian:
where \(\left(v_x, v_y, v_z, x, y, z \right)\) is Mercury’s Cartesian state vector with respect to the Sun, \(\varepsilon = 1/c^2\), \(v=\sqrt{v_x^2+v_y^2+v_z^2}\), \(r=\sqrt{x^2+y^2+z^2}\) and \(\mu=GM_\odot\) is the gravitational parameter of the system. See the tutorial about Mercury’s relativistic precession for more details about this dynamical system.
Let us begin by defining the Hamiltonian as a symbolic expression in heyoka.py’s expression system:
import heyoka as hy
# Create the symbolic variables.
vx, vy, vz, x, y, z = hy.make_vars("v_x", "v_y", "v_z", "x", "y", "z")
# mu and eps constants.
mu = 0.01720209895 * 0.01720209895 * 365 * 365
eps = 2.5037803127808595e-10
# Auxiliary quantities.
v2 = vx*vx + vy*vy + vz*vz
r2 = x*x + y*y + z*z
r = hy.sqrt(r2)
# Initial conditions corresponding (roughly)
# to Mercury's current orbit.
ic = (-10.461611630114545,
6.667253667139184,
0.818635813947965,
0.16614243942411336,
0.2568228239702581,
0.0315338776710321)
# Define the Hamiltonian.
Ham = 1./2*v2 - mu/r + eps * (mu**2/(2*r2) - 1/8.*v2*v2 - 3./2.*mu*v2/r)
As usual, we are using Solar masses, astronomical units and years as units of measure.
Let us now convert the heyoka.py Hamiltonian to a SymPy expression via the to_sympy() function and pretty-print it to screen:
hy.to_sympy(Ham)
The SymPy expression can be converted back to heyoka.py using the from_sympy() function:
hy.from_sympy(hy.to_sympy(Ham))
((3.8912886205581639e-07 / ((2.0000000000000000 * x**2.0000000000000000) + (2.0000000000000000 * y**2.0000000000000000) + (2.0000000000000000 * z**2.0000000000000000))) + (0.50000000000000000 * v_x**2.0000000000000000) + (0.50000000000000000 * v_y**2.0000000000000000) + (0.50000000000000000 * v_z**2.0000000000000000) - (39.422903948847882 / (x**2.0000000000000000 + y**2.0000000000000000 + z**2.0000000000000000)**0.50000000000000000) - ((2.5037803127808595e-10 * ((59.134355923271826 * v_x**2.0000000000000000) + (59.134355923271826 * v_y**2.0000000000000000) + (59.134355923271826 * v_z**2.0000000000000000))) / (x**2.0000000000000000 + y**2.0000000000000000 + z**2.0000000000000000)**0.50000000000000000) - (2.5037803127808595e-10 * (v_x**2.0000000000000000 + v_y**2.0000000000000000 + v_z**2.0000000000000000) * ((0.12500000000000000 * v_x**2.0000000000000000) + (0.12500000000000000 * v_y**2.0000000000000000) + (0.12500000000000000 * v_z**2.0000000000000000))))
In order to formulate Hamilton’s equations for this system we need, as usual, to differentiate the Hamiltonian with respect to the state variables. Let us do it with heyoka.py’s diff() function:
# Formulate the equations of motion
# using heyoka.py's expression system.
ta = hy.taylor_adaptive(
# Hamilton's equations.
[(vx, -hy.diff(Ham, x)),
(vy, -hy.diff(Ham, y)),
(vz, -hy.diff(Ham, z)),
(x, hy.diff(Ham, vx)),
(y, hy.diff(Ham, vy)),
(z, hy.diff(Ham, vz))],
# Initial conditions.
ic
)
heyoka.py’s diff() function applies elementary calculus rules in order to differentiate a symbolic expression, performing only very basic simplifications in the process. We can get an estimate of the complexity of the resulting expressions by taking a look at the size of the Taylor decomposition for this dynamical system (this is the number of elementary subexpressions in which the equations of motion have been decomposed):
# NOTE: subtract 12 to avoid counting the 6 definitions of the state variables
# and the 6 definitions of their differential equations.
len(ta.decomposition) - 12
68
Let us now try to simplify Hamilton’s equations of motion using SymPy. We first define a small wrapper to invoke SymPy’s simplify() function on a heyoka.py expression:
import sympy
def simplify(ex):
return hy.from_sympy(sympy.simplify(hy.to_sympy(ex)))
Then, we create a new integrator with the simplified equations of motion:
ta_spy = hy.taylor_adaptive(
# Hamilton's equations.
[(vx, simplify(-hy.diff(Ham, x))),
(vy, simplify(-hy.diff(Ham, y))),
(vz, simplify(-hy.diff(Ham, z))),
(x, simplify(hy.diff(Ham, vx))),
(y, simplify(hy.diff(Ham, vy))),
(z, simplify(hy.diff(Ham, vz)))],
# Initial conditions.
ic
)
Let us take a look at the size of the Taylor decomposition for the new integrator:
len(ta_spy.decomposition) - 12
46
In other words, thanks to SymPy’s expression simplification capabilities, we reduced by \(\sim 30\%\) the symbolic complexity (i.e., the number of elementary subexpressions) of our ODE system.
What is the performance impact of these simplifications? Let us find out:
time ta_spy.propagate_until(10000.)
CPU times: user 2.33 s, sys: 0 ns, total: 2.33 s
Wall time: 2.33 s
(<taylor_outcome.time_limit: -4294967299>,
3.484957425280681e-06,
0.017389902845583543,
1030387,
None,
None)
time ta.propagate_until(10000.)
CPU times: user 2.38 s, sys: 0 ns, total: 2.38 s
Wall time: 2.38 s
(<taylor_outcome.time_limit: -4294967299>,
0.006522499652827606,
0.017389869779425565,
1012008,
None,
None)
The performance increase is measurable but not as large as the reduction in symbolic complexity.
As a word of warning, users should be aware that the general-purpose simplify() function from SymPy can be quite expensive from a computational point of view, especially for large expressions. SymPy provides various specialised simplification primitives that can be considerably more efficient.
More details on the conversion process#
When converting a heyoka.py expression to SymPy, the following conventions are adopted:
variables are converted to SymPy
Symbols,constants are converted to SymPy numbers,
runtime parameters are converted to SymPy
Symbols conventionally called"par[n]", wherenis the parameter index,functions are mapped to the corresponding SymPy functions.
Similarly, when converting a SymPy expression to heyoka.py, the following conventions are adopted:
SymPy
Symbols are converted to heyoka.py variables, unless the name of the symbol is"par[n]", in which case theSymbolis converted to a runtime parameter,SymPy numbers are converted to heyoka.py constants, if the conversion is lossless (otherwise an error is raised),
SymPy functions are mapped to heyoka.py functions.
When converting a SymPy expression to heyoka.py, it is possible to customise the conversion process by passing a substitution dictionary to the from_sympy() function. Here’s an example:
# Create a few SymPy symbols.
x, y, z = sympy.symbols("x y z")
# Convert the expression x*(y+z) to heyoka.py,
# but provide a custom substitution rule for
# the subexpression y+z.
hy.from_sympy(x * (y+z), {y+z: hy.expression("t")})
(t * x)
Note that the same effect could be achieved by doing the substitution in SymPy before the conversion to heyoka.py. However, for large expressions, using the custom substitution dictionary in the from_sympy() function will perform better.
Handling rational numbers#
heyoka.py is very conservative when converting from SymPy expressions containing rational numbers. In particular, if a rational number cannot be converted exactly to a floating-point value, heyoka.py will raise an error. In order to avoid this, you can use the sympy.N function to forcibly convert all rational numbers to floating-point values:
hy.from_sympy(sympy.N(x / 3))
(0.33333333333333331 * x)