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}) :=
Block({grad, normal},
grad = Grad(surface, vars);
normal = Normalize(grad //. Thread(Rule(vars, pos)));
-(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 + (Xi)^2 + (Zeta)^2, {x(t),
y(t)}, {Derivative(1)(x)(t),
Derivative(1)(y)(t)}, {(Xi), (Zeta)})((1)),
Derivative(1)(y)(t) ->
reflect(-1 + (Xi)^2 + (Zeta)^2, {x(t),
y(t)}, {Derivative(1)(x)(t),
Derivative(1)(y)(t)}, {(Xi), (Zeta)})((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")
)
```