# Two bouncing balls in 1 dimension, issues with two different methods?

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]]]]  • Could you add some comments to your code to help people understand what you're doing? – Chris K Oct 29 '17 at 19:15 • @ChrisK, I will as soon as possible, sorry for the oversight – Yuriy S Oct 29 '17 at 19:21 • By removing the N and letting $MaxExtraPrecision = Infinity then probably waiting for a few hours you can get exact results. – b3m2a1 Oct 29 '17 at 20:03
• @b3m2a1, thank you for the advice, I think I'll just increase precision until the difference between the results is small enough for me – Yuriy S Oct 29 '17 at 20:11
• – b3m2a1 Oct 29 '17 at 21:24

Your first block of code never makes an approximation. Mathematica will do exact arithmetic when passed exact parameters.

Your second block of code makes many approximations. By removing the N calls you can recover the exact results from the first code, but it takes way longer to run (because it needs massive amounts of memory to store all the intermediate symbolic results).

On the other hand, I'm not entirely clear as to why you want exact results.

I had some fun and wrote a quick little compiled, generalized version of your code which let me speed things up some.

Here's what you get out of that:

AbsoluteTiming[
trajData =
bounceBalls[
20,
{
<|"Height" -> 157/100|>,
<|"Height" -> 2|>
},
"TimeStep" -> 1./30000,
"Gravity" -> 10*OptionValue[bounceBalls, "Gravity"]
];
]

{12.4596, Null}

traj = Thread[ {#["Time"], #["Heights"]}] & /@ trajData // Transpose;
ListLinePlot[traj,
PlotStyle -> AbsoluteThickness[1],
AspectRatio -> 2/10
]


pdiff = {#["Time"], Abs[Subtract @@ #["Heights"]]} & /@ trajData;
ListLinePlot[pdiff,
PlotStyle -> AbsoluteThickness[1],
AspectRatio -> 2/10
]


The difference between this and the exact result is small:

Y1[[All, 2]] - traj[[1, All, 2]] // Mean

-0.000439656

Y2[[All, 2]] - traj[[2, All, 2]] // Mean

0.0000223441


Here's the code for it:

bounceBallsCore =
Compile[{
{ballInit, _Real, 2},
{g, _Real},
{t0, _Real},
{timeStep, _Real},
{steps, _Integer}
},
(* balls are specified by {y, v, t, r, e} *)

Module[{ balls, t = t0, dt},
(* turn this into {y, v, y0, v0, t0, r, e}, where y0, v0,
and t0 are fed into the eq. of motion *)
balls =
Map[Join[#[[;; 2]], #] &, ballInit];
Table[
t += timeStep;
(* move all balls *)
Do[
dt = (t - balls[[i, 5]]);
balls[[i, 1]] =
balls[[i, 3]] + balls[[i, 4]]*dt - 1/2*g*dt^2;
balls[[i, 2]] =
balls[[i, 4]] - g*dt,
{i, Length@balls}
];
(* handle collisions *)
Do[
(* we'll assume the balls are pre-sorted by height *)
(*
ground collision *)
If[i == 1,
If[balls[[i, 1]] < balls[[i, 6]],
balls[[i, 3]] = balls[[i, 1]] = Max@{balls[[i, 1]], 0};
balls[[i, 4]] = -balls[[i, 2]]*balls[[i, 7]];
balls[[i, 5]] = t;
]
];
(* intra-ball collision *)
If[Length@balls > i,
Do[
If[
balls[[j, 1]] <
balls[[i, 1]] ||
(balls[[j, 1]] - balls[[i, 1]] <
balls[[i, 6]] + balls[[j, 6]]),
balls[[i, 3]] = balls[[i, 1]] =
Max@{Min@{balls[[i, 1]], balls[[j, 1]]}, 0};
balls[[j, 3]] = balls[[j, 1]] =
Max@{Max@{balls[[i, 1]], balls[[j, 1]]}, 0};
balls[[i, 4]] = balls[[j, 2]]*balls[[j, 7]];
balls[[j, 4]] = balls[[i, 2]]*balls[[i, 7]];
balls[[i, 5]] = t;
balls[[j, 5]] = t;
],
{j, i + 1, Length@balls}
]
],
{i, Length@balls}
];
(* return balls *)
balls,
steps
]
]
];

Options[bounceBalls] =
{
"Gravity" ->
QuantityMagnitude[
Quantity[1., "StandardAccelerationOfGravity"],
("Meters"/"Seconds"^2)
],
"TimeStep" -> 1./100,
"Elasticity" -> 1,
"Velocity" -> 0,
"Interpolate" -> False
};
bounceBalls[
t0 : _?NumericQ : 0,
tf_?NumericQ,
ballSpecs : {KeyValuePattern[{"Height" -> _?NumericQ}] ..},
ops : OptionsPattern[]
] :=
Module[{
g, dt, r, e, v,
steps,
balls,
pos
},
{g, dt, r, e, v} =
"Velocity"}];
balls =
Fold[
If[MatchQ[#, KeyValuePattern[{#2[[1]] -> _?NumericQ}]],
#,
Append[#, #2]
] &,
#,
{"InitialTime" -> t0, "Radius" -> r, "Elasticity" -> e,
"Velocity" -> v}
] & /@
ballSpecs;
steps = Floor[(tf - t0)/dt];
MapIndexed[
<|
"Time" -> t0 + dt*#2[[1]],
"Heights" -> #[[All, 1]],
"Velocities" -> #[[All, 2]],
"Elasticities" -> #[[All, 7]]
|> &,
bounceBallsCore[
SortBy[
Lookup[
"Elasticity"}],
#["Height"] &
],
g,
t0,
dt,
steps
]
]
]


And here's a fun result of the generalization. We'll slowly decrease the elasticity of the collision for the bottom balls:

frames =
Table[
With[{
trajData2 =
bounceBalls[
20,
{
<|"Height" -> 157/100, "Elasticity" -> e|>,
<|"Height" -> 2|>
},
"TimeStep" -> 1./500,
"Gravity" -> 10*OptionValue[bounceBalls, "Gravity"]
]
},
ListLinePlot[
Thread[ {#["Time"], #["Heights"]}] & /@ trajData2 // Transpose,
PlotStyle -> AbsoluteThickness[1],
AspectRatio -> 2/10
]
],
{e, 1, .99, -.0005}
];
ListAnimate[frames]


We can see a nice downward trend in the trajectory and some changing of the bounce period (when it doesn't get quashed completely)

• I seriously need to learn how to properly use Mathematica, thanks. The generalization is nice.  I don't need exact results, but I need to be able to measure period of the distance and determine if it's periodic or not, and that would require to know the function with a good precision and certainty – Yuriy S Oct 29 '17 at 20:38
• Makes sense. I think you could also apply some sense of fuzzy-periodicity, too. Like is it periodic within tolerance tol. For that I'd suggest FindTransientRepeat[trajVals, 3, SameTest -> Function[ Abs[Subtract[##]] < tol] ] or something. Post a question about it if you go that route, as someone more expert than I probably has a better way. – b3m2a1 Oct 29 '17 at 20:46
• Parametric plots for positions and/or velocities (i.e. phase space trajectories) are also rather interesting to study, especially in relation to chaos theory, which is why I wanted the explicit functions – Yuriy S Oct 29 '17 at 21:06

General comments on dealing with impacts Your are dealing with a low number of contact points (two). For low number of contact constraints (typically, <10), event driven methods are known to be quite efficient. For larger numbers, time-stepping methods are more recommended. In the present case, my answer relies on an event-driven method, while your is more like a time-stepping method. If you want to dig further into contact dynamics, I recommend using the keywords nonsmooth dynamics.

I know you specifically asked for no alternative method, but then it seems more like a numerics question than a Mathematica question, IMHO. Additionally, I could not refrain myself from posting an answer with WhenEvent because it is so simple :).

The idea of my solution is to use the efficiency of NDSolve for the very simple governing ODEs and use WhenEvent to detect when the bottom mass reaches the ground, or when the two balls meet. Note: I do not know why if I remove the 1; in the definition of events, this no longer works; I asked this question here.

g = 9.81;
eqs = {x1''[t] + g == 0, x2''[t] + g == 0};
ci = {x1[0] == 1, x2[0] == 2, x1'[0] == 0, x2'[0] == 0};
events := {WhenEvent[x1[t] == 0, x1'[t] -> -x1'[t]],
WhenEvent[x2[t] == x1[t], v1 = x1'[t]; x1'[t] -> x2'[t]],
WhenEvent[x2[t] == x1[t], 1; x2'[t] -> v1]}
sol = NDSolveValue[eqs~Join~ci~Join~events, {x1, x2}, {t, 0, 5}]

Plot[{sol[[1]][t], sol[[2]][t]}, {t, 0, 5}]


To view the animated system:

show[{x1_, x2_}] :=
Graphics[{Blue, PointSize[.2], Point[{0, x1}], Red, Point[{0, x2}],
Black, Line[{{-.1, 0}, {.1, 0}}]},
PlotRange -> {{-.1, .1}, {-0.1, 2}}]
tab = Table[show[Through[sol[t]]], {t, Subdivide[0, 2, 99]}];
Animate[show[Through[sol[t]]], {t, 0, 2}]


Periodic solutions I don't think the solution you provide corresponds to a periodic solution. Indeed, assuming all impacts perfectly alternate between the ground, the balls, the ground, etc. then the problem reduces to finding two "bouncing" parabolas which, after a certain amount of time, recover the initial state. In short, look at the above plot and think of it in black and white: it is simply two parts of parabolas replicated.

Assuming for simplicity that the initial velocities are zero, the first time when each parabola meets the ground is $t=\sqrt{2x_i(0)}{g}$. The solution is periodic if $\sqrt{x_2(0)/x_1(0)}\in\mathbb{Q}$. I am not sure if this sufficient condition is also necessary, because of the assumptions.

For example, for $x_1(0)=1$ and $x_2(0)=4$ ($\sqrt{4}\in\mathbb{Q}$), we get:

For $x_1(0)=16$ and $x_2(0)=25$ ($\sqrt{25/16}\in\mathbb{Q}$), we get another periodic solution:

However, for your initial conditions, $\sqrt{2/1.5} = 2/\sqrt{3}\not\in\mathbb{Q}$.

• It's really is much more simple. thank you. I just didn't know how to use WhenEvent correctly. – Yuriy S Oct 29 '17 at 21:18
• @YuriyS My pleasure. Out of curiosity, may I ask why you are looking for periodic solutions? – anderstood Oct 29 '17 at 21:38
• @ anderstood, simple curiosity, I haven't found any literature on this exact problem, and it's quite different compared to well known double pendulum – Yuriy S Oct 29 '17 at 21:39
• @YuriyS That's an interesting problem. I happen to have worked on a pretty similar problem, but only one contact, and no external forces. The work is available here and you can get an idea of it by checking the animations in the lower part of the webpage. – anderstood Oct 29 '17 at 21:42
• @YuriyS I edited the answer to add thoughts about periodic solutions. – anderstood Oct 31 '17 at 16:17