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?
WhenEvent
. $\endgroup$