0
$\begingroup$

Here is a whole block of code that computes a function at the point vx0=2.2 and returns 61.9113, which is much lower than the point found by FindMinimum at vx0=3, where the function is 98771.5. The Sow and Reap shows that FindMinmum is not doing a proper scan of the function, other methods and NMinimize always return vx0=3.

Can FindMinimum and NMinimize be so wrong? how to get it to find the actual minimum (that is close to 2.43)?

g[a_, initial_] := (ThrowBall[a, initial][[1]] - 25)^2

Clear[ThrowBall]
Options[ThrowBall] = {spacing -> 0.01, debug -> False}

ThrowBall[\[Alpha]_, initial_, opt : OptionsPattern[]] := 
 Module[{debug = OptionValue[debug], sol, traj, 
   trajs, $xdot, $z, $zdot, t1, timpact, ts, $timpact}, 
  sol = DSolve[{x''[t] == 0, z''[t] == -g, x[0] == initial["x"], 
     z[0] == initial["z"], x'[0] == initial["xdot"], 
     z'[0] == initial["zdot"]}, {x[t], z[t]}, t];
  If[debug, Print[x[t] /. sol[[1]]];
   Print[D[x[t] /. sol, t]];
   Print[D[z[t] /. sol, t]];];
  x$expr = ReplaceAll[x[t], sol[[1]]];
  $x[t_] = x$expr;
  z$expr = ReplaceAll[z[t], sol[[1]]];(*/.g->gN*)$z[t_] = z$expr;
  ts = {0};
  trajs = {};
  If[debug, Print[Simplify[$x[ts[[-1]]]/m]];];
  (**)For[i = 1, Simplify[$x[ts[[-1]]]/m] < 25, 
   i++, $timpact = Solve[ReplaceAll[z[t] == 0, sol], t];
   timpact = 
Assuming[{s > 0, m > 0}, Simplify[t /. $timpact //. g -> gN]];
   (*t1=Select[timpact,(D[#,s]>D[ts[[-1]],
   s])&][[1]];*)(*prone to numerical confusion between the last time \
of rebounce and the smallest zero fo the new parabola*)
   t1 = MaximalBy[Select[timpact, (D[#, s] > D[ts[[-1]], s]) &], 
      D[#, s] &][[1]];
   AppendTo[ts, t1];
   If[debug, Print[ts];];
   x$expr = x[t] /. sol[[1]] /. g -> gN;
   $x[t_] = x$expr;
   z$expr = z[t] /. sol[[1]] /. g -> gN;
   $z[t_] = z$expr;
   x$dot = D[x[t] /. sol, t] /. g -> gN;
   $xdot[t_] = x$dot;
   z$dot = D[z[t] /. sol, t] /. g -> gN;
   $zdot[t_] = z$dot;
   If[debug, Print[{i, \[Alpha]}]];
   sol = 
    DSolve[{x''[t] == 0, z''[t] == -g, x[t1] == $x[t1], z[t1] == 0 m, 
  x'[t1] == $xdot[t1], z'[t1] == -\[Alpha]*$zdot[t1]}, {x[t], 
  z[t]}, t];
   traj = 
Table[{Simplify[$x[t*s]/m], Simplify[$z[t*s]/m]}, {t, ts[[-2]]/s, 
  ts[[-1]]/s, OptionValue[spacing]}];
   AppendTo[trajs, traj];];
  $t25 = t /. Solve[$x[t] == 25 m, t][[1, 1]];
  z25 = Simplify[$z[t]/m /. g -> gN /. t -> $t25];
  If[debug, Print[z25]];
  {z25, Flatten[trajs, 1]}]

td = TruncatedDistribution[{0.84, 1}, NormalDistribution[0.9, 0.02]]

gN = 9.8 m/s^2

Inaccuracy[x0_, vx0_, z0_] := 
 0.16*Mean[
   Table[
    g[a, <|"x" -> x0 m, "z" -> z0 m, "xdot" -> vx0 m/s, 
       "zdot" -> 0|>]*PDF[td, a], {a, 0.84, 1, 0.001}]]

Inaccuracy[0, 2.2, 50]
pts = Reap[
  NMinimize[{Inaccuracy[0, vx0, 50], vx0 < 3}, {vx0}, 
   StepMonitor :> Sow[{vx0}]]]
```
$\endgroup$
1
  • $\begingroup$ Comments are not for extended discussion; this conversation has been moved to chat. $\endgroup$
    – Kuba
    Jan 10, 2022 at 7:25

1 Answer 1

6
$\begingroup$

I think I understand the point of most of the code. It seems to me that you're throwing a ball—with some initial $(x_0,y_0)$ position, velocity $(v_{x,0},v_{y,0})$, and acceleration $(0,a_y=-g)$—that bounces off the ground with a coefficient of restitution $\alpha$. When the ball has travelled $x_\textrm{target}=25\textrm{ m}$ horizontally, you want to know the height of the ball above the ground. Then, you have some objective function (that I don't understand, but I don't think I need to) and wish to minimize the objective function.

The Problem

I believe the issue is that Mathematica often tries to do some fancy symbolic analysis before beginning the minimization routine. Because of this, it first feeds the vx0 into your Inaccuracy function symbolically. If you just run Inaccuracy[0, vx0, 50] Mathematica will return 0.9936758426873792*(25. - 3062.5/vx0^2)^2. It then feeds this to the minimizer, which quickly returns an inaccuracy of 98771.5 at vx0 = 3. As far as it knows, it has successfully minimized your equation.

To fix this, we just have to define

Inaccuracy[x0_, vx0_?NumericQ, z0_] := 
 0.16*Mean[
   Table[
    g[a, <|"x" -> x0 m, "z" -> z0 m, "xdot" -> vx0 m/s, 
       "zdot" -> 0|>]*PDF[td, a], {a, 0.84, 1, 0.001}]]

as shown on What are the most common pitfalls awaiting new users. ?NumericQ prevents the function from evaluating unless the argument is explicitly a number. With this fixed, it takes 1267 seconds to return the correct answer on my machine. I had almost given up.

I think there's a number of ways to approach this, but my main consideration is that the way your Inaccuracy function is calculated, the throw has to be simulated 161 times for each guess of $v_{x,0}$ given the $\alpha$ values used in this example. This is probably going to be the most time consuming portion, so we want it as fast as possible. Right now, it takes about 17 seconds for one call of Inaccuracy to return on my computer.

Using NDSolve

If we're willing to use high-level functions like DSolve and NDSolve, then I would recommend programming the bounce into them directly. For brevity, I'll just show NDSolve. DSolve is similar, but has a few caveats like needing to prevent execution while $\alpha$ is still just a symbol. Since you use $a_x = 0\textrm{ m/s}^2$ in yours, I have done so as well though it wouldn't be too difficult to generalized. Technically, I don't really need x[t] due to the way I calculate the time at which the ball hits the target, but thought I would leave it in case you decided you need non-zero $x$-acceleration.

(* returns the height above the ground when the ball has
travelled a horizontal distance of xtarget *)
throwND[xtarget_, \[Alpha]_, x0_, y0_, vx0_, vy0_, ay_] := 
  (NDSolveValue[{
     x''[t] == 0,
     y''[t] == ay,
     x'[0] == vx0,
     y'[0] == vy0,
     x[0] == x0,
     y[0] == y0,
     WhenEvent[y[t] == 0, y'[t] -> -\[Alpha] y'[t]]
     },
    y[t],
    {t, 0, (xtarget - x0)/vx0}
    ]) /. t -> (xtarget - x0)/vx0

(* Set the above function to be listable so I can send all 161
alpha values simultaneously *)
SetAttributes[throwND, Listable]

(* Precalculate alpha values and the distribution; no need to
redo it on every call of the minimizer *)
\[Alpha]vals = Range[0.84, 1, 0.001];
td = PDF[
  TruncatedDistribution[
    {0.84, 1}, 
    NormalDistribution[0.9, 0.02]
  ], 
\[Alpha]vals];

(* The objective function. I've preprogrammed some of the numbers
to match your own version, but you can adjust this *)
inaccND[x0_, vx0_?NumericQ, y0_] := 
  0.16*
  Mean[(throwND[25, \[Alpha]vals, x0, y0, vx0, 0, -9.8] - 25)^2*td]

(* Call FindMinimum *)
FindMinimum[inaccND[0, vx0, 50], vx0, Method -> "ConjugateGradient"]
FindMinimum[{inaccND[0, vx0, 50], 2 < vx0 < 3}, {vx0, 2.2}]

(* Both return {0.0830755, {vx0 -> 2.43019}} *)

The conjugate gradient method takes about 39 seconds, and the constrained method takes about 27 seconds.

I don't know why it has issues with your code, as both inaccuracy functions return the same results given the same input. The most confounding thing is that it returns 98771.5 as the minimum when $v_{x,0}=3\textrm{ m/s}$ when Inaccuracy[0, 3., 50] returns 101.7 or so. However, I would recommend rethinking your approach to the problem both for speed and ease of reading.

Using Algebra

I can get a further 100x speed-up by working through the math and finding an algebraic representation, and my FindMinimum runs in about 0.27 seconds. I could probably shave further time off of this by compiling the functions (EDIT: I just checked and I can get it down to 0.027 seconds for the minimization after compiling), but my main interest was the math behind it so I didn't take it that far. This post is already pretty long, but I'll post the code and try to give a brief rundown. Basically, if you work it out, the formula for the $n$th parabola is:

$y_n(t) = \frac{1}{2} a_y \left(t + \frac{v_{y,0}}{a_y} - \sqrt{\frac{-2y_{p,0}}{a_y}}\sum_{i=1}^n (\alpha^{i-1}+\alpha^i)\right)^2 + \alpha^{2n}y_{p,0}$,

where $y_{p,0} = y_0 - \frac{v_{y,0}^2}{2 a_y}$ is the height of the $n$th peak, and the sum has a nice closed form. Under my notation, the parabola that contains the point of release is $n=0$. We need a helper function to determine what our $n$ is at a given time. Basically, I set everything inside the brackets of $\frac{1}{2} a_y (\textrm{this stuff here})^2$ equal to zero and solve for $n$. This gives us a function with a bunch of logs. One important note is that these fail for $\alpha = 1$, but I solve for the limit and they have a solution. So if $\alpha = 1$, I just use that solution. The code, with minimal commentary, is:

(* Function that determines n *)
myn[t_, \[Alpha]_, y0_, vy0_, ay_] := If[\[Alpha] == 1, 
     (ay*(Sqrt[(vy0^2 - 2*ay*y0)/ay^2] + t) + vy0)/(2*ay*
     Sqrt[(vy0^2 - 2*ay*y0)/ay^2]), 
     (Log[\[Alpha] + ((\[Alpha] - 1)*(ay*t + vy0))/(ay*
          Sqrt[(vy0^2 - 2*ay*y0)/ay^2]) + 1] - Log[2])/
       Log[\[Alpha]]]

(* Function to calculate the nth parabola *)
myy[n_, \[Alpha]_, y0_, vy0_, ay_, t_] := 
 Module[{yp0 = y0 - vy0^2/(2*ay)}, 
     If[\[Alpha] == 
    1, ((ay*(t - 2*n*Sqrt[(vy0^2 - 2*ay*y0)/ay^2]) + vy0)^2 + 
      2*ay*y0 - vy0^2)/
         (2*ay), (1/2)*
     ay*(t + vy0/ay - 
        Sqrt[(-2*yp0)/
           ay]*(((\[Alpha] + 1)*(\[Alpha]^n - 1))/(\[Alpha] - 1)))^
             2 + \[Alpha]^(2*n)*yp0]]

SetAttributes[myn, Listable]
SetAttributes[myy, Listable]

(* The rest of this code is the same as before, I just include
it here for easy copy-and-paste *)
throw[xtarget_, \[Alpha]_, x0_, y0_, vx0_, vy0_, ay_] := Module[
  {t = (xtarget - x0)/vx0, n},
  n = myn[t, \[Alpha], y0, vy0, ay];
  myy[Floor[n], \[Alpha], y0, vy0, ay, t]
  ]
\[Alpha]vals = Range[0.84, 1, 0.001];
td = PDF[
   TruncatedDistribution[{0.84, 1}, 
    NormalDistribution[0.9, 0.02]], \[Alpha]vals];
inacc[x0_, vx0_?NumericQ, y0_] := 
 0.16*Mean[(throw[25, \[Alpha]vals, x0, y0, vx0, 0, -9.8] - 25)^2*td]
FindMinimum[inacc[0, vx0, 50], vx0, Method -> "ConjugateGradient"]

(* {0.0830755, {vx0 -> 2.43019}} *)

Conclusion

One nice thing about either of these methods is that the inaccuracy function is fast enough that you can readily plot the function to confirm visually that it seems to be near the minimum.

Plot[
  inacc[0, vx0, 50], 
  {vx0, 2, 3}, 
  Axes -> False, 
  Frame -> True, 
  FrameLabel -> {"\!\(\*SubscriptBox[\(v\), \(x, 0\)]\)(m/s)", "Inaccuracy"}, 
  GridLines -> {{2.43019}, None}, 
  PlotRange -> {{2, 3}, {-5, 510}}
]

Plot of inaccuracy function vs vx0.

$\endgroup$
4
  • $\begingroup$ Thanks! I see the problem InaccuracySymb[0, vx, 25] evalues to something incorrect. I will need to go deeper as of what is going on internally, my function was clearly set up to solve everything numerically and detect the bounce and then the next bounce etc. Testing this semi-symbolic input was the key to figure out where the nonsense was stemming from! I have learned something. $\endgroup$
    – Rho Phi
    Jan 15, 2022 at 20:53
  • $\begingroup$ Regarding the version of the code with ?NumericQ I am a bit puzzled by the fact that a single evaluation takes 12 seconds, then I ask for the FindMinimum and it takes 900 seconds. I am puzzled because the points that are collected by Reap after Saw are just 8 {{{2.33}, {2.46036}, {2.43279}, {2.43011}, {2.43019}, {2.43019}, {2.43019}, {2.43019}}}}}. I had run with MaxIterations -> 10, not clear to me if this 10 is related to the number of points I reap from sow, or it's yet another thing. $\endgroup$
    – Rho Phi
    Jan 15, 2022 at 20:59
  • $\begingroup$ Yes, this case is simple enough analytical values can be found, the real world case will not enjoy that. In any case I will study NDSolveValue as it seems pretty useful! $\endgroup$
    – Rho Phi
    Jan 15, 2022 at 21:06
  • $\begingroup$ Just a comment, Inaccuracy was given as an Mean of evalutation of the objective function as to mimic a rough NIntegrate of g. That was my workaround for NIntegrate shooting other nonsense ... there I was using symbolds and numbers by design, e.g. x0 contained a symbol m to represent units of measurement. I will check what happens without that, probably the mix of symbols and numbers caused some trouble to NIntegrate there. In case I will open a new question. $\endgroup$
    – Rho Phi
    Jan 15, 2022 at 21:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.