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$ – Henrik Schumacher Jun 18 '18 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$ – Nishchal Dwivedi Jun 19 '18 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$ – 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]

enter image description here

| improve this answer | |
  • $\begingroup$ Yes, PolarPlot specifically is looking to plot the radius as a function of $\theta$, not a parametric curve. $\endgroup$ – SPPearce Jun 19 '18 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 '18 at 20:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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