# How to reflect from an ellipse while solving a system of differential equations?

While solving a system of differential equations to plot trajectories, we can use NDSolve and WhenEvent like shown here.

I am trying to solve this for a system of equations (Hamiltons Equations) for a free particle, so that the reflection is from a surface of ellipse (like elliptical pool table).

Where r and th are r and theta coordinates.

In the r, theta cordinates, the equation of ellipse reads as:

The initial conditions I have given (in the variable of time i.e. r(t), theta(t), etc) are r(t=0)= 2.29 (which is r(theta=Pi)), theta(t=0) = Pi, pr(t=0)=1, pth(t=0) = 0. Where pr and pth are radial and angular momentum given by Hamilton's Equations

This is how I am implementing it, but the reflection is occurring from a circle and not an ellipse. Please help.

The code is:

b=1;
eps=0.9;
m=1;
R=b/(Sqrt[1-(eps*Cos[th])^2]);
H = pr[t]^2/(2*m) + pth[t]^2/(2*m*r[t]^2);
xxx = NDSolve[{r'[t] == D[H, pr[t]],th'[t] == D[H, pth[t]], pr'[t] == -D[H, r[t]], pth'[t] == -D[H, th[t]],r[0] == 2.29, th[0] == Pi, pr[0] == 1, pth[0] == 0,WhenEvent[r[t]^2 == (b/Sqrt[1 - (eps*Cos[th[t]])^2])^2, pr[t] -> -pr[t]]}, {r[t], th[t], pr[t], pth[t]},{t, 0,100}]
Show[PolarPlot[{xxx[[1]][[1]][[2]]}, {t, 0, 100}],PolarPlot[b/(Sqrt[1 - (eps*Cos[th])^2]), {th, 0, 2*Pi},PlotStyle -> Red]]


The output I get is a very pretty but wrong image:

The trajectories must lie inside the ellipse (red), because the reflection condition given by WhenEvent defines the boundary of the circle.

What am I doing wrong?

• I guess this would be simpler in Cartesian coordinates... – Henrik Schumacher Jun 18 '18 at 9:38
• @HenrikSchumacher May be, but the next extension of the problem which I wish to solve will be easier to be solved in Polar Coordinates. I wish to stick to this coordinate system! – Nishchal Dwivedi Jun 19 '18 at 9:23
• Then you have to express the normal of the boundary ellipsoid and apply a reflection of the speed vector with respect to this normal in the WhenEvent. – Henrik Schumacher Jun 19 '18 at 9:29

So, the problem here is the use of PolarPlot. Once we use PolarPlot, the function takes the parameter 't' as the theta variable and hence generates faulty graph.
Here the interpolating functions generated are of $r(t)$ and $\theta(t)$. We need to do a ParametricPlot of $r(t)cos(\theta(t))$ vs $r(t)sin(\theta(t))$ and it will give the correct graph.
Show[ParametricPlot[{xxx[[1]][[1]][[2]]*Cos[xxx[[1]][[2]][[2]]], xxx[[1]][[1]][[2]]*Sin[xxx[[1]][[2]][[2]]]}, {t, 0, 100}], PolarPlot[b/Sqrt[1 - (eps*Cos[th])^2], {th, 0, 2*Pi}],PlotStyle -> Red], PlotRange -> All]

• Yes, PolarPlot specifically is looking to plot the radius as a function of $\theta$, not a parametric curve. – KraZug Jun 19 '18 at 10:49