Problem with For Loop when evaluating a recursive function at multiple points

Question

I am running into problems with my For loop when trying to evaluate my expression. Fix some parameter values :

{a = 0.6958877284595301, pl = 0.153415660318751`, 
ph = 0.5548313351573813`, U0 = 1, d = 0.9, 
Uh = (a (1 + d (-1 + ph)) U0)/(1 + d (-1 + a ph)), 
Ul = (a (1 + d (-1 + pl)) U0)/(1 + d (-1 + a pl)), R = 3, c = 1, 
Pl = R - c Ul, Ph = R - c Uh, u = Log[ph/pl], 
dow = -Log[(1 - ph)/(1 - pl)], do = -Log[(1 - a ph)/(1 - pl)]};

(those parameters are chosen such that u, dow, and do are multiples of one another .)

and further define :

f1 = Table[{l, Max[(ph Exp[l] + pl) Ph, (pl) Pl]}, {l, -50, +50, 
0.321378298772996`}];

I have a recursive expression which I am evaluating at different values of the variable l . When I set the bounds in the For - Loop such that only one value of l satisfies the bounds given, the For Loop gives me the correct answer .

For[l = -50. + dow , l <= -50. + dow, l = l + 0.321378298772996`, 
Print[Max[(pl + Exp[l] ph) Ph + 
d pl f1[[Part[Flatten[Position[f1, l + u]], 1], 2]] + 
d (1 - pl) f1[[Part[Flatten[Position[f1, l - dow]], 1], 2]], 
pl (R - 
   c (((1 - d + a d pl) Ul)/(1 - d) - 
      d pl/(1 - d)
        Ul + (1 - 
         a) d (Piecewise[{{pl/(1 - d) Ul , 
           l + u < Log[pl/ph (Pl - Ph)/Ph]}, { 
           pl Uh + d pl/(1 - d) Ul, 
           l + u > Log[pl/ph (Pl - Ph)/Ph]}}]))) + 
d pl pl Pl (1 - d^1)/(1 - d) + 
d (1 - pl) f1[[Part[Flatten[Position[f1, l + do]], 1], 2]]]]]

0.628187

On the other hand, if widen the bounds for which l (which would define the second round of the recursion, call it f2), I am getting the following error messages :

For[l = -50. + dow , l <= 50. - u, l = l + 0.321378298772996`, 
Print[Max[(pl + Exp[l] ph) Ph + 
d pl f1[[Part[Flatten[Position[f1, l + u]], 1], 2]] + 
d (1 - pl) f1[[Part[Flatten[Position[f1, l - dow]], 1], 2]], 
pl (R - 
   c (((1 - d + a d pl) Ul)/(1 - d) - 
      d pl/(1 - d)
        Ul + (1 - 
         a) d (Piecewise[{{pl/(1 - d) Ul , 
           l + u < Log[pl/ph (Pl - Ph)/Ph]}, { 
           pl Uh + d pl/(1 - d) Ul, 
           l + u > Log[pl/ph (Pl - Ph)/Ph]}}]))) + 
  d pl pl Pl (1 - d^1)/(1 - d) + 
  d (1 - pl) f1[[Part[Flatten[Position[f1, l + do]], 1], 2]]]]]

0.628187

0.628187

Part::partw: Part 1 of {} does not exist.

Part::pkspec1: The expression {}[[1]] cannot be used as a part specification.

Part::pkspec1: The expression {}[[1]] cannot be used as a part specification.

Part::partw: Part 1 of {} does not exist.

Part::pkspec1: The expression {}[[1]] cannot be used as a part specification.

General::stop: Further output of Part::pkspec1 will be suppressed during this calculation.

Part::partw: Part 1 of {} does not exist.

General::stop: Further output of Part::partw will be suppressed during this calculation.

I tried to split it up in the different expressions to see where the problem occurs - see below . I am not able to make sense of this . Can anyone help? Thanks so much!

For[l = -50. + dow , l <= -50. + 3 dow, l = l + 0.321378298772996`, 
Print[f1[[Part[Flatten[Position[f1, l + do]], 1], 2]]]]

0.330625

0.330625

0.330625

0.330625

0.330625

 For[l = -50. + dow , l <= -50. + 4 dow, l = l + 0.321378298772996`, 
 Print[f1[[Part[Flatten[Position[f1, l + do]], 1], 2]]]]

0.330625

0.330625

0.330625

0.330625

0.330625

Part::partw: Part 1 of {} does not exist.

Part::pkspec1: The expression {}[[1]] cannot be used as a part specification.

Part::partw: Part 1 of {} does not exist.

Part::pkspec1: The expression {}[[1]] cannot be used as a part specification.

For[l = -50. + dow , l <= -50. + 3*do, l = l + 0.321378298772996`, 
Print[f1[[Part[Flatten[Position[f1, l - dow]], 1], 2]]]]

0.330625

0.330625

For[l = -50. + dow , l <= -50. + 6*do, l = l + 0.321378298772996`, 
Print[f1[[Part[Flatten[Position[f1, l - dow]], 1], 2]]]]

0.330625

0.330625

Part::partw: Part 1 of {} does not exist.

Part::pkspec1: The expression {}[[1]] cannot be used as a part specification.

Part::partw: Part 1 of {} does not exist.

Part::pkspec1: The expression {}[[1]] cannot be used as a part specification.

Part::partw: Part 1 of {} does not exist.

General::stop: Further output of Part::partw will be suppressed during this calculation.

Part::pkspec1: The expression {}[[1]] cannot be used as a part specification.

General::stop: Further output of Part::pkspec1 will be suppressed during this calculation.

For[l = -50. + dow , l <= -50. + 3*do, l = l + 0.321378298772996`, 
Print[f1[[Part[Flatten[Position[f1, l + u]], 1], 2]]]]

0.330625

0.330625

For[l = -50. + dow , l <= -50. + 10*do, l = l + 0.321378298772996`, 
Print[f1[[Part[Flatten[Position[f1, l + u]], 1], 2]]]]

0.330625

0.330625

Part::partw: Part 1 of {} does not exist.

Part::pkspec1: The expression {}[[1]] cannot be used as a part specification.

Part::partw: Part 1 of {} does not exist.

Part::pkspec1: The expression {}[[1]] cannot be used as a part specification.

Part::partw: Part 1 of {} does not exist.

General::stop: Further output of Part::partw will be suppressed during this calculation.

Part::pkspec1: The expression {}[[1]] cannot be used as a part specification.

General::stop: Further output of Part::pkspec1 will be suppressed during this calculation.

Answer 1

I'm sorry, I have to confess I have absolutely no idea what's going on in your program, so I'm not sure if I can help or not. That said, my guess is that you're running into problems because of rounding/accuracy inconsistency; when your code can't find the number it expects in f1, that's when you get that error.

I tried changing your second 'widened' version just a little:

For[l = -50. + dow, l <= 50. - u, l = l + 0.321378298772996`, 
 Print[Short[
   Max[(pl + Exp[l] ph) Ph + 
     d pl f1[[Part[
         Flatten[Position[Round[f1, 0.0001], Round[l + u, 0.0001]]], 
         1], 2]] + 
     d (1 - pl) f1[[Part[
         Flatten[Position[Round[f1, 0.0001], Round[l - dow, 0.0001]]],
          1], 2]], 
    pl (R - c (((1 - d + a d pl) Ul)/(1 - d) - 
           d pl/(1 - d) Ul + (1 - 
              a) d (Piecewise[{{pl/(1 - d) Ul, 
                l + u < Log[pl/ph (Pl - Ph)/Ph]}, {pl Uh + 
                 d pl/(1 - d) Ul, 
                l + u > Log[pl/ph (Pl - Ph)/Ph]}}]))) + 
     d pl pl Pl (1 - d^1)/(1 - d) + 
     d (1 - pl) f1[[Part[
         Flatten[Position[Round[f1, 0.0001], Round[l + do, 0.0001]]], 
         1], 2]]]]]]

Just rounding to the nearest 0.0001 every time it's checking that array. That is not a great solution, just a workaround/patch/test, but that makes it so there are no longer any errors, and the results seem to start increasing, up to around $3 \times 10^{21}$. I have no idea if this is the behavior you're looking for.

If you'd like more helpful help, it makes life much easier on commenters/reviewers if there are variable names, or comments, or something else to figure out what's going on. Just a friendly suggestion for next time. Anyway, hopefully this helped!