WhenEvent and Energy conservation

Question

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 + ξ^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"]
 ]

Answer 1

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[
  {grad, normal}, grad = Grad[surface, vars];
  normal = Normalize[grad //. Thread[Rule[vars, pos]]];
  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"]]

Answer 2

• 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]] /. 
    Line[a_] -> {Arrowheads[{{.01, 
         1, {Graphics[{Red, Opacity[1], Disk[]}], 0}}}], Arrow[a]}, 
   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}, 
  ControlPlacement -> Bottom, AnimationRate -> 5]