Regularization

Regularization is an optional, advanced feature. The core integrator runs unregularized by default; opt in by wrapping the base algorithm in a RegularizedIntegrator and passing it to solve.

When to use it

The Weber potential contains a Coulomb singularity at zero separation. For most simulations with well-separated particles this is not an issue. Enable regularization when:

• Two or more particles are expected to have close encounters (separation approaching zero).
• You observe energy drift or NaN values near close passages without it.

Regularization handles the Coulomb/Kepler singularity only. Weber's velocity-dependent correction is not regularized — the LC and KS transforms apply to the conservative part of the force.

Enabling regularization

prob = HamiltonianProblem(
    sys, tspan, q0, p0;
    masses  = [1.0, 1.0],
    charges = [1.0, -1.0],
    c       = 10.0,
    dt      = 0.01,
)

alg = RegularizedIntegrator(SymmetricProjectionIntegrator())   # opt in
sol = solve(prob, alg)

All keyword arguments below go on the RegularizedIntegrator(...) constructor itself. Its kwargs mirror RegularizationOptions (with enabled = true set implicitly).

Backends

For 3D problems, :lifted_pair automatically falls back to :adaptive_cartesian (with an optional warning controlled by warn_on_fallback).

# Explicit 3D choice
alg = RegularizedIntegrator(
    SymmetricProjectionIntegrator();
    backend = :adaptive_cartesian,
)
sol = solve(prob, alg)

Activation hysteresis

Regularization switches on when any pair separation drops below r_on and switches off once all separations in the active component exceed r_off. The thresholds are computed from the initial minimum separation by default:

r_on  = r_on_factor  × min_initial_separation   (default factor: 0.15)
r_off = r_off_factor × min_initial_separation   (default factor: 0.25)

You can override them directly:

alg = RegularizedIntegrator(
    SymmetricProjectionIntegrator();
    r_on  = 0.05,
    r_off = 0.10,
)
sol = solve(prob, alg)

Chain mode

When three or more particles form a connected close-encounter cluster, regularization falls back to chain mode (adaptive Cartesian substeps for the whole cluster). Chain mode is enabled by default; disable with RegularizedIntegrator(...; chain_enabled = false).

Collision bounce

For head-on (ℓ = 0) collisions between like-charge pairs, a reflection boundary can be applied before each macro-step. This avoids the non-regularizable ℓ ≠ 0 singularity (where particles reach r = 0 at infinite speed).

Collision bounce is a CollisionBounce callback; pass it through the callbacks kwarg of solve (or init):

sol = solve(prob, SymmetricProjectionIntegrator();
            callbacks = CollisionBounce(0.02))   # reflect at r < 0.02

As a convenience, RegularizedIntegrator also accepts a collision_bounce_radius kwarg and synthesises a matching callback automatically:

alg = RegularizedIntegrator(SymmetricProjectionIntegrator();
                            collision_bounce_radius = 0.02)
sol = solve(prob, alg)

Collision bounce works best without Levi-Civita regularization (the unregularized symplectic integrator keeps energy bounded across the bounce).

Diagnostics

Every HamiltonianSolution carries a regularization::RegularizationDiagnostics field with step counts, backend used, minimum encounter distance, and more.

sol = solve(prob, SymmetricProjectionIntegrator())
d = sol.regularization
println("Backend used: ", d.used_backend)
println("Activation count: ", d.activation_count)
println("Min encounter distance: ", d.min_encounter_distance)

Key fields and what they mean:

API reference

RegularizationOptions(; kwargs...)

RegularizationDiagnostics