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"] ] `