Re: Optimizing fixed point iteration
- To: mathgroup at smc.vnet.net
- Subject: [mg83104] Re: Optimizing fixed point iteration
- From: mcmcclur at unca.edu
- Date: Sun, 11 Nov 2007 02:56:40 -0500 (EST)
- References: <fh1ckb$csj$1@smc.vnet.net>
On Nov 9, 5:26 am, Yaroslav Bulatov <yarosla... at gmail.com> wrote:
> Can anyone see the trick to speed up the following code significantly?
> depth;step=.02;
> Table[(c = 0;
> NestWhile[-4 (# - a) (# - b) &, .51, (c++; Abs[#1 - #2] > .1) &, 2,
> depth]; c), {a, -1., 0., step}, {b, 0., 1., step}] // ArrayPlot
Here's one solution that iterates up to 50 times
over a square with resolution 1000 by 1000.
convergenceCount = Compile[{{a, _Real}, {b, _Real}},
Length[FixedPointList[-4 (# - a) (# - b) &, .51, 50,
SameTest -> ((Abs[#] > 1 + Abs[a + b] ||
Abs[#1 - #2] < 0.001) &)]]];
M = Table[convergenceCount[a, b],
{a, -1, 0, 0.001}, {b, 0, 1, 0.001}]; // Timing
{12.5899, Null}
A few comments:
* The computation was performed by V6 on my Macbook
Pro running at 2.16 GHz.
* Clearly Compile helps quite a lot. Sadly,
however, NestWhileList is not fully compilable.
Hence the use of FixedPointList with the clever
SameTest.
* The orbit diverges to rapidly infinity for most
parameter values a and b. For such parameters,
your code computes the entire orbit up to depth.
It's pretty easy to prove though that the orbit
diverges if the value z ever satisfies
|z| > max(|a*b|, 1+|a+b|).
Thus, computation can stop once this threshhold
is reached. This is particularly important since
the orbit can quickly exceed $MaxMachineNumber,
forcing high precision arithmetic.
Note that my code simply stops when the
threshhold is exceeded, returning a small
number. You might need to fiddle with it a
bit to return a large number, like your code.
The picture looks pretty good, though.
Mark