# WhenEvent and Energy conservation

I'm using NDSolve to do a simple Sinai billiard simulation. But the performance isn't very good using WhenEvent to reflect the ball off the walls. Apparently the energy is not conserved. How can I use NDSolve to correctly find the solution? (Maybe there is some method it can be supplied to ensure energy conservation?)

Here is my code. The billiard table is a square with a circular pillar in the center. reflect gives the new velocity after the ball hits the pillar. However, even though I have numerically verified that the reflect function conserves energy, when NDSolve uses it in WhenEvent the energy is not conserved.

reflect[surface_, pos_, vel_, vars_ : {x, y}] :=
-(2 (normal . vel)*normal - vel)
]

sinaiSolver[initialData_, duration_ : 100] :=
Block[{eqns, sol, pillar},
{{xi, yi}, {vxi, vyi}} = initialData;
eqns = {x''[t] == 0, y''[t] == 0,
x[0] == xi, y[0] == yi,
x'[0] == vxi, y'[0] == vyi,
WhenEvent[x[t]^2 + y[t]^2 - 1 == 0,
{Derivative[1][x][t] ->
reflect[-1 + ξ^2 + ζ^2, {x[t],
y[t]}, {Derivative[1][x][t],
Derivative[1][y][t]}, {ξ, ζ}][[1]],
Derivative[1][y][t] ->
reflect[-1 + ξ^2 + ζ^2, {x[t],
y[t]}, {Derivative[1][x][t],
Derivative[1][y][t]}, {ξ, ζ}][[2]]}
],
WhenEvent[x[t] == -2, x'[t] -> -x'[t]],
WhenEvent[x[t] == 2, x'[t] -> -x'[t]],
WhenEvent[y[t] == -2, y'[t] -> -y'[t]],
WhenEvent[y[t] == 2, y'[t] -> -y'[t]]};
sol = NDSolve[eqns, {x, y}, {t, 0, duration}]
]


Here is a plot that shows energy isn't conserved. The trajectories still look fine if you plot them, except that the particle changes speed on reflection with the pillar.

sol1 = sinaiSolver[{{1, 1}, {0.345, 0.493}}, 1000];
Block[{energy},
energy = 1/2 (x'[t]^2 + y'[t]^2);
Plot[energy /. sol1, {t, 0, 100}, PlotRange -> All,
PlotTheme -> "Scientific",
FrameLabel -> "Energy vs Time"]
]


• First observation: you have a typo in your definition of energy (square is outside brackets instead on the $y'$). However, this doesn't solve the problem yet ... Commented Aug 15, 2021 at 17:49
• @Domen Oh thanks, I'll fix that and attach the new plot. The typo was incurred when writing code for the question and wasn't in the code where I originally made the observation. Commented Aug 15, 2021 at 18:42
• Well, reflect is a problem. It is a bit hard to debug because the purpose of vars is not clear and because the physics is also not clear. Furthermore, I do not understand why you are using 2nd order equations when 1st order should be sufficient. Commented Aug 15, 2021 at 18:44
• @yarchik reflect implements specular reflection: en.wikipedia.org/wiki/Specular_reflection#Vector_formulation . vars is needed so that reflect can compute the surface normal using Grad given an equation for the surface. Although maybe that can be hard-coded at the beginning of the solver to save time. I don't think it would help with the energy conservation though. Commented Aug 15, 2021 at 18:51

You have an interesting problem in your code. To debug it, we can include a Print function inside reflect:

reflect[surface_, pos_, vel_, vars_ : {x, y}] := Block[
Print[surface, "   ", pos[[1]]^2 + pos[[2]]^2, "   ", pos, vel,
vars, -(2 (normal . vel)*normal - vel)];
-(2 (normal . vel)*normal - vel)]


By observing the output of your solver, we can see that for each reflection at the central circle, the reflection function is being called twice and the final result is therefore wrong! The reason for this is that in your WhenEvent, you are calling the reflect twice, separately for each component. You should call it only once and then use the components to update the velocities.

sinaiSolver[initialData_, duration_ : 100] :=
Block[{eqns, sol, pillar}, {{xi, yi}, {vxi, vyi}} = initialData;
eqns = {
x'[t] == vx[t], y'[t] == vy[t], x[0] == xi, y[0] == yi,
vx[0] == vxi, vy[0] == vyi,
WhenEvent[
x[t]^2 + y[t]^2 - 1 == 0, {
reflection = reflect[-1 + ξ^2 + ζ^2, {x[t], y[t]}, {vx[t], vy[t]}, {ξ, ζ}],
vx[t] -> reflection[[1]], vy[t] -> reflection[[2]]
}
],
WhenEvent[x[t] == -2, vx[t] -> -vx[t]],
WhenEvent[x[t] == 2, vx[t] -> -vx[t]],
WhenEvent[y[t] == -2, vy[t] -> -vy[t]],
WhenEvent[y[t] == 2, vy[t] -> -vy[t]]};
sol = NDSolve[eqns, {x, y}, {t, 0, duration},
DiscreteVariables -> {vx, vy}]]

tf = 100;
sol1 = sinaiSolver[{{1, 1}, {0.345, 0.493}}, tf];
Show[Graphics@Disk[],
Graphics[{EdgeForm[Black], Transparent,
Rectangle[{-2, -2}, {2, 2}]}],
ParametricPlot[{x[t], y[t]} /. sol1, {t, 0, tf}, AspectRatio -> 1,
PlotRange -> {{-3, 3}, {-3, 3}}]]


Note that I have made some simplification in the code (such as the reduction from the 2nd order to the 1st order of ODEs), but overall this is not really necessary and you can leave it your way -- the key fix-up is calling the reflect only once.

This now conserves the energy up to the accuracy of the solver.

Block[{energy}, energy = 1/2 (x'[t]^2 + y'[t]^2);
Plot[Evaluate[{energy} /. sol1], {t, 0, tf}, PlotRange -> All,
PlotTheme -> "Scientific", FrameLabel -> "Energy vs Time"]]


• Oof! Yes, your solution nails it thank you. In particular I didn't know I had the power to write reflection =  inside that block, and when you try to assign the new {x'[t],y'[t]} directly to reflect[...] it gets confused about evaluation order in a manner that causes a stall. Very nice solution! Commented Aug 15, 2021 at 18:58
• We product the 5 equations to construct (x - 2) (x + 2) (y - 2) (y + 2) (x^2 + y^2 - 1)==0 as condition to simplified some codes.
reflect[vector_,
normal_] = -(vector - 2 (vector - Projection[vector, normal]));
f[x_, y_] = (x - 2) (x + 2) (y - 2) (y + 2) (x^2 + y^2 - 1);
plot = ContourPlot[f[x, y] == 0, {x, -2.5, 2.5}, {y, -2.5, 2.5},
PlotPoints -> 50];
sol = NDSolveValue[{x''[t] == 0,
y''[t] ==
0, {x[0], y[0]} == {1, 1}, {x'[0], y'[0]} == {0.345, 0.493},
WhenEvent[
f[x[t], y[t]] ==
0, {{x'[t], y'[t]} ->
reflect[{x'[t], y'[t]}, {Derivative[1, 0][f][x[t], y[t]],
Derivative[0, 1][f][x[t], y[t]]}]}]}, {x[t], y[t]}, {t, 0,
150}, MaxStepSize -> .01];
ani = Animate[
Show[plot,
ParametricPlot[sol, {t, 0, c}, PlotPoints -> 200,
PlotStyle -> Directive[Opacity[.2], Gray]] /.
PlotRange -> 2.5], {c, $MachineEpsilon, 150}, ControlPlacement -> Bottom, AnimationRate -> 1]  • Test another curves. Clear["Global*"]; reflect[vector_, normal_] = -(vector - 2 (vector - Projection[vector, normal])); funs = {f[x_, y_], g[x_, y_], h[x_, y_]} = {x^2/8^2 + y^2/6^2 - 1, x^6 - 5 x^4 y + 3 x^4 y^2 + 10 x^2 y^3 + 3 x^2 y^4 - y^5 + y^6 - 2, (x - 2)^2 + (y - 2)^2 - 1}; plot = ContourPlot[funs == 0, {x, -10, 10}, {y, -10, 10}, PlotPoints -> 50, AspectRatio -> Automatic]; sol = NDSolveValue[{x''[t] == 0, y''[t] == 0, {x[0], y[0]} == {1, 1}, {x'[0], y'[0]} == {5, 6}, WhenEvent[#[x[t], y[t]] == 0, {{Derivative[1][x][t], Derivative[1][y][t]} -> reflect[{Derivative[1][x][t], Derivative[1][y][t]}, {Derivative[1, 0][#][x[t], y[t]], Derivative[0, 1][#][x[t], y[t]]}]}] & /@ {f, g, h}}, {x[t], y[t]}, {t, 0, 80}, MaxStepSize -> 0.001]; ani = Animate[ Show[plot, ParametricPlot[sol, {t, 0, c}, PlotPoints -> 300, PlotStyle -> Directive[Opacity[.2], Gray]] /. Line[a_] -> {Arrowheads[{{.01, 1, {Graphics[{Red, Opacity[1], Disk[]}], 0}}}], Arrow[a]}, PlotRange -> All], {c,$MachineEpsilon, 80},
`