# 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}] :=   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"]  ] ``