I'm trying to simulate 2 balls with the same mass and diameter bouncing one on top of another under gravity, see the illustration below (not ideal, but this is the best result I've got so far, the numbers are the time in seconds, the height of the box is 2 meters):
However, I'm not interested in the pretty animation, but rather in the variation of the distance with time. My main motivation - I want to see how the period of this function, if it's indeed periodic, depends on the initial conditions.
For the initial values of $h_2=2$ and $h_1=1.5$, I obtained a periodic dependence with my first method:
The first plot is for the positions $y_1$ and $y_2$, while the second plot is for the distance $y_2-y_1$.
This looks promising, however the method is approximate: I'm using a constant finite time step $dt$ and check the distances between the blue ball and the ground as well as between the blue and red balls. To get the correct simulation, I had to take a very small $dt$, as you can see at the top of the image.
I have tried a better way (in my opinion), of projecting the collision times $t_a$ (between the blue ball and the ground) and $t_b$ (between the balls) and then comparing them and changing the parameters accordingly. But for the same initial conditions, the results break down as shown below, compared to the first method:
My code for the first method is:
g = 981/10; (*gravitational acceleration in meters per second squared*)
dt = 1/20000; (*the time step*)
h10 = 15/10; (*the initial heights in meters*)
h20 = 2;
h1 = h10;
h2 = h20;
u1 = 0; (*the initial velocities*)
u2 = 0;
Tm = 20; (*the time limit*)
Nm = Floor[Tm/dt]; (*number of iterations in the loop*)
d = 5/100; (*diameter of the balls, in meters*)
Y1 = Y2 = S = Table[{1, 1}, {j, 1, Nm}]; (*lists for keeping results*)
t1 = 0; (*initial times for both balls, since for each collision they need to be updated*)
t2 = 0;
t = 0; (*the time variable*)
j = 1; (*the counter*)
y1 = h1; (*initial positions and velocities, as variables*)
y2 = h2;
v1 = u1;
v2 = u2;
While[j <= Nm, (*the loop condition*)
t = t + dt; (*time propagation*)
y1 = h1 + u1 (t - t1) - 1/2 g (t - t1)^2; (*updating positions and velocities*)
y2 = h2 + u2 (t - t2) - 1/2 g (t - t2)^2;
v1 = u1 - g (t - t1);
v2 = u2 - g (t - t2);
If[y1 <= d/2, u1 = -v1; t1 = t; h1 = y1]; (*check for collision of the blue ball with the ground, the velocity is reflected, the equation parameters are updated*)
If[y2 - y1 <= d, u1 = v2; u2 = v1;
t1 = t; t2 = t; h1 = y1; h2 = y2]; (*check for collision of the two balls, the velocities are exchanged, the equation parameters are updated*)
Y1[[j]] = {t, y1}; (*writing down the results*)
Y2[[j]] = {t, y2};
S[[j]] = {t, y2 - y1};
j++];
ListPlot[{Y1, Y2}, PlotRange -> Full, ImageSize -> Full,
AspectRatio -> 1/10] (*plotting the positions*)
ListPlot[S, PlotRange -> Full, Joined -> True, ImageSize -> Full,
AspectRatio -> 1/10] (*plotting the distance between balls*)
And for the second method:
$MaxExtraPrecision = 10000;
g = 981/10;
d = 5/100;
h10 = 15/10;
h20 = 2;
h1 = h10;
h2 = h20;
u1 = 0;
u2 = 0;
t1 = 0;
t2 = 0;
Tm = 20;
T1 = {t1}; (*initial results lists, for building the plots after the computation*)
T2 = {t2};
H1 = {h1};
H2 = {h2};
U1 = {u1};
U2 = {u2};
While[t1 < Tm, (*loop condition, now for the time*)
ta = N[t1 + (u1 + Sqrt[u1^2 + 2 (h1 - d/2) g])/g, 10000]; (*the time of the collision between the blue ball and the ground*)
tb = If[u1 - u2 + g (t1 - t2) != 0,
N[(1/2 g (t1^2 - t2^2) + u1 t1 - u2 t2 - (h1 - h2) - d)/(
u1 - u2 + g (t1 - t2)), 10000], 2 Tm]; (*the time of the collision between the balls*)
If[ta <= tb || tb <= t2, (*if the collision with the ground comes first or the collision between the balls is impossible (i.e. it happened in the past and they are flying away), then we are updating only the first ball accordingly *)
h1 = d/2;
u1 = -u1 + g (ta - t1);
t1 = ta;
T1 = Append[T1, N[t1, 10]];
H1 = Append[H1, N[h1, 10]];
U1 = Append[U1, N[u1, 10]],
h1 = h1 + u1 (tb - t1) - 1/2 g (tb - t1)^2; (*else, the collision between the balls happens first and we are updating both of them accordingly*)
h2 = h1 + d;
v1 = u2 - g (tb - t2); (*this dummy variable is needed because u1 will be used to update u2*)
u2 = u1 - g (tb - t1);
u1 = v1;
t1 = t2 = tb;
T1 = Append[T1, N[t1, 10]];
H1 = Append[H1, N[h1, 10]];
U1 = Append[U1, N[u1, 10]];
T2 = Append[T2, N[t2, 10]];
H2 = Append[H2, N[h2, 10]];
U2 = Append[U2, N[u2, 10]]]];
y[t_, h_, u_, t0_] := h + u (t - t0) - 1/2 g (t - t0)^2; (*the general position function*)
L1 = Length[T1]; (*number of collisions for the blue ball*)
L2 = Length[T2]; (*number of collisions for the red ball*)
Y1[t_] :=
Piecewise[
Table[{y[t, H1[[k]], U1[[k]], T1[[k]]],
T1[[k]] <= t < T1[[k + 1]]}, {k, 1, L1 - 1}]];
Y2[t_] :=
Piecewise[
Table[{y[t, H2[[k]], U2[[k]], T2[[k]]],
T2[[k]] <= t < T2[[k + 1]]}, {k, 1, L2 - 1}]]; (*positions over time as piecewise functions*)
Plot[Y2[t] - Y1[t], {t, 0, Tm}, AspectRatio -> 1/10,
ImageSize -> Full, PlotRange -> {0, h20}]
I thought the problem was numerical approximation to $t_a$ and $t_b$, however after greatly increasing precision, I have not seen any change in results.
Why does the second method break down? Shouldn't it be more exact than the first one? Where's my error?
To clarify: I'm not that interested in other methods, such as solving differential equations, I think this problem doesn't need them, however if there's some improvement possible for my methods while keeping the general idea, I will be very grateful.
Update
I have restarted Mathematica and the second code works just fine (even though it doesn't give the perfect periodic dependence of the first case).
So my new question: what's the problem with the second code, why is it unreliable? I'd really love to hear any thoughts on this.
Here are the results for $h_1=1$ and $h_2=2$, they are quite different, after enough collisions, so another question would be: which method do I trust more?
Update 2 - Important!
I have $g$ 10 times larger than normal above, so the time dependence needs to be scaled down.
Also, the second algorithm had been breaking down for $d=0$ in periodic cases, because when both balls were reaching the ground at the same time, their velocities needed to change signs before exchanging them, thus the second ball had been going below the ground. I fixed the problem in this part of the algorithm:
If[ta == tb,
h2 = h2 + u2 (tb - t2) -
1/2 g (tb -
t2)^2;(*first, an exotic case of both balls reaching the \
ground at the same time, only possible for d=0, I had to make it \
though to avoid errors*)
h1 = h1 + u1 (tb - t1) - 1/2 g (tb - t1)^2;
v1 = -u2 + g (tb - t2);
u2 = -u1 + g (tb - t1);
u1 = v1;
t1 = tb;
t2 = tb;
T1 = Append[T1, N[t1, 10]];
H1 = Append[H1, N[h1, 10]];
U1 = Append[U1, N[u1, 10]];
T2 = Append[T2, N[t2, 10]];
H2 = Append[H2, N[h2, 10]];
U2 = Append[U2, N[u2, 10]]]; If[
ta < tb ||
tb <= t2,(*if the collision with the ground comes first or the \
collision between the balls is impossible (i.e.it happened in the \
past and they are flying away),then we are updating only the first \
ball accordingly*)
h1 = d/2;
u1 = -u1 + g (ta - t1);
t1 = ta;
T1 = Append[T1, N[t1, 10]];
H1 = Append[H1, N[h1, 10]];
U1 = Append[U1, N[u1, 10]],
If[ta > tb && tb > t2,
h2 = h2 + u2 (tb - t2) - 1/2 g (tb - t2)^2;(*else,
the collision between the balls happens first and we are updating \
both of them accordingly*)
h1 = h1 + u1 (tb - t1) - 1/2 g (tb - t1)^2;
v1 = u2 -
g (tb - t2);(*this dummy variable is needed because u1 will be \
used to update u2*)
u2 = u1 - g (tb - t1);
u1 = v1;
t1 = tb;
t2 = tb;
T1 = Append[T1, N[t1, 10]];
H1 = Append[H1, N[h1, 10]];
U1 = Append[U1, N[u1, 10]];
T2 = Append[T2, N[t2, 10]];
H2 = Append[H2, N[h2, 10]];
U2 = Append[U2, N[u2, 10]]]]
Nand letting$MaxExtraPrecision = Infinitythen probably waiting for a few hours you can get exact results. $\endgroup$