6
$\begingroup$

I've written a piece of code that determines if a piece toast with butter lands on the butter-side or not, depending on its initial velocity and the table height. If you are interested in the physics behind this, see here.

The code works and returns the results that I want/would expect, but its pretty slow. I tried improving the code using some tips from here but I'm not really sure if I did it correctly. I already asked a question on CodeReview SE, but I didn't get any responses see here.

Since then I made some improvements on the code and it runs already a lot faster than in the beginning, but I'm still wondering if and how it can be improved.

Please not that the code needs about 60 seconds on my machine with the current settings.

Any other kind of suggestions on how to improve the code are of course welcome as well.


ClearAll["Global`*"]
(* Added a Progress Bar to the Map function *)
MapMonitored[f_, args_List] := 
 Module[{x = 0}, 
  Monitor[MapIndexed[(x = #2[[1]]; f[#1]) &, args], 
   ProgressIndicator[x/Length[args]]]]

LowestPositionOfBread[sol_, t_, tst_, vtstar_] := 
  Module[{bb = 0.02, LL = 0.1},
   uu[time_] = Flatten[u[time] /. sol];
   uu0 = Abs[uu[tst] - uu'[tst] tst + uu[tst]];
   angle = uu'[tst] t - uu[tst] + uu0;
   xb = Wurfparabel[t, sol, tst, 
      vtstar] + {-bb/2.0 Sin[angle], -bb/2.0 Cos[angle]};
   xleftdown = xb + {-LL/2.0 Cos[angle], LL/2.0 Sin[angle]};
   xleftup = xleftdown + {bb Sin[angle], bb Cos[angle]};
   xrightup = xleftup + {LL Cos[angle], -LL Sin[angle]};
   xrightdown = xrightup + {-bb Sin[angle], -bb Cos[angle]};
   Return[Abs[Min[xleftdown, xleftup, xrightup, xrightdown]]]
   ];

AngleAtGround[sol_, tst_, tt_] := Module[{},
   uu[t_] = Flatten[u[t] /. sol];
   uu0 = Abs[uu[tst] - uu'[tst] tst + uu[tst]];
   angle = uu'[tst] tt - uu[tst] + uu0;

   Return[Mod[angle[[1]], 2.0 Pi]]
   ];

BreadOrButterSide[phi_] := 
  If[ (* Returns 1.0 for ButterSide and 0.0 if not *)
   Pi/2.0 <= phi,
   If[phi <= 3.0 Pi/2.0,
    1.0,
    0.0],
   0.0
   ];

Wurfparabel[tt_, sol_, tst_, vtstar_] := 
  Module[{g = 9.81}, beta = Abs[Flatten[u[tst] /. sol]];
   xspace = vtstar[[1]] tt Cos[beta];
   xx0 = Abs[trajectory[tst][[1]] - vtstar[[1]] tst Cos[beta]];
   yspace = -(Norm[v[tst]] tt Sin[beta] + g tt^2/2.0);
   y0 = Abs[
     trajectory[tst][[2]] + (Norm[v[tst]] tst Sin[beta] + 
        g tst^2/2.0)];
   Return[{xspace + xx0, yspace + y0}]
   ];

(* Implementation of the physics behind the Problem 
if interested in details see: \
https://physics.stackexchange.com/questions/474412/block-sliding-down-\
a-table *)
x[t_] = {{s[t] Cos[u[t]] + b/2 Sin[u[t]]}, {-s[t] Sin[u[t]] + 
     b/2 Cos[u[t]]}};
m = 0.1; l = 0.1; g = 9.81;
J = 1.0/12.0 m (l^2 + b^2);

T[t_] = Simplify[1.0/2.0 (m) Flatten[x'[t]].Flatten[x'[t]]] + 
   1.0/2.0 J u'[t]^2;
V[t_] = m g x[t][[2, 1]];
L[t_] = Simplify[T[t] - V[t]];

<< VariationalMethods`
FunctionWrapper[StartingValuesList_, Lagrangian_: L, 
  PositionVector_: x] := Module[{g = 9.81},
  b = 0.02 ;
  eoms[t_] = Simplify[VariationalD[Lagrangian[t], s[t], t]];
  eomu[t_] = Simplify[VariationalD[Lagrangian[t], u[t], t]];
  TabelHeight = StartingValuesList[[2]];
  solution = 
   NDSolve[{eoms[t] == 0.0, eomu[t] == 0.0, u[0] == 0.0, s[0] == 0.01,
      s'[0] == StartingValuesList[[1]], u'[0] == 0.1}, {u, s}, {t, 
     0.0, 10.0}];
  f[t_] = D[x[t][[1, 1]], {t, 2}] /. solution;
  tstar = t /. FindRoot[f[t], {t, 0.1}]; (* 
  Time at which Bread loses contact to the table *)
  trajectory[t_] := Flatten[x[t] /. solution]; (* 
  Trajectory of the Center of mass before the contact is lost *)
  v[t_] = Flatten[D[trajectory[t], t]]; (* 
  Velocity vector of center of mass *)
  tGroundTime = 
   t /. Flatten[
     FindRoot[{LowestPositionOfBread[solution, t, tstar, v[tstar]] - 
        TabelHeight}, {t, 0.5}]];
  Return[BreadOrButterSide[
    AngleAtGround[solution, tstar, tGroundTime]]]]

velocities = Array[# &, 75, {0.01, 1.8}];
TabelHeights = Array[# &, 40, {0.2, 2.5}];
mesh = Tuples[{velocities, TabelHeights}];

plotmesh = 
 MapMonitored[FunctionWrapper, mesh] // 
  AbsoluteTiming; Print["Computation time ", plotmesh[[1]], " sec"]

PointsOfInterest = GroupBy[
    Thread[{mesh, plotmesh[[2]]}],
    Last -> First][[2]];
SmallestValues = 
  Join @@ Values@GroupBy[PointsOfInterest, First, MinimalBy[Last]];
BiggestValues = 
  Join @@ Values@GroupBy[PointsOfInterest, First, MaximalBy[Last]];
LowerBound = Fit[SmallestValues, {1, x, x^2, x^3, x^4}, {x}];
UpperBound = 
  Fit[Select[BiggestValues, #[[2]] != TabelHeights[[-1]] &], {1, x, 
    x^2, x^3, x^4}, {x}];

<< MaTeX`

Show[Plot[{LowerBound, UpperBound}, {x, velocities[[1]], 
   velocities[[-1]]}, Filling -> {1 -> {{2}, {Green, None}}}, 
  PlotRange -> {{0.01, 1.8}, {0.0, TabelHeights[[-1]]}}, 
  FillingStyle -> Opacity[0.2]], 
 Plot[LowerBound, {x, velocities[[1]], velocities[[-1]]}, 
  Filling -> Axis, PlotStyle -> Red], 
 Plot[UpperBound, {x, velocities[[1]], velocities[[-1]]}, 
  Filling -> Top, PlotStyle -> Red], 
 Frame -> True, FrameLabel -> MaTeX /@ {"v", "h"}]
$\endgroup$
  • 4
    $\begingroup$ Note that there's GeneralUtilities`MonitoredMap which also provides a time estimate. $\endgroup$ – Chip Hurst May 4 at 21:01
  • 15
    $\begingroup$ One way to get a faster simulation would be to butter both sides... $\endgroup$ – Daniel Lichtblau May 4 at 21:14
  • 1
    $\begingroup$ In the CodeReview post, you reference getting it down to 0.6 seconds, is this true here, now? Very excited to tear this code apart and see how it ticks :) $\endgroup$ – CA Trevillian May 9 at 14:20
  • 1
    $\begingroup$ @CATrevillian I‘m not at home right now so I‘m not sure about the processor but 16GB of ram :) $\endgroup$ – Sito May 9 at 14:28
  • 1
    $\begingroup$ @CATrevillian just checked it, 4 cores and 16 gb. $\endgroup$ – Sito May 9 at 18:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.