Non conservation of energy in ModelingToolkit with flexible bouncing ball example

BambOoxX · July 22, 2026, 4:00pm

Question

Based on the examples in the ModelingToolkit documentation I’m trying to expand a bit and model the impact behavior between two masses, one being connected to a base through a spring.

To make it simple, an example with only two masses colliding yields expected results, but as soon as I add the spring, this leads to a non-conservative behavior.

The only obvious thing I can think about is that the continuous event is only acting at the very specific gap == 0 condition while it should act on the gap <= 0 inequality condition.

What would be the proper way to account for these inequality conditions ?

Side uncertain note: this seemed to behave as expected in previous versions but I can’t remember which.

Environment

[e9467ef8] GLMakie v0.13.13
[961ee093] ModelingToolkit v11.35.2
[16a59e39] ModelingToolkitStandardLibrary v2.29.5
[1dea7af3] OrdinaryDiffEq v7.1.2

on julia 1.12.6 windows.

Two masses without spring

using ModelingToolkit
import ModelingToolkit: connect, t_nounits as t, D_nounits as D
using ModelingToolkitStandardLibrary
using GLMakie
using OrdinaryDiffEq

@named mass1 = ModelingToolkitStandardLibrary.Mechanical.TranslationalPosition.Mass(m=1.)
@named mass2 = ModelingToolkitStandardLibrary.Mechanical.TranslationalPosition.Mass(m=0.5)
@named fix = ModelingToolkitStandardLibrary.Mechanical.TranslationalPosition.Fixed(s_0=0.)
@named spring = ModelingToolkitStandardLibrary.Mechanical.TranslationalPosition.Spring(k=10., l=1.)
@variables gap(t) p(t) E(t) l(t)
@parameters e_newton = 1.


@mtkcompile twomass = System(
    [
        gap ~ mass2.s - mass1.s,
        p ~ mass1.v * mass1.m + mass2.v * mass2.m,
        E ~ 1/2 * mass1.m * mass1.v^2 + 1/2 * mass2.m * mass2.v^2,
    ],
    t;
    systems=[mass1, mass2],
    continuous_events=[gap ~ 0] => [
        mass1.v ~ ((mass1.m - 1 * mass2.m) * Pre(mass1.v) + (1 + 1) * mass2.m * Pre(mass2.v)) / (mass1.m + mass2.m)
        mass2.v ~ ((mass2.m - 1 * mass1.m) * Pre(mass2.v) + (1 + 1) * mass1.m * Pre(mass1.v)) / (mass1.m + mass2.m)
    ]
)

u0 = Dict{Num,Any}([
    mass1.s => 1.,
    mass1.v => 0.,
    mass2.s => 2.,
    mass2.v => -1.,
]
)

prob = ODEProblem(twomass, u0, (0.0, 10.0), fully_determined=true)

sol = solve(prob, Tsit5())
begin
    f = Figure()
    ax = Axis(f[1, 1])
    plot!(sol, idxs=[twomass.mass1.s, twomass.mass2.s, twomass.gap])
    axislegend(ax)
    ax = Axis(f[2, 1])
    plot!(sol, idxs=[twomass.mass1.v, twomass.mass2.v])
    axislegend(ax)
    ax = Axis(f[3, 1])
    plot!(sol, idxs=[twomass.E])
    axislegend(ax)
    f
end

which gives

Two masses with spring

@mtkcompile twomassspring = System(
    [
        gap ~ mass2.s - mass1.s,
        p ~ mass1.v * mass1.m + mass2.v * mass2.m,
        l ~ spring.flange_a.s - spring.flange_b.s,
        E ~ 1/2 * mass1.m * mass1.v^2 + 1/2 * mass2.m * mass2.v^2 + 1/2 * spring.k * (l - spring.l)^2,
        connect(fix.flange, spring.flange_b),
        connect(mass1.flange, spring.flange_a),
    ],
    t;
    systems=[mass1, mass2, fix, spring],
    continuous_events=[gap ~ 0] => [
        mass1.v ~ ((mass1.m - 1 * mass2.m) * Pre(mass1.v) + (1 + 1) * mass2.m * Pre(mass2.v)) / (mass1.m + mass2.m)
        mass2.v ~ ((mass2.m - 1 * mass1.m) * Pre(mass2.v) + (1 + 1) * mass1.m * Pre(mass1.v)) / (mass1.m + mass2.m)
    ]
)


u0 = Dict{Num,Any}([
    mass1.s => 1.,
    mass1.v => 0.,
    mass2.s => 2.,
    mass2.v => -1.,
]
)

prob = ODEProblem(twomassspring, u0, (0.0, 10.0), fully_determined=true)

sol = solve(prob, Tsit5())
begin
    f = Figure()
    ax = Axis(f[1, 1])
    plot!(sol, idxs=[twomassspring.mass1.s, twomassspring.mass2.s, twomassspring.gap, twomassspring.l])
    axislegend(ax)
    ax = Axis(f[2, 1])
    plot!(sol, idxs=[twomassspring.mass1.v, twomassspring.mass2.v])
    axislegend(ax)
    ax = Axis(f[3, 1])
    plot!(sol, idxs=[twomassspring.E])
    axislegend(ax)
    f
end

which yields

Obviously, the energy conservation is not respected, however events are properly catched.

zdenek_hurak · July 22, 2026, 10:42pm

First, note that there must be some issue in how you plot the state variables even in the first case, because the velocities of both balls do not change after the collision. But I am not much into MTK to be able to identify the issue.

Second, the fact that general RK methods such as Tsit5() do not conserve energy is not suprrising, is it? These methods have not been designed with that goal in mind. If conservation of energy is a priority for you, have a look at some symplectic solvers: Dynamical, Hamiltonian, and 2nd Order ODE Solvers · DifferentialEquations.jl (perhaps give it a try with McAte4()). But you will have to reformulate your model using the partitioned format as described at the beginning of the page Dynamical, Hamiltonian and 2nd Order ODE Problems · DifferentialEquations.jl.

BambOoxX · July 23, 2026, 9:52am

First, regarding the plots, the plot call signatures in the example are rigorously the most basic way to plot the solutions in ModelingToolkit, so I dont’t think there can be any issue with that, although something bad could internally happen.
Then, while I agree that RK methods are not generally energy-conservative, and that most integrators are not made to be energy conservative while simulating unilateral contact problems (see e.g. Making sure you're not a bot!), the behavior here is far beyond any such considerations. The energy jumps here are just plain wrong, as well as the satisfaction of unilateral contact conditions, and not integrator related.
Moreover, there is absolutely no requirement to simulate such systems as second order ones as there are plenty of first order algorithms (see Moreau-Jean or Schatzman-Paoli schemes).
Last but not least, here the contact is handled through events, so the integrator itself is not specialized for contact, rather it stops on the event, then restarts the solution with new initial conditions at the event, thanks to the impact law and momentum conservation law.

I investigated the issue by switching from ModelingToolkitStandardLibrary.Mechanical.TranslationalPosition to ModelingToolkitStandardLibrary.Mechanical.Translational components.

This results in the following two mass simulation

using ModelingToolkit
import ModelingToolkit: connect, t_nounits as t, D_nounits as D
using ModelingToolkitStandardLibrary
using GLMakie
using OrdinaryDiffEq, OrdinaryDiffEqRosenbrock

@named mass1 = ModelingToolkitStandardLibrary.Mechanical.Translational.Mass(m=1.)
@named mass2 = ModelingToolkitStandardLibrary.Mechanical.Translational.Mass(m=0.5)
@named fix = ModelingToolkitStandardLibrary.Mechanical.Translational.Fixed()
@named spring = ModelingToolkitStandardLibrary.Mechanical.Translational.Spring(Val(:relative),k=10.,)
@variables gap(t) p(t) E(t) l(t)
@parameters e_newton = 1.

@mtkcompile twomass = System(
    [
        gap ~ mass2.s - mass1.s,
        p ~ mass1.v * mass1.m + mass2.v * mass2.m,
        E ~ 1/2 * mass1.m * mass1.v^2 + 1/2 * mass2.m * mass2.v^2,
    ],
    t;
    systems=[mass1, mass2],
    continuous_events=[gap ~ 0] => [
        mass1.v ~ ((mass1.m - 1 * mass2.m) * Pre(mass1.v) + (1 + 1) * mass2.m * Pre(mass2.v)) / (mass1.m + mass2.m)
        mass2.v ~ ((mass2.m - 1 * mass1.m) * Pre(mass2.v) + (1 + 1) * mass1.m * Pre(mass1.v)) / (mass1.m + mass2.m)
    ]
)


u0 = Dict{Num,Any}([
    mass1.s => 1.,
    mass1.v => 0.,
    mass2.s => 2.,
    mass2.v => -1.,
]
)

prob = ODEProblem(twomass, u0, (0.0, 10.0), fully_determined=true)

sol = solve(prob, Tsit5(), saveat = 0:0.1:10.)
begin
    f = Figure()
    ax = Axis(f[1, 1])
    plot!(sol, idxs=[twomass.mass1.s, twomass.mass2.s, twomass.gap])
    axislegend(ax)
    ax = Axis(f[2, 1])
    plot!(sol, idxs=[twomass.mass1.v, twomass.mass2.v])
    axislegend(ax)
    ax = Axis(f[3, 1])
    plot!(sol, idxs=[twomass.E])
    axislegend(ax)
    f
end

Note that the velocity plot is correct in this case.

For the two mass + spring case I get this now


@mtkcompile twomassspring = System(
    [
        gap ~ mass2.s - mass1.s,
        p ~ mass1.v * mass1.m + mass2.v * mass2.m,
        l ~ spring.delta_s,
        E ~ 1/2 * mass1.m * mass1.v^2 + 1/2 * mass2.m * mass2.v^2 + 1/2 * spring.k * (l)^2,
        connect(fix.flange, spring.flange_b),
        connect(mass1.flange, spring.flange_a),
    ],
    t;
    systems=[mass1, mass2, fix, spring],
    continuous_events=[gap ~ 0] => [
        mass1.v ~ ((mass1.m - 1 * mass2.m) * Pre(mass1.v) + (1 + 1) * mass2.m * Pre(mass2.v)) / (mass1.m + mass2.m)
        mass2.v ~ ((mass2.m - 1 * mass1.m) * Pre(mass2.v) + (1 + 1) * mass1.m * Pre(mass1.v)) / (mass1.m + mass2.m)
    ]
)

u0 = Dict{Num,Any}([
    mass1.s => 1.,
    mass1.v => 0.,
    mass2.s => 2.,
    mass2.v => -1.,
]
)

prob = ODEProblem(twomassspring, u0, (0.0, 10.0))

sol = solve(prob, Tsit5())
begin
    f = Figure()
    ax = Axis(f[1, 1])
    plot!(sol, idxs=[twomassspring.mass1.s, twomassspring.mass2.s, twomassspring.gap, twomassspring.l])
    axislegend(ax)
    ax = Axis(f[2, 1])
    plot!(sol, idxs=[twomassspring.mass1.v, twomassspring.mass2.v])
    axislegend(ax)
    ax = Axis(f[3, 1])
    plot!(sol, idxs=[twomassspring.E])
    axislegend(ax)
    f
end

Which show as you said that the scheme is not energy conservative, but really not in the same proportions as before.

Increasing the mass 2 (mass2.m = 5.) yields multiple impacts

And changing the Newton restitution coefficient to 0 (fully plastic impact) however yields wrong results again…

In itself it’s normal here that energy dissipates at impact, but incorrect that the contact conditions are not propagated…

Given all this, I will open one issue on ModelingToolkit.

ChrisRackauckas · July 23, 2026, 11:55am

Yeah, we probably need to just delete Mechanical.TranslationalPosition

BambOoxX · July 23, 2026, 12:05pm

Well I personally find it more convenient to handle proper positions but that is debatable.
However, if it behaves incorrectly, it is better to remove it indeed.
Did you expect the unilateral contact fail with the e=0 restitution coefficient at the very end due to maintained contact ?

Jake · July 23, 2026, 12:09pm

Slightly derailing the conversation but I am curious how Gugan’s approach using the Hertzian contact equations compare, if you have used them as well?

BambOoxX · July 23, 2026, 12:16pm

Nope, I have not used this.
I just played with this toy example, to test for a more complex case.
The advantage with such simple systems is that nothing is hidden under the rug of complexity !

ChrisRackauckas · July 23, 2026, 12:38pm

It’s a property of the approach. It’s fundamental that approach can be lossy in the definition of its physics, so then it can drift with numerical error.

BambOoxX · July 23, 2026, 12:51pm

Maybe I was not very clear, but to put it in a simple manner, when I parametrize things, I like to be able to properly say where some mass is positioned, or what length does the spring have, or if two points are co-located. It feels much more natural than handling that at the velocity level.
Then for impacts, it surely has advantages to directly modify velocities.
And internally, at the solver level everything is possible ^^.

BambOoxX · August 3, 2026, 8:54am

Thanks to @cryptic.ax I got clarifications for this problem in an issue Bouncing ball like example does not respect unilateral contact conditions · Issue #4822 · SciML/ModelingToolkit.jl · GitHub.

To summarize : continuous_events cannot handle maintained conditions such as the ones that may appear during prolonged contacts. For these, one needs to use discrete_events and pay attention to the solver/timestep size. For simple “impacts” continuous_events work fine.

Topic		Replies	Views
Help with ModelingToolkit Modelling & Simulations modelingtoolkit	2	518	December 1, 2021
Surprising error on example code from ModelingToolkit documentation! Modelling & Simulations modelingtoolkit	2	439	May 13, 2022
2D arrays with ModelingToolkit Modelling & Simulations question , package , modelingtoolkit	22	1052	December 19, 2023
Understanding ModelingToolkit connectors and components General Usage question , modelingtoolkit	9	1279	December 26, 2023
ModelingToolkit compiler does not simplify as expected Modelling & Simulations question , modelingtoolkit	9	259	January 12, 2026

Non conservation of energy in ModelingToolkit with flexible bouncing ball example

Question

Environment

Two masses without spring

Two masses with spring

Related topics