q = {{x[t]}, {y[t]}};
dq = D[q, t];
ddq = D[dq, t];
dx = D[x, t];
dy = D[y, t];
ddx = D[dx, t];
ddy = D[dy, t];
KE = 0.5 m (x'[t]^2 + y'[t]^2);
PE = m*g*y[t];
L = KE - PE;
(*The constraint*)
\[Phi] = y;
d\[Phi] = D[\[Phi], q];
Eq = (Thread[D[D[L, dq\[Transpose]], t] - D[L, q\[Transpose]] == 0]);
ELtemp = Solve[Eq[[1]] && Eq[[2]], Flatten[ddq]];
EL = {x''[t] == ELtemp[[1, 1, 2]], y''[t] == ELtemp[[1, 2, 2]]};
p = D[L, dq\[Transpose]];
H = p.dq - L;
m = 1;
g = 9.8;
ICnt = 5;
xzero = 0; xdzero = 3; yzero = 10; ydzero = 0; Tzero = 0;
x1 = {}; y1 = {}; 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[{x'[Tplus]}, {y'[Tplus]}, {\[Lambda]}]];
For[i = 0, i < ICnt, i++,
InitCon = {\[Theta][Tzero] == xzero &&
yzero, \[Theta]'[Tzero] == -xdzero && -ydzero};
sol = NDSolve[Join[EL, Initial], {x, y}, {t, Tzero, 30},
Method -> {"EventLocator", "Event" -> y,
"EventAction" :>
Throw[{Tend = t, yend = y[t], ydotend = y'[t]},
"StopIntegration"]}];
xsol = x[t] /. sol[[1]]; ysol = y[t] /. sol[[1]];
Hsol = H /. sol[[1]];
x1 = Append[x1, xsol]; y1 = Append[y1, ysol];
Time = Append[Time, Tend];
Hset = Append[Hset, {Hsol, Tzero < t <= Tend}];
P1 = Append[P1, {xsol && ysol, Tzero < t <= Tend}];
yzero = yend; ydzero = ydotend;
Tzero = Tend
];
\[Theta]finalsol = Piecewise[P1];
Hfinalsol = Piecewise[Hset];
(*InitCon={x'[0]\[Equal]3, x[0]==0, y[0]==10, y'[0]==0};
sol=NDSolve[Join[EL,InitCon],{x,y},{t,0,30}] ;
xsol=x[t]/.sol[[1]];
ysol=y[t]/.sol[[1]];*)
(*PLOTS*)
Plot[xsol, {t, 0, 2}, AxesLabel -> {t, x}]
Plot[ysol, {t, 0, 2}, AxesLabel -> {t, y}]
ParametricPlot[{xsol, ysol}, {t, 0, 2}, AxesLabel -> {x, y}]
(*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];
(*ANIMATION*)
X[\[Tau]_] := xsol /. \[Theta]finalsol /. t -> \[Tau];
Y[\[Tau]_] := ysol /. \[Theta]finalsol /. t -> \[Tau];
Animate[Graphics[{PointSize[0.05], Point[{ X[\[Tau]], Y[\[Tau]]}],
Line[{{-10, 0}, {10, 0}}], Line[{{0, -10}, {0, 10}}]}], {\[Tau], 0,
2}]
I tried this to get a ball to bounce to 5 times and then stop. But it gives me certain errors like equations are not a quantified system of equalities and some replacement rules are not valid.
Any suggestions on how to proceed?
Thanks
Varun
Attachments: