SOLVE_NIVP

I had the pleasure of peer-reviewing the research software package solve_nivp, published in the Journal of Open Source Software (JOSS), DOI: 10.21105/joss.09775, authored by David Riley and Ioannis Stefanou.

I was invited as co-reviewer alongside Ductho Le and WhenXuan by Editor Daniel S. Katz. This post is not a review log, it is an explanation of what solve_nivp actually does, since the problem it solves is genuinely interesting and under-discussed.

What is Wrong with Non-Smooth Systems?

Smooth ODEs

A classical ordinary differential equation (ODE) looks like this:

where is the state and is a smooth function. “Smooth” here means differentiable, well-behaved, and (crucially) friendly to solvers like Runge-Kutta, DOPRI, or SciPy’s solve_ivp^[1]. One takes small steps, estimates the derivative, continues by moving forward, repeat. It works beautifully when “behaves.”

The issue is that a large class of real systems does not behave in the “ideal” proposed way; leaving them seemingly ignored.

Impacts, Switching, Constraints in Nonsmooth Systems

Consider, for example, a ball bouncing on a floor. Between bounces, the dynamics are perfectly smooth:

where is height and is gravitational acceleration. But the moment hits zero, the velocity flips:

where is the velocity just before impact, is just after, and is the coefficient of restitution where (1 => perfectly elastic and 0 => dead drop). This is the Newton impact law; and it is not differentiable. The velocity field has a discontinuity at the “constraint” surface .

Other examples:

• A relay circuit where a diode switches on at a threshold voltage
• A sliding-mode controller that changes sign when the error crosses past zero
• A dry-friction joint where stick and slip have entirely different dynamics
• Financial models with hard constraints (prices must be non-negative, quantities cannot be fractional)

All of these belong to the family of nonsmooth dynamical systems. Specifically, solve_nivp targets systems expressible as differential inclusions:

where is a set-valued map, at smooth points, is a singleton and you recover the classical ODE. At discontinuities however, returns a set of feasible velocities (the Filippov convex hull at a switching surface, or a single value constrained by an inequality).

Why Classical Solvers Fail

Stiffness at a given Discontinuity

If you take a smooth solver (Runge-Kutta, Adams-Bashforth) and naively apply it to a nonsmooth system, two unfortunate results occur:

1. Event detection is now brittle! The solver doesn’t know when an impact occurs, this leads it to overshoot the constraint surface, finds itself on the wrong side, and then attempts to integrate backwards toward the event. This requires rootfinding at every step near a contact surface and gets extremely expensive.

2. Regularisation makes everything stiff. One of the most prevalent workarounds is to replace the hard constraint with a smooth approximation, which is a steep sigmoid instead of a step (This can be thought of as taxing the use of an illegal item, instead of banning it.) This introduces a large Lipschitz constant, which forces explicit solvers to use tiny timesteps or blow up in terms of memory, time, etc. You end up taking thousands of steps to resolve what is fundamentally a single instance.

It is worse for a Filippov system (switching at a surface): classical solvers exhibit chattering near the sliding surface, flipping between modes at each step, which leads to never having proper convergence.

Interesting Approach: Encode the Nonsmoothness

The cleaner solution, which solve_nivp implements, is to build the nonsmooth rules directly into the time-stepping scheme from the initial, rather than trying to smooth them away or detect them as events opposed to the ideal system.

The Time-Stepping Scheme

Implicit Euler as the Foundation

The implicit Euler method updates state via:

Unlike its explicit counterpart , this requires solving a (possibly nonlinear) system at each step. However the stability desired is derived fundamentally from it’s unconditional stableness for stiff problems.

Augmenting

Now, suppose the system has a contact constraint: state must satisfy (e.g., position above floor), and there exists a contact force that only acts when the constraint is active.

The complementarity condition encodes this:

The last conditions says that either the constraint is inactive i.e. and the force is zero, or the constraint is active i.e. and the force is whatever is needed to prevent penetration of system. These 3 conditions altogether form a Linear Complementarity Problem (LCP) when the dynamics are linear, or a Nonlinear Complementarity Problem (NCP) otherwise.

The augmented implicit step at each time is to find any given satisfying

where is the constraint Jacobian. This is a Mixed Complementarity Problem (MCP), and this is what solve_nivp solves at each step instead of tradition.

Handling Impacts

For our previous bouncing ball from earlier, the impact law is applied exactly at the transition step. When the solver detects that the constraint surface is crossed (i.e. in the unconstrained prediction), it applies the restitution law:

and then continues from the new post-impact state. There is surprisingly no extra-smoothing present here! =D

What solve_nivp Provides

The API

The library’s entry point is solve_ivp_ns. One specifies a given constraint type via a projection argument (e.g. 'unilateral' for a one-sided contact, 'identity' for unconstrained), and choose a nonlinear solver ('VI' for variational-inequality fixed-point, or 'semismooth_newton'). Our lovely bouncing ball now looks like this:

1
2
3
4
5
6
7
8
9
10
11
12
13
14

import numpy as np
import solve_nivp

t, y, h, fk, info = solve_nivp.solve_ivp_ns(
    fun=rhs,                          # right-hand side f(t, x)
    t_span=(0.0, 5.0),
    y0=np.array([0.0, 1.0]),          # [velocity, height]
    method='trapezoidal',
    projection='unilateral',          # encodes q >= 0 contact constraint
    solver='VI',
    adaptive=True,
    atol=1e-6,
    rtol=1e-3,
)

The function returns a 5-tuple: (t, y, h, fk, info); that being: time array, state trajectory, step sizes, constraint forces, and a solver info dictionary respectively. The solver is then setting up and solving the complementarity system with the given arguments at each step.

Which Systems?

solve_nivp covers three main categories including the ones mentioned earlier to this:

• Unilateral constraints with impacts, the bouncing ball, rigid-body contact, granular media
• Switching systems (Filippov), sliding-mode controllers, relay circuits, piecewise-linear dynamics
• DAEs with inequality constraints, mechanical systems with joints that can separate or lock

For each category there are worked examples in the repository^[2].

Documentation and Tests

The solve_nivp package is clearly shipped with the following on their github releases:

• A full Sphinx documentation site with mathematical background
• Jupyter notebook tutorials for each system class
• A test suite covering the main integrator cases and edge conditions (zero-velocity impact, constraint grazing, etc.

Conclusion

Nonsmooth systems are not edge cases in engineering, they are everywhere and they must be responded to and considered for. Any system with contact, switching, or hard constraints is nonsmooth, and the classical ODE toolkit cannot handle them out-of-the-box. solve_nivp brings the complementarity approach to python with a clean API and proper documentation.

What stands out most to me is that most numerical software in this space is written as Fortran-wrapped Matlab, roughly readable by no one (including me!). This was rather refreshing python.

Thank you to David Riley and Ioannis Stefanou for the work on a genuinely useful package, and to Daniel S. Katz for editing the submission and inviting me.

The published article: DOI: 10.21105/joss.09775