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).

H=\frac{pr(t)^2}{2 m}+\frac{pth(t)^2}{2 m r(t)^2}

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:

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:enter image description here

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?

  • $\begingroup$ I guess this would be simpler in Cartesian coordinates... $\endgroup$ Jun 18, 2018 at 9:38
  • $\begingroup$ @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! $\endgroup$ Jun 19, 2018 at 9:23
  • $\begingroup$ 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. $\endgroup$ Jun 19, 2018 at 9:29

1 Answer 1


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]

enter image description here

  • $\begingroup$ Yes, PolarPlot specifically is looking to plot the radius as a function of $\theta$, not a parametric curve. $\endgroup$
    – SPPearce
    Jun 19, 2018 at 10:49
  • $\begingroup$ It looks like the angles of incidence do not equal the angles of reflection in your plot. I suspect the Sin and Cos in your equation would work better for a circle than an ellipse. $\endgroup$
    – Bill Watts
    Jun 20, 2018 at 20:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.