This page describes a type of fractal derived from the Newton-Raphson method, which is more normally used as an approximate method of solving equations.
This is not a new idea to me; I was given the idea by a colleague at work, and several other people have web pages about it too. I'm putting up yet another one because it contains some additions to the concept which I haven't seen anywhere else.
Newton-Raphson iteration should be familiar to anyone who has studied calculus; it's a method for finding roots of a function by using the derivative of the function to improve an approximation to the root.
To perform a Newton-Raphson approximation, suppose you have a
f(x), with derivative
you have an approximation
a to a root of the function.
The Newton-Raphson procedure is to calculate
a' = a - f(a)/f'(a), which is a
closer approximation to the root. Typically you would then iterate
this again, and again, until the successive values were
extremely close together, at which point you would conclude that you
had a very good approximation to the actual value
f(r) equals zero.
The Newton-Raphson method is useful in practice because of its extremely fast convergence. The distance from the root to each approximation is roughly squared at each iteration; so assuming the distance is already small enough that this makes it smaller rather than larger, you expect to double the number of correct decimal digits in each approximation. So if you can find a reasonably good approximation to begin with, Newton-Raphson can very quickly give you an excellent one.
What the Newton-Raphson formula is essentially doing is drawing a tangent to the curve at the point of the original approximation, then following that tangent to where it crosses the x-axis. Since any differentiable function looks close to a straight line when viewed at sufficient magnification, this explains why it works so well: the curve itself does not diverge from the tangent line by very much, and the points where they cross the x-axis are very close to each other. Hence, this technique massively improves an already good approximation.
However, if you start with a really bad approximation, much
more interesting things happen. Suppose the function curves up from
one intersection with the x-axis and back down to another, like a
parabola; and suppose your initial approximation is somewhere near
the top of this arc. Now drawing a tangent to the curve and
following it to the x-axis will land you a huge distance away from
the roots of the function - and as your initial approximation
a crosses the maximum point of the curve, the
Newton-Raphson second approximation
a' will flip from
one side of the roots to the other. In fact, as
the relatively short distance across the maximum of the curve,
a' will cover most of the real line.
This sort of behaviour, expanding a small area into a large one, is exactly the sort of behaviour we expect to give rise to self-similar fractals. So if we were to start a Newton-Raphson iteration at each point on the real line, run each iteration until it converged to within a given tolerance level of a root, and then colour the starting point according to which root it ended up at, we might well expect to see fractal shapes.
Fractals on one line are not very interesting, however; so let's work in the complex plane. The Newton-Raphson iteration still works perfectly well there, so that's where I'll be generating my fractals.
Here's an example fractal, generated from the
z^4-1 - so the four roots of the function
are at -1, +1,
In this image, we see a large boring area surrounding each root of the function - as we would expect, since any point near a root will converge rapidly to that root and do nothing interesting. But between the areas of boring well-behaved convergence, we see some beautiful fractal shapes. Let's zoom in on one of those boundary areas:
Just as we predicted - each of the heart-shaped blobs making up the boundary line is itself composed of boundary lines made up of further heart-shaped blobs. This pattern is a true fractal.
Here's a rather different example. This time the function being
(z-3)(z-2)(z-1)z(z+1)(z+2)(z+3), so it has
seven roots strung out in a long line:
In this case, the fractal shapes are much smaller compared to the overall structure of the image. But they're not completely absent. If we zoom in on a couple of the little blobs on the boundary lines, we see this:
Each blob is divided up into coloured areas similar to those covering the whole plane, and on each dividing line we see more blobs looking much the same as the larger blobs.
These images are reasonably pretty, but they're a bit garish to be turned into desktop wallpaper in their current form. Is there anything we can do to make them less stark?
Yes, there is. One obvious thing we can do, as well as noticing which root of the function the iteration ended up at, is to count how many iterations it took to get there. We can then colour each pixel a different shade of the colour assigned to that root depending on the number of iterations. So, using the obvious approach of setting the pixel shade to the number of iterations modulo the number of available shades (so that each colour cycles through those shades), we see something like this:
This is not only prettier, but it also shows us exactly where each root of the function is - instead of just knowing the roots are somewhere in the large coloured areas, we can now positively identify each root as the centre of the bright spot in each area.
The cyclic behaviour is not quite optimal, though; it works well enough if the number of available colours is limited, but it means there are sudden edges (like the ones at the very centre of each region) where a dark colour suddenly becomes a bright colour again. Perhaps if we have true colour available, it would be better to have the pixel shading be monotonic - always getting darker the more iterations are needed, but fading out by less and less and never actually reaching blackness:
Now that's starting to look much nicer, I think. But it
would be even better if the visible boundaries between different
shades of the same colour could be removed. I suspect that doing
this rigorously requires some really horrible maths and a lot of
special cases, but I've found that a good ad-hoc approximation is
obtained simply by looking at the last iteration, in which the point
first comes within the specified distance of a root. We look at
D0 from the previous point to the root and
D1 from the new point to the root, and we
know that the threshold distance
T is somewhere in
between the two. I've found that simply looking at
(log T-log D0)/(log D1-log D0), in other words whether
the log of the threshold radius was near to the start or the end of
the inward distance travelled by the point (on a logarithmic scale),
produces a perfectly acceptable result which we can use to smooth
out those boundaries:
In order to write a program to generate these images, it's necessary
to know both the function being used (typically a polynomial) and
the exact locations of all its roots. Finding the exact roots of a
general polynomial is not easy (cubics and quartics are just about
solvable, but quintics and beyond fall foul of Galois theory), so it
makes much more sense to start by deciding where we want
the roots to be, and using that to compute the polynomial by
multiplying together a series of
(x-a) terms. This does
not restrict the range of polynomials we can end up with, since in
the complex plane any polynomial can be fully factorised.
So the actual parameters you would pass to the fractal program consist of the coordinates of a set of points, together with a colour for each point. This led to an interesting idea: suppose we imagine a small number of coloured points drifting gently around the plane, and at each instant of time we compute a Newton-Raphson fractal for the current positions of the points. This should lead to a sequence of fractals which flow naturally on from each other, and as the bright central point in each coloured region moves, the regions move with them and the fractal phenomena on the region boundaries swirl continuously.
Here is such an animation. To create this, I've set up three points,
each moving along the Lissajous curve
(sin t, sin 2t),
and each one third of the way further around the curve than the
[download cascade.mpeg, 320x256 MPEG animation, 633KB]
The image quality isn't amazing (due to the MPEG compression), but one feature that's just about visible is the lines of additional bright spots moving within each coloured region, which appear to become the blobs on the connecting line when the other two regions come together to squash the line of spots.
Looking at the above rainbow-coloured plot with seven roots in a line, I found myself thinking that it's not very pretty, because the fractal blobs are so small. Why are they so small? Could anything be done to make them bigger?
I speculated that perhaps the reason the blobs are so small
(implying that Newton-Raphson converges particularly well for this
function) might be because polynomials of a high degree are
generally very steep, and curve over very sharply at the top of
peaks, so there isn't a large region of the plane in which the
fractal self-similarity can be observed. This suggested than an
obvious way to soften the process might be to reduce the overall
degree of the function: to make it more like
x^7. How to do that while preserving the
locations of the roots? Why, take the square root, of course!
Thinking about this a bit further, that does sound like a good idea. In the real numbers, the square root transformation squashes things towards the x-axis, and it squashes them more the further away they are; so I would indeed expect it to turn a very sharp switchback into a more gentle curve.
Moreover, it's trivially easy to do in the Newton-Raphson iteration.
Suppose we have a function
f, and we define a function
g to be
f^k for some power
g' = k f^(k-1) f'; so when we compute
f^(k-1) term cancels out, and we
are left with simply
(1/k) . (f/f'). In other words,
the sole effect of raising the entire function to the power
k is to multiply a constant factor of
into the distance moved per iteration.
So this had to be worth a try. Here are the results from taking the rainbow plot above, and raising the entire function to a variety of overall powers:
As I had hoped, the small fractal blobs have expanded into larger fractal blobs, which is good; and they've become more complex and spiky in shape as well, because the smaller fractal blobs on them have also expanded. What I wasn't expecting, however, was that the dividing lines between the main convergence regions have entirely changed shape, becoming more wiggly.
It turns out that actually taking the square root of the function
(raising it to the power 0.5) is not feasible, since the
algorithm converges more and more slowly the closer to 0.5 you get.
This is because raising the function to the power 0.5 causes us to
double the distance moved per iteration; so once we're near to a
f/f' gives you an accurate estimate of the
distance to the root, and we go twice as far, causing us to
end up just as far away from the root on the far side, so no wonder
the iteration fails to converge. 0.52 was the smallest value I could
use without having to increase my iteration limit.
This tendency to overshoot the root we're aiming at when using
k<1 also offers an alternative intuitive explanation
for the fractal shapes becoming bigger. The fractal shapes arise
because near to a boundary between convergence regions there is a
tendency for a single iteration step to travel a long distance into
remote regions of the complex plane; so if we're moving even further
in a single iteration, then we would indeed expect not to have to be
so close to the boundary before experiencing this phenomenon.
Raising the value of
k to something greater
than 1 has the opposite effect: we now converge toward our root in
smaller steps, more slowly but more surely, and so we're less likely
to overshoot and the fractal effects on the boundary lines become
less pronounced. I'm not going to exhibit any pictures of that,
because they're boring. However, I'll come back to this later.
As well as raising the entire polynomial to a power, it's
also interesting to see what happens if we try raising only one of
(z-a) factors to a power. To do this we must first
think about whether this can be efficiently implemented in software.
It can, it turns out. Observe that if you have a function made up of
the product of a lot of smaller factors
f = abcd, then the product rule gives its derivative as
the sum of a list of terms
f' = a'bcd + ab'cd + abc'd +
abcd', in which each term looks very similar to
itself. In fact we can divide both sides by
f to give
f'/f = a'/a + b'/b + c'/c + d'/d, which is precisely
the thing whose reciprocal we subtract from
z in the
N-R formula. Now when each factor
a is of the form
(z-k), then the resulting term
1/(z-k); and as we observed in the previous
section, if the factor is raised to a power so that it's of the form
(z-k)^e then the resulting term
only altered by a constant, so it becomes
So if we have a function which is the product of a number of terms
of the form
(z-a1)^b1, then we can efficiently perform
its Newton-Raphson iteration without ever computing the entire
function itself. We simply compute
b1/(z-a1) for each
term, sum them, take the reciprocal, and subtract that from
A useful point about this procedure is that it means we can conveniently raise each factor of our "polynomial" (which isn't very polynomial any more, really) to not only an arbitrary real power, but to a complex power if we so wish. And it turns out that we do so wish, because some very impressive fractals turn up if we do. Here's a map of 121 small fractal images, all generated from functions with the same three roots, namely the complex roots of 1. But the red root (1 itself) has been raised to a different power in each picture: across the map the real part of the power runs from 0 to 1, and downwards the imaginary part runs from 0 to 1.
As mentioned above, raising a root to a power of 0.5 or less inhibits convergence of the iteration to that root at all. But in the presence of other roots to which we can still converge, the region of non-convergence - shown in black on the above map - forms complex and interesting fractal shapes. Meanwhile, applying an imaginary part to the power of the red root causes a twisting effect: the higher the power, the more the shape is less straight and more spirally; and it appears that it's the real part of the power which must be greater than 0.5 to manage to converge.
Some of the images in that map are well worth expanding to a larger size to admire in more detail. Here are four particularly good ones:
If you look again at the images in the above section, you may notice that some of the shapes - particularly the ones that occur when things are raised to real powers near 0.5 - are beginning to look interestingly like the shape of the Mandelbrot set. This suggested to me that there might be some sort of close relationship between these fractals and Mandelbrot/Julia sets.
So an obvious question to ask along these lines is, is there a
function to which we could apply the Newton-Raphson formula and end
up with the Julia set iteration
z := z^2+c?
Well, if you write down the equation
z - f(z)/f'(z) =
z^2+c, or equivalently
z - f/(df/dz) = z^2+c,
then it looks a lot like a differential equation. If only we were
working in the positive reals instead of the complex numbers, we
could solve it quite easily by separating variables into
= -dz/(z^2-z+c), factorising the denominator of the RHS and
turning it into partial fractions, and integrating to get
(z-a)^(1/(b-a)) . (z-b)^(1/(a-b)), where
b are the roots of the quadratic
c=1/4, which is a special case arising from
z^2-z+c being a perfect square and looks rather
different; but I'll ignore that, because it's not a very interesting
Julia set anyway.)
And it turns out that even in the complex numbers this answer works
plausibly. The function
(z-b)^(1/(a-b)) is of precisely the form discussed in the
previous section: it's the product of linear terms raised to
arbitrary powers. Hence we can conveniently do the Newton-Raphson
iteration for this function by the method described above: compute
1/((b-a)(z-a)) + 1/((a-b)(z-b)) and subtract its
z. If you do that, remembering that by
ab=c, you will
indeed find that a lot of algebraic mess cancels out and you end up
with the iteration
z := z^2+c.
So there is a family of functions to which the application of the Newton-Raphson formula yields a Julia set iteration, and moreover those functions aren't very far away from the ones we've already been considering here. To prove it works, here are a couple of plots of Julia sets, with their corresponding Newton-Raphson plot:
|Julia set for
||Julia set for
|Newton-Raphson fractal for
||Newton-Raphson fractal for
You will have noticed, of course, that the colouring is very different. The two types of fractal are using matching iteration formulae, but a Julia set plotter concentrates on how long it takes the iteration value to land outside a critical circle, whereas the Newton-Raphson plotter is actually waiting for the iterates to converge to a point, and is only incidentally observing what happens when they don't. So you wouldn't actually want to throw away your dedicated Julia-set plotting programs; but it's interesting, nonetheless, that something very like Julia sets are a special case of Newton-Raphson fractals.
Another noticeable thing about some of the above fractals is that
some have much more fractal content than others. The original
fractal at the top of this page with roots at -1, +1,
+i has four big regions of flat colour
meeting at the origin, and the self-similar nature of the fractal
causes that quadruple contact point to be replicated in other parts
of the plane. This gives rise to qualitatively more interesting
fractal phenomena than the plot with seven roots in a line, in which
all the large convergence regions are separated by boring curves
which never meet, and the fractal blobs reflect this structure.
So I wondered, is there a way that we can predict how the convergence regions are going to be shaped, and thereby construct polynomials which have triple or quadruple meeting points exactly where we want them?
Yes, as it turns out, there is.
I observed above that raising a polynomial to a real power greater than 1 has the effect of causing the iteration to take smaller steps towards the root, which in turn causes convergence to be slower but surer and reduces the incidence of overshoot leading to fractal phenomena. As a result of this, the boundary lines between regions become clearer and smoother and simpler. This seemed like a good thing if we wanted to know the overall structure of the fractal; but merely clearer isn't good enough. I wanted clearest.
So I wondered what would happen "in the limit", as the iteration
speed slowed down more and more. In other words, I was interested in
what I'm going to describe as continuous Newton-Raphson,
in which you start your "iteration" at an arbitrary point
z0 on the plane and then let
continuously according to the differential equation
dz/dt = -f(z)/f'(z).
As in the previous section, we can "solve" this equation by
separating variables and naïvely integrating, trying hard to
ignore the question of whether this is rigorous or even meaningful
in the complex numbers. This time we end up with
1/f df/dz dz
= -dt, and then we can apply the integration-by-substitution
formula to the LHS to give us
df/f = -dt, which
integrates to give us
f(z) = Ae^-t for some constant
A. We can check by differentiating with respect to
t that this does indeed turn out to be a plausible
answer to the equation even in the complex numbers: the chain rule
d/dt f(z) = f'(z) dz/dt, so we end up with
f'(z) dz/dt = -Ae^-t = -f(z) as required. And the
A is clearly equal to
f at the point where we started our integration.
"Very nice", I hear you protest, "but what does that mean?"
t is a positive real, it means that whatever
z follows from our starting point
must have the property that at all times
f(z) is equal
k f(z0), where
k = Ae^-t is a real
value which continuously and monotonically decreases from 1 toward
0. So we normally expect the limit of such an "iteration" to be a
point at which
f(z) is actually equal to zero, i.e. a
But that's not the only thing that can happen. It's also possible to
imagine that we might encounter a point along this path at which
there isn't a clear direction we should head in: either there is no
direction in which we can head from
z which will cause
f(z) to continue decreasing as a real multiple of its
initial value, or perhaps there is more than one such direction.
How do we find such points? Well, for most points in the complex
plane the differential equation we started with will usually give us
a unique direction in which we can head to linearly decrease
f(z): we compute
f(z)/f'(z), and head in
the direction pointed to by that complex number. This fails if
f(z)=0, in which case there is no clear direction to
head in; we expect that, of course, because if
then we've already reached a root. But it can also fail if
f'(z)=0. So this suggests that the roots of the
f might be worth investigating.
So consider some root
f', at which
f itself is non-zero. If we can find any continuous
paths in the complex plane which start at
r and have
f increasing as a real multiple of
f(r), then any point on one of those paths will be a
point at which starting a continuous Newton-Raphson process will
cause it to head backwards along the same path and terminate at
r rather than at a root. In other words, we expect
those paths to be precisely the boundaries between the
convergence regions for the various roots (since that's the obvious
set of points which we expect not to converge sensibly to a root);
and moreover, on every such boundary we expect to find a
f' (because a continuous Newton-Raphson process
started on any boundary has to end up somewhere).
That's a lot of maths to endure without a break for a pretty picture, and it's also a long and rather handwavey chain of reasoning to endure without some sort of reassurance that what I'm saying still makes sense. So, here I present a sample polynomial Newton-Raphson fractal, with the roots of the polynomial's derivative marked as black blobs. Observe that each dividing line between convergence regions has a blob somewhere on it, and that in particular the point where three regions (and three such lines) meet has a blob.
So now we know, at least in theory, how to find the lines dividing the different convergence regions. Now, what distinguishes a simple dividing line, separating only two regions, from a point at which three or more regions meet?
Let's consider our root
We're now interested in how many directions we can head
r in, such that
f increases as a
real multiple of
f(r). Equivalently, we're interested
in directions we can head in such that
for some real positive
k; in other words, we want small
e such that
f(r+e)/f(r)-1 is real
So now let's consider the Taylor series expansion for the function
f(r+e) is equal to
f(r) + ef'(r) + e^2 f''(r) + e^3 f'''(r) + ..., and we
f'(r) = 0 by construction, so this gives us
f(r+e)/f(r)-1 = e^2 f''(r)/f(r) + e^3 f'''(r)/f(r) +
In the usual case,
f''(r) will be non-zero, so for
small values of
f(r+e)/f(r)-1 will be
approximately equal to
Ke^2 for some (complex) constant
K. This means that we expect to have two
directions in which we can head in order to make
f(r+e)/f(r)-1 real and positive: one direction
corresponding to each square root of
1/K. This fits
with what we expect, because in the usual case a root of
f' gives rise to a region boundary extending away from
it in two opposite directions.
f''(r) is zero as well, then the
e^3 term will now be the dominating one in the Taylor
expansion; so for small values of
f(r+e)/f(r)-1 will be approximately equal to
Ke^3 for some
K. And now we expect to have
three region boundaries coming out of our root
r, one for each cube root of
1/K. If the
third derivative of
f is zero at
then we can have four region boundaries, and so on.
And there it is: that's the result I've been working towards for
this entire section. A multiple-region meeting point occurs when a
f' is also a root of
f'', and the more higher derivatives are zero at
r the more regions meet there. For polynomials in
particular, this translates into
r being a repeated
f', and the more repeated the merrier.
To demonstrate beyond reasonable doubt that this actually works, I
will now construct from first principles a Newton-Raphson fractal
containing two points at each of which five regions meet.
For a five-way meeting point we need
f' to have a
four-times-repeated root; so let's set
f' to be
(z-1)^4(z+1)^4, which has two repeated roots as
desired. We "integrate" this in the naïve way (actually what
we're doing is finding an anti-derivative of it, integration in the
complex plane being a generally messy concept) to obtain the
(35z^9 - 180z^7 + 378z^5 - 420z^3 +
315z)/315. We can discard the constant factor of
1/315 since it makes no difference to the convergence,
and this gives us
35z^9 - 180z^7 + 378z^5 - 420z^3 +
315z. We find the roots of this (using Newton-Raphson for its
more conventional purpose) and feed them into our fractal plotter,
and get this picture:
And, exactly as we asked for, this has two five-way meeting points at locations +1 and -1. Bingo!
Of course, integration lets us add an arbitrary constant term
without affecting the derivative of our result. Here's what we get
if we add a constant term of
50+200i to the above
The roots of the polynomial have moved around, and the picture is distorted, but the region meeting points are still exactly where we asked for them.
Here's another curiosity. The following fractal is plotted using a
f, placing the roots at
The interesting feature of this picture is the black areas, which are regions in which Newton-Raphson failed to converge to a root within the program's limit of 256 iterations. At first sight one might assume that this is an artifact of having a finite iteration limit at all, but one would (it turns out) be wrong. Even if you crank up the iteration limit to a much larger number, those black areas stay black, because they represent genuine non-convergence.
What's actually happening here is that for this particular
polynomial the Newton-Raphson method gives rise to period-2 cyclic
behaviour. There's a pair of points
b, one in the middle of each of the two main black
areas, which have the property that a single iteration of
Newton-Raphson starting from
a takes you to
b, and vice versa. But this is more than just a cycle
between those two points: it's an attracting cycle. If you
start from a point somewhere near
iteration of Newton-Raphson will take you to somewhere even
b, from which another iteration will
land you nearer still to
a again. Hence, there's a
sizeable region around each point of the cycle which all converges
to the cycle itself, and hence never settles down to a root of the
This is interesting to me because it happens so rarely; I've plotted quite a lot of these fractals, including a lot of animated ones in which the points wander continuously around the plane, and this particular arrangement of roots is the only case I've found in which a polynomial gives rise to a non-converging area. (It's not the only actual polynomial: you can obviously rotate or translate the root positions to obtain other polynomials with the same behaviour, and you can also move the roots around by a small amount relative to each other before the black areas go away. But this general pattern - four roots in a diamond layout of roughly this shape - is the only pattern I know of that does this.)
On one level, it's reasonably easy to see "why" it happens.
If you define
g(x) = x - f(x)/f'(x) to be the
Newton-Raphson iteration for this polynomial, then find roots of the
g(g(x))=x (probably using a computer algebra
package, since it's a pretty ugly mess) you will observe that it has
roots which are not also roots of the similar first-order equation
g(x)=x (and hence of
f). That proves that
there are cycles at all; then, to prove there are
attracting cycles, compute the derivative of
G(x)=g(g(x)) at those roots and find the ones where it has
modulus less than 1 (since for small
G(r+h) will be approximately
G(r) = r and
|G'(r)| < 1 then
this will be closer to
r than we started out).
But that doesn't really explain the phenomenon, or predict
what other sorts of polynomial will give rise to cyclic behaviour of
this type. Do there exist polynomials which exhibit
n cycles for all
n > 1, for
example? How often does this happen? How would one go about
constructing polynomials with non-converging regions to order?
I don't currently know the answers to these questions, but I'd be interested to hear them if anyone else does.
If you want to generate some of these fractal images yourself, you can download a program to generate them here: C source code and precompiled Windows executable. The program will produce PNG, Windows .BMP, or PPM files, and image processing software should be able to convert those to other formats if you prefer.
Just typing "
newton" should give some help about what
all the command-line options do. If you want a quick start, here's a
selection of sample command lines you might like to try:
newton -o simple.png -s 256x256 -x 2 -c 1,0,0:1,1,0:0,0.7,0:0,0.5,1 -- -1 +1 -i +i
z^4-1. To zoom in on one of the boundary lines as I did above, replace "
-x 2" with "
-x 0.4 -X 1.05 -Y 1.05".
newton -o rainbow.png -s 320x256 -x 5 -c 1,0,0:1,0.7,0:1,1,0:0,1,0:0,0.7,1:0,0,1:0.5,0,1 -- -3 -2 -1 0 1 2 3
newton -o powermap.png -s 256x256 -x 2 -c 1,0,0:1,1,0:0,1,0 -- 1/power -0.5-0.8660254i -0.5+0.8660254i
0.5to get the Mandelbrot-like picture, and with
0.4+0.9ito get the disconnected yellow-and-green plot.
newton -o holes.png -s 320x256 -y 2 -c 1,0,0:1,1,0:0,0.7,0:0,0.5,1 -- -2.3 +2.3 -i +i
-C noto turn off cyclic shading (so the colours get uniformly darker as more iterations are needed),
-B yesto turn on blurring of the iteration boundaries, and
-f 64if you want to increase the number of shades used.
-poption, for example
-p 0.52. This will slow down convergence, so you will probably also need to increase the number of shades of colour (the
-foption, as discussed above) in order to make the result look nice.
2+i/1-0.5ispecifies a factor of
newtonon its own to list them.
Here are some pre-generated larger versions of the above images:
[bigsimple] [bigzoomed] [bigrainbow] [bigrainbowz1] [bigrainbowz2]
Newton-Raphson fractals are not a new idea of mine, although most other pages I've seen don't go so far into the maths. Here are a few such pages; you can probably find more by googling for "Newton-Raphson fractal".