4
\$\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.

The code is written in Mathematica and to be honest, it was my first time writing a piece of code in this language (which wasn't a simple calculation of an integral or something alike). It works and produces results I would expect but it incredibly slow and I don't even know where to start to correct it.

So what I'm looking for are some general hints on how to imporve when writing Mathematica code, style conventions (is there something like pep8 for Python?) and how to speed the code up.


ClearAll["Global`*"]

MapMonitored[f_, args_List] := Module[{x = 0}, 
  Monitor[MapIndexed[(x = #2[[1]]; f[#1]) &, args], 
   ProgressIndicator[x/Length[args]]]]

DrawBread[sol_,t_]:=Module[{bb=0.02,LL=0.1},ut=Flatten[u[t]/.sol];

  xb:=trajectory[t]/.{b\[Rule]0.02,L\[Rule]0.1}+{-bb/2.0 \
  Sin[ut],-bb/2.0 Cos[ut]};

  xleftdown:=xb+{-LL/2.0 Cos[ut],LL/2.0 Sin[ut]};

  xleftup:=xleftdown+{bb Sin[ut],bb Cos[ut]};

  xrightup:=xleftup+{LL Cos[ut],-LL Sin[ut]};

  xrightdown:=xrightup+{-bb Sin[ut],-bb Cos[ut]};

  p=Block[{b=0.02,L=0.1},
    Show[
      ListLinePlot[
        {Flatten[{xb,xleftdown,xleftup,xrightup,xrightdown,xb}]},
        PlotStyle\[Rule]Blue
      ],
      ListLinePlot[
        {Flatten[xleftup],Flatten[xrightup]},
        PlotStyle[Rule]Red
      ],
      ListLinePlot[
        {Flatten[xrightup],Flatten[xrightdown],Flatten[xb]},
        PlotStyle\[Rule] Blue
      ],
      ListPlot[
        {0.0,trajectory[t]},
        PlotStyle\[Rule]Red
      ],
      AspectRatio\[Rule]1.0,
      PlotRange\[Rule]{{-0.1,1.0},{0.05,-1.5}}
    ]
  ];

  Return[p]];

DrawBreadAfterChange[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]};

  p=Block[{b=0.02,L=0.1},
   Show[
     ListLinePlot[
       {Flatten[xb],Flatten[xleftdown],Flatten[xleftup],Flatten[xrightup],
       Flatten[xrightdown],Flatten[xb]},
       PlotRange\[Rule]{{-0.1,1.0},{0.05,-1.5}},
       PlotStyle\[Rule]Blue
     ],
     ListLinePlot[
       {Flatten[xleftup],Flatten[xrightup]},
       PlotRange\[Rule]{{-0.1,1.0},{0.05,-1.5}},
       PlotStyle\[Rule]Red
     ],
     ListLinePlot[
       {Flatten[xrightup],Flatten[xrightdown],Flatten[xb]},
       PlotRange\[Rule]{{-0.1,1.0},{0.05,-1.5}},
       PlotStyle\[Rule]Blue
     ],
     ListPlot[{Flatten[wurfparabel[t,sol,tst,vtstar]]},
       PlotStyle\[Rule]Red],
       AspectRatio\[Rule]1.0
     ]
   ];
  Return[p]
];

DrawBreadFinal[sol_,t_,tst_,vtstar_]:=If[t\[LessEqual]tst,
  DrawBread[sol,t],
  DrawBreadAfterChange[sol,t,tst,vtstar]
];

CalculateEdgesOfBread[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[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}]];


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`

GetTheThingsINeed[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}];

  trajectory[t_] := Flatten[x[t] /. solution];

  v[t_] = Flatten[D[trajectory[t], t]];

  tGroundTime = 
   t /. Flatten[
     FindRoot[{CalculateEdgesOfBread[solution, t, tstar, v[tstar]] - 
        TabelHeight}, {t, 0.5}]];

  Return[BreadOrButterSide[
    AngleAtGround[solution, tstar, tGroundTime]]]];

velocities = Range[0.1, 1.8, 0.2];
TabelHeights = Range[0.3, 1.5, 0.8];
mesh = Tuples[{velocities, TabelHeights}];

plotmesh = 
 MapMonitored[GetTheThingsINeed, mesh] // 
  AbsoluteTiming; plotmesh[[1]]

ListPlot[Lookup[
  GroupBy[
   Thread[{mesh, plotmesh[[2]]}],
   Last -> First],
  {1.0, 0.0}],
 PlotStyle -> {Green, Red}, AxesLabel -> {"v", "h"}
 ]

Edit: I cleared the code up a bit.. Didn't impact performance at all.

Edit II: On a whim I just changed the line

tGroundTime = 
   t /. Flatten[
     NSolve[{CalculateEdgesOfBread[solution, t, tstar, v[tstar]] == 
        TabelHeight, t > 0.0}, t]];

to

tGroundTime = 
  t /. Flatten[
    FindRoot[{CalculateEdgesOfBread[solution, t, tstar, v[tstar]] - 
       TabelHeight}, {t, 0.5}]];

and the Code is now approximately 1000 times faster.. Not sure why this makes such a huge change..

EDIT III: Cleared the code up again and moved as much as possible out of the GetTheThingsINeed out. Turns out this improves the performance again, now the code runs in about 0.8 sec compared to the 12 sec in the beginning.

\$\endgroup\$
  • 1
    \$\begingroup\$ NSolve was searching for all possible positive solutions, whereas FindRoot uses something like Newton's method to find a single solution. \$\endgroup\$ – KraZug Apr 27 at 19:15

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.