3
$\begingroup$

I have the following recurrence equation

p[t_, s_] := 
 p[t, s] = ((1 + s)*p[t - 1, s]^2 + (1 + s)*
    p[t - 1, s]*(1 - p[t - 1, s]))/((1 + s)*p[t - 1, s]^2 + (1 + s)*2*
    p[t - 1, s]*(1 - p[t - 1, s]) + (1 - p[t - 1, s])^2)
p[0, s_] := 0.00001

and I would like to solve for when

p[50, s] == 0.8

One approach that worked for me was to first construct a list

listvals = Table[{s, p[50, s]}, {s, 0, 0.5, 0.0025}];

Then do a polynomial fit on these

fittedfunc[s_] = Fit[listvals, Table[s^k, {k, 0, 130, 1}], s];

and then solve for when

FindRoot[fittedfunc[s] == 0.8, {s, 0.4, 0.5}]
{s -> 0.414704}

I was just wondering though if there is perhaps some other, more direct way to do this, without having to approximate the curve first? I had a try with ?NumericQ but couldn't quite get it work. Anybody any ideas?

$\endgroup$
4
$\begingroup$

Interpolation can also be used and I present this just this somewhat 'over the top' to illustrate a number of ways to check, hone in initial values etc.

I have 'simplified' the recursive function.

fun[s_, x_] := x/(2 x - x^2 + (1 - x)^2/(s + 1))
r[s_, n_] := Nest[FullSimplify[fun[s, #]] &, 0.00001, n]
tab = Table[{j, r[j, 50]}, {j, 0, 0.5, 0.001}];
if = Interpolation[tab];
ans = t /. First@NSolve[if[t] == 0.8, t];
p3d = ListPointPlot3D[
   Flatten[Table[{j, k, r[j, k]}, {j, 0, 0.5, 0.01}, {k, 0, 100, 1}], 
    1]];
Manipulate[
 Column[{DiscretePlot[r[s, x], {x, 0, 100}, Frame -> True, 
    PlotRange -> {0, 1}, GridLines -> {{50}, {0.8}}],
   Plot[if[x], {x, 0, 0.5}, 
    Epilog -> {Point[{s, r[s, 50]}], 
      Text[Style["x", Red], {ans, r[ans, 50]}]}, 
    GridLines -> {{ans}, {0.8}}, Frame -> True],
   Show[p3d, 
    Graphics3D[{Red, PointSize[0.04], Point[{ans, 50, r[ans, 50]}]}], 
    ListPointPlot3D[Table[{j, 50, r[j, 50]}, {j, 0, 0.5, 0.01}], 
     PlotStyle -> Black],
    ListPointPlot3D[Table[{ans, j, r[ans, j]}, {j, 0, 100, 1}], 
     PlotStyle -> Black],
    Graphics3D[{Yellow, PointSize[0.04], Point[{s, 50, r[s, 50]}]}]
    ], Row[{"Answer:", ans}]
   }], {s, 0, 0.5, Appearance -> "Labeled"}]

enter image description here

$\endgroup$
  • $\begingroup$ Many thanks for this - that's also a very elegant approach! $\endgroup$ – Tom Wenseleers Oct 5 '15 at 19:35

Your Answer

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

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