Defining a function to determine when convergence is guaranteed

Question

I'm investigating how many iterations are needed for a particular orbit to reach a fixed point of a function. Since I have many functions to test, I want to define a function in Mathematica that will take as its input

1. the function I want to test
2. the seed or starting point of the iteration
3. total number of iterations to carry out.

Since the orbit never actually reached the fixed point, I am considering any point that comes within a millionth of the actual fixed point to have reached the fixed point.

I have defined the following function

iterationsteps[function_, seed_, steps_] := Length[
 Select[FixedPointList[function, seed, steps] - FixedPoint[function, seed], # >= 1*10^(-6)]]

But whenever I try to run it on some test values, it just ends up running for ages and never returns anything. I'd be extremely grateful is someone can point me in the right direction.

Thank you!

\Edit:

Some troubleshooting and the helpful answers below revealed that the problem was caused by FixedPoint.

iterationsteps[function_, seed_, steps_] := 
Length[Select[
FixedPointList[function, seed, steps] - 
FixedPoint[function, seed, 
SameTest -> (Abs[#1 - #2] < 1*^-6 &)], # >= 1*10^(-6) &]]

Answer 1

You're missing a & in your question, which might be simply a typo.

iterationsteps[function_, seed_, steps_] := 
 Length[Select[FixedPointList[function, seed, steps] - FixedPoint[function, seed], # >= 1*10^(-6) &]]

iterationsteps[Cos, 1., 100]
(* -> 16 *)

iterationsteps[Cos, 1, 100]
(* -> runs forever *)

The second runs forever because Cos[Cos[Cos[...Cos[1]...]]] always changes when another Cos is applied to, so FixedPoint never stops. (The number 1 is treated as exact so Cos[1] stays unevaluated.) That is a common mistake, and I'm just guessing you might have made it.

Answer 2

Your solution runs forever because the FixedPointList and FixedPoint don't actually terminate (assuming what you say about your function is true), so the evaluation doesn't enter Select yet.

Instead, a simpler way would be to use the SameTest option of FixedPointList to decide the termination. As a simple example (since you didn't give a test function):

Length@FixedPointList[1 + Floor[#/2] &, 1000, SameTest -> (Abs[#2 - #1] < 3 &)]
(* 10 *)

Here I've forced termination when the difference between two successive values is less than three (instead of the default SameQ).

Answer 3

Here's a groovy trick to help find an interval about a fixed point inside of which, points are guaranteed to converge. The basic observation is based on the fact that a fixed point is (by definition) attractive if and only if the value of its derivative is less that 1 in absolute value. Thus, assuming the function is continuously differentiable, there exist a unique, largest, open interval about the fixed point with the property that $|f'(x)|<1$ for all $x$ in the interval. Once inside that interval, points must converge to the fixed point under iteration.

Let's illustrate with an example polynomial chosen from the logistic family, namely $p(x)=2.3x(1-x)$. We can define the function and classify its fixed points like so.

Clear[p];
p[x_] = 2.3 x (1 - x);
fixedPoints = x /. NSolve[p[x] == x, x]
p' /@ fixedPoints

(* Out: {0., 0.565217} *)
(* Out: {2.3, -0.3} *)

We see that the origin is a repulsive fixed point and that there is an attractive fixed point a bit bigger than $1/2$. Here's the interval for this example. Note that I cut out an annoying and unimportant warning message issued by Reduce

Reduce[Abs[p'[x]] < 1, x, Reals]

(* Out: 0.282609 < x < 0.717391 *)

I claim that if we start inside this interval, even near an endpoint, then we are guareanteed convergence. Typically, in fact, the domain of convergence will extend a bit past the endpoints.

NestList[p, 0.72, 11]

(* Out: {0.72, 0.46368, 0.571966, 0.563088, 0.565846, 0.565028, 
         0.565274, 0.5652, 0.565222, 0.565216, 0.565218, 0.565217}
*)

Of course, it's also nice to know when a point might escape to infinity. For this, you can apply the same idea to

P[x_] = Simplify[1/p[1/x]]

(* Out: (0.43478*x^2)/(-1. + x) *)

This function P is sort of a flipped version of p; the behavior of P near zero mirros that of p near infinity. Thus, the following tells us that $x<0$ or $x>1/0.4995\approx 2$, then iterates of $p$ diverge to infinity.

(* Reduce[Abs[P'[x]] < 1, x, Reals] *)

(* x < 0.449518 || x > 1.55048 *)

Thus, for this particular function, we could iterate as follows. I prefer NestWhileList as a somewhat more natural alternative to FixedPointList.

orbit[x0_] := NestWhileList[p, x0, 
  !(# < 0 || 0.29 < # < 0.71 || # > 2) &, 1, 100];
orbit[0.01]

(* Out: {0.01, 0.02277, 0.0511785, 0.111686, 0.228189, 0.405073} *)

Things are little trickier when dealing in the complex realm and functions with poles can introduce problems with infinity trick but, with care, the same basic ideas apply. Points of period $n$ can be dealt with by considering the $n^{\text{th}}$ iterate of the map.