How do I simulate the impact of a ball with a wall a finite number of times (say 5)?
I could simulate the impact of a pendulum with a wall at x=0 (5 times). I used the Event Locator method for the same. How do I do it similarly for a ball hitting the floor at y=0?
q = {{\[Theta][t]}};
dq = D[q, t];
ddq = D[dq, t];
x = R*Sin[\[Theta][t]]; (*X and Y Coordinates*)
y = -R*Cos[\[Theta][t]];
dx = D[x, t];
dy = D[y, t];
KE = (1/2)*m*(dx^2 + dy^2);
PE = m*g*y;
L = KE - PE; (*Lagrangian is the difference of Kinetic and \
Potential Energies*)
(*The constraint*)
\[Phi] = x;
d\[Phi] = D[\[Phi], q];
Eq = (Thread[
D[D[L, dq\[Transpose]], t] - D[L, q\[Transpose]] ==
0]); (*Solving the Euler Lagrange Equations*)
ELtemp = Solve[Eq[[1]], Flatten[ddq]];
EL = {\[Theta]''[t] == ELtemp[[1, 1, 2]]};
(*Hamiltonian*)
p = D[L, dq\[Transpose]];
H = p.dq - L;
R = 1;
m = 1;
g = 9.81;
ICnt = 5;(*The impact is to be taken 5 times*)
\[Theta]zero =
Pi/2; \[Theta]dzero = 0; Tzero = 0;(*Setting the initial conditions*)
\
\[Theta]1 = {}; Time = {}; P1 = {}; Hset = {};
(*The calculations of the Lagrangian just before and after impact \
with the wall*)
FullSimplify[D[L, dq]*dq - L];
(*This is the condition just before impact and is constrained to hit \
the wall along x=0*)
FullSimplify[
eq1 = (Thread[(D[L, dq] /. t -> Tplus) - (D[L, dq] /.
t -> Tminus) == \[Lambda]*d\[Phi]])];
(*This is the condition just after inpact and there is no constraint*)
FullSimplify[
eq2 = (Thread[((D[L, dq]*dq - L) /.
t -> Tplus) - ((D[L, dq]*dq - L) /. t -> Tminus) == 0])];
ImpactTemp =
Solve[eq1[[1]] && eq2[[1]], Join[{\[Theta]'[Tplus]}, {\[Lambda]}]];
For[
i = 0, i < ICnt, i++,
Initial = {\[Theta][Tzero] == \[Theta]zero, \[Theta]'[
Tzero] == -\[Theta]dzero};
sol = NDSolve[Join[EL, Initial], {\[Theta]}, {t, Tzero, 10},
Method -> {"EventLocator", "Event" -> R*Sin[\[Theta][t]],
"EventAction" :>
Throw[{Tend =
t, \[Theta]end = \[Theta][t], \[Theta]dotend = \[Theta]'[
t]}, "StopIntegration"]}];
\[Theta]sol = \[Theta][t] /. sol[[1]]; Hsol = H /. sol[[1]];
\[Theta]1 = Append[\[Theta]1, \[Theta]sol];
Time = Append[Time, Tend];
Hset = Append[Hset, {Hsol, Tzero < t <= Tend}];
P1 = Append[P1, {\[Theta]sol, Tzero < t <= Tend}];
\[Theta]zero = \[Theta]end; \[Theta]dzero = \[Theta]dotend;
Tzero = Tend
];
\[Theta]finalsol = Piecewise[P1];
Hfinalsol = Piecewise[Hset];
(*Hamiltonian and its Plot*)
(*Hamiltonian and its Plot
p=D[L,Transpose[dq]];
H=p.dq-L/.sol[[1]]/.Hsol[[1]];
Plot[H,{t,0,10}]*)
(*Hamiltonian and its Plot*)
Plot[\[Theta]finalsol, {t, 0, Tend},
AxesLabel -> {t, HoldForm[\[Theta]]}]
Plot[Hfinalsol, {t, 0, Tend}, AxesLabel -> {t, HoldForm[H]}]
(*Plot[xsol,{t,0,10}, AxesLabel\[Rule]{t,x}]
Plot[ysol,{t,0,10}, AxesLabel\[Rule]{t,y}]*)
X[\[Tau]_] := R*Sin[\[Theta]finalsol] /. t -> \[Tau];
Y[\[Tau]_] := -R*Cos[\[Theta]finalsol] /. t -> \[Tau];
Animate[Graphics[{PointSize[0.02], PlotRange -> {{-3, 3}, {-3, 3}},
Line[{{-1.5, 0}, {1.5, 0}}], Line[{{0, 1.5}, {0, -1.5}}],
Line[{{0, 0}, {X[\[Tau]], Y[\[Tau]]}}],
Point[{X[\[Tau]], Y[\[Tau]]}]}], {\[Tau], 0, Tend}]
Attachments: