# Defining a function to determine when convergence is guaranteed

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) &]]

-
What is the nature of your function? Polynomial? Rational? Transcendental? Numerically defined? If you're just looking for fixed points, is there some reason you don't just try to solve the equation $f(x)=x$? If you know the fixed point and you're really just interested in the rate of convergence, as your question suggests, then that is governed by the value of the derivative at the fixed point. Are you really just interested in fixed points, or more general periodic behavior? All these questions, and more, can be answered with more details. – Mark McClure Mar 14 '13 at 2:02
@MarkMcClure I have about 15 different functions that I need to investigate for a Chaos Theory class. Most of them are polynomial, but there are a few transcendental ones in the mix. I am basically trying to investigate how many iteration steps it takes to reach the fixed point. I am trying to define a procedure that takes as its input a function and returns the number of iterations necessary to reach the fixed point (or "close enough") to the fixed point. – Winterflower Mar 14 '13 at 18:56
@Winterflower You can simply use my solution to compute the steps necessary to reach "close enough". You needn't do it twice and select the elements (unless you have that for a different purpose). The 10 in the output in my answer is the number of steps for that function and seed (there may be a ± 1 involved). – R. M. Mar 14 '13 at 19:05
I edited your title, as "rate of convergence" is really something else. – Mark McClure Mar 14 '13 at 20:05

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.

-
Thank you for your reply! The missing "&" was a typo. Some troubleshooting showed that the problem was in the fact that FixedPoint never finished, ie. my test function testfun[x_]:=x^2+0.25 never actually reaches the fixed point. – Winterflower Mar 14 '13 at 18:44

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).

-
Super! SameTest fixed the problem. Thanks! – Winterflower Mar 14 '13 at 18:47
If floating-point values are involved, this condition should read Abs[#2 - #1] <= 3 + eps Abs[#2 - #1], where eps is some small value as appropriate for the working precision. – Oleksandr R. Mar 14 '13 at 20:15

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.

-