[FEATURE] component/variation-level tolerance specification

[alangfor]

Feature Request

Is your feature request related to a problem?
Propagation of higher-order variations can drive the integrator to reduce its time-step in order to satisfy a uniform numerical tolerance, overshadowing the baseline state’s needs.

Describe the solution you’d like
The ability to specify unique integration tolerances for different components or variation orders. For example, one might want a tighter tolerance (e.g., 1e-14) for the baseline state but a looser tolerance (e.g., 1e-8) for first-order variations.

Describe alternatives you’ve considered
Lowering the overall tolerance lets the integrator take larger steps for the variational equations but reduces the accuracy of the baseline state.

Additional context
SciPy’s solve_ivp allows per-component tolerances (via an array for atol), but this may not map easily to Heyoka. A dedicated approach in hy.var_ode_sys could allow setting different tolerances across variation levels.

Thank you for considering this feature!

[bluescarni]

@alangfor this would be a good feature to have, but I am afraid it is not implementable in heyoka.

The problem is that when you set a higher tolerance (e.g., 1e-8), you will see that heyoka will take a step size similar to the lower tolerance:

>>> ta = hy.taylor_adaptive(hy.model.pendulum(), [.1, 0.])
>>> ta.step()
(<taylor_outcome.success: -4294967297>, 0.6084356584597945)
>>> ta = hy.taylor_adaptive(hy.model.pendulum(), [.1, 0.], tol=1e-8)
>>> ta.step()
(<taylor_outcome.success: -4294967297>, 0.5693063644864683)

This happens because heyoka changes the integrator order (rather than the step size) when you change tolerance. This means that in order to have per-variable tolerances, we would need to implement different Taylor series orders for different variable, but I do not think that this is feasible.

I am aware of potential issues surrounding timestepping with the variational equations, but I am a bit surprised to hear that they lead to shorter timesteps being taken. The problem I am aware of has more to do with the variational quantities exploding in magnitude (e.g., in chaotic systems), which could in principle lead to inaccurate integration due to the way heyoka interprets the tolerance setting. Would you mind if possible to show an example of the timestepping issues you are talking about?

Something else to keep in mind about the variational equations, is that you have the option to periodically rescale (by hand) the variational quantities if they get too big/too small. This is quite easy to do at the first order, but it gets more complicated for higher variational orders. Perhaps this could also help in your situation?

[alangfor]

Thanks for the quick response. No worries, it is definitely just a 'nice to have', heyoka already does a lot!

For reference, here is a script that shows the smaller step size for propagating systems with variational EOMs compared to a non-variational system at the same Taylor order. This may be driven by the chaotic dynamics near the libration points in the CR3BP, since the variational equations have terms that grow ~exponentially.

import heyoka as hy
import numpy as np

# make variables
x, y, z, xd, yd, zd = hy.make_vars('x', 'y', 'z', 'xd', 'yd', 'zd')
mu = hy.par[0]

# EOMs for planar lagrangian CR3BP
n = 1.0
r1 = hy.sqrt((x + mu)**2 + y**2)
r2 = hy.sqrt((x - (1 - mu))**2 + y**2)
Ust = n**2 * 0.5 * (x**2 + y**2) + (1 - mu) / r1 + mu / r2
xdd = hy.diff(Ust, x) + 2 * n * yd
ydd = hy.diff(Ust, y) - 2 * n * xd

# wrap EOMs for hy
state = [x, y, xd, yd]
dstate = [xd, yd, xdd, ydd]
sys = [(s, dsdt) for s, dsdt in zip(state, dstate)]
vsys = hy.var_ode_sys(sys = sys, args = state, order = 1)

mu_em = 1.21505856096e-2
ta_var = hy.taylor_adaptive(vsys,
                        pars = [mu_em],
                        tol = 1e-16,
                        compact_mode = True,
                        parallel_mode = True)

ta = hy.taylor_adaptive(sys,
                        pars = [mu_em],
                        tol = 1e-16,
                        compact_mode = True,
                        parallel_mode = True)

# lyapunov IC from JPL
x0 = np.array([6.8868283959576881e-1, 0.0, 6.6345557835782965e-1, 0.0])
T = 6.0146936600708560

# set and prop lyapunov state
ta_var.state[:vsys.n_orig_sv] = x0 
ta.state[:] = x0 
out = ta.propagate_for(T) # w/out variations
out_var = ta_var.propagate_for(T) # w/ variations

xf_var = ta_var.state[:vsys.n_orig_sv]
xf = ta.state[:]

print('\nHeyoka version', hy.__version__)

print('\nNon-variational Integration')
print('Taylor order:', ta.order)
print('Max step size:', out[2])
print('Total steps:', out[3])

print('\nVariational Integration')
print('Taylor order:', ta_var.order)
print('Max step size:', out_var[2])
print('Total steps:', out_var[3])

print('\nFinal state difference:')
print(xf_var - xf)

Heyoka version 7.0.1

Non-variational Integration
Taylor order: 20
Max step size: 0.31221466887610744
Total steps: 55

Variational Integration
Taylor order: 20
Max step size: 0.2852236973856958
Total steps: 65

Final state difference:
[1.76803017e-14 2.88657986e-15 7.91033905e-15 2.08166817e-15]

[bluescarni]

@alangfor thanks for the example and sorry for the late reply...

If you want to do some experiments with rescaling the variational variables (since you are operating at first order), the basic idea is as follows.

The variational equations for an ODE of type $\frac{dx}{dt}=f\left( x \right)$ are formulated as follows:

$$ \frac{d\delta_{x_0}}{dt} = \frac{\partial f}{\partial x}\delta_{x_0} $$

where $\delta_{x_0} = \frac{\partial x}{\partial x_0}$ are the variational variables. Trivially, you can introduce rescaled versions of $\delta_{x_0}$

$$ \Delta_{x_0} = \alpha \delta_{x_0}, $$

with $\alpha \in \mathbb{R}$, and they will obey the same differential equation:

$$ \frac{d\Delta_{x_0}}{dt} = \frac{\partial f}{\partial x}\Delta_{x_0}. $$

So in heyoka the basic idea would be that you monitor the magnitude of the variational variables and, once they get too large, you rescale their current values by a factor of your choice, and you store this factor in a list somewhere. As the integration progresses, you can keep on rescaling as needed (while always keeping track of the list of factors that you have used so far).

At the end of the integration, in order to recover the value of the original (unscaled) variational variables, you can divide their current values by the product of the factors.

Let me know if this is of any help in controlling the step size...