Multiple Critical Points

[jdebord]

I read "1+abx" in both the book and the web page. There should be no problem

[element90]

I read "1+abx"

The + wasn't as well formed as it is elsewhere in the example code, looking again with the knowledge that is a + I can see that it is. That makes much more sense.

When I've got time I'll have a go at finding the roots of a polynomial using the example code as I want to allow the use of complex coefficients.

Thanks for help so far.

[element90]

I've now got laguere implemented as LaguerreMethod but I keep on getting "too many iterations". I changed the throw to just output the message to console so that I can continue.

For

z^5 + z^4 + z^3 + z^2 + z + 1 = 0

Code:

========
Laguerre
========
(-0.5,0.866025403784439)
(-0.5,-0.866025403784439)
(-0.621744414124811,0.440596998965663)
(0.5,0.866025403784439)
2.83554605884655e-12
too many iterations in LaguerreMethod
(0.12174441412198,-1.30662240275026)
----------
(1.11022302462516e-16,-1.66533453693773e-16)
(2.22044604925031e-16,-1.66533453693773e-16)
(0,5.55111512312578e-17)
(0,0)
(0,0)
----------
(0.5,-0.866025403784439)
(0.5,0.866025403784439)
(-0.5,0.866025403784439)
(-0.5,-0.866025403784439)
(-0.5,0.866025403784439)

The first set of five values are the unpolished values out of laguere (the number before the "too many iterations message is the absoulute difference between the final two values of the "root"), the middle five values are the polished roots plugged back into the polynomial and the final set are the polished roots. There are only 3 distinct values in the five roots, does anybody know if these are correct?

I also tried

6z^5 + 5z^4 + 4z^3 + 3z^2 + 2z + 1 = 0

this produced

Code:

Laguerre
========
(-0.37569519922526,0.570175161011412)
(-0.663279565832258,0.22112299756456)
(-0.232299697642872,-0.569492124305167)
(0.10671690809259,0.986823103537636)
2.70791049882879e-12
too many iterations in LaguerreMethod
(0.331224221271789,-1.20862913780804)
3.59210439948278e-88
too many iterations in LaguerreMethod
----------
(0.294194556360142,-0.668367097443301)
(0.294194556360142,0.668367097443301)
(-0.37569519922526,-0.570175161011412)
(-0.670332047603097,0)
(-0.37569519922526,0.570175161011412)
----------
(-2.22044604925031e-16,0)
(-2.22044604925031e-16,0)
(1.11022302462516e-16,-5.55111512312578e-17)
(0,0)
(0,5.55111512312578e-17)

Note the second "too may iterations" message, that occurred when polishing the fourth root. Note: I've moved the replacement of very small imaginary values to after the polishing of the values. Curiously the best root found is the one that produced the "too many iterations" messages.

Finally I tried

5z^4 + 4z^3 + 3z^2 + 2z + 1 = 0

Code:

========
Laguerre
========
(-0.53783227490299,0.358284686345128)
(-0.646905558112761,-0.318596356712242)
(0.14518558903178,1.08792289420352)
2.49093028900545e-12
too many iterations in LaguerreMethod
(0.239552243981517,-1.12761122383684)
----------
(0.13783227490299,-0.678154389105336)
(0.13783227490299,0.678154389105336)
(-0.53783227490299,-0.358284686345128)
(-0.53783227490299,0.358284686345128)
----------
(0,2.77555756156289e-16)
(0,-2.77555756156289e-16)
(0,-1.66533453693773e-16)
(-2.22044604925031e-16,-5.55111512312578e-17)

Again there is the "too many iterations" message, the roots found here match the roots I found using Newton-Raphson and trial and error for the initial values.

[jdebord]

I got this with my FreeBasic program:

Code:

Polynomial:

 +    1.000000
 +    1.000000 X
 +    1.000000 X^2
 +    1.000000 X^3
 +    1.000000 X^4
 +    1.000000 X^5


 1 real root(s):

Z(1) =   -1.000000

 4 complex roots:

Z(2) =    0.500000 +    0.866025
Z(3) =    0.500000 -    0.866025
Z(4) =   -0.500000 +    0.866025
Z(5) =   -0.500000 -    0.866025



Polynomial:

 +    1.000000
 +    2.000000 X
 +    3.000000 X^2
 +    4.000000 X^3
 +    5.000000 X^4
 +    6.000000 X^5


 1 real root(s):

Z(1) =   -0.670332

 4 complex roots:

Z(2) =    0.294195 +    0.668367
Z(3) =    0.294195 -    0.668367
Z(4) =   -0.375695 +    0.570175
Z(5) =   -0.375695 -    0.570175



Polynomial:

 +    1.000000
 +    2.000000 X
 +    3.000000 X^2
 +    4.000000 X^3
 +    5.000000 X^4


 4 complex roots:

Z(1) =   -0.537832 +    0.358285
Z(2) =   -0.537832 -    0.358285
Z(3) =    0.137832 +    0.678154
Z(4) =    0.137832 -    0.678154

I have checked the roots with Maple.

I use the eigenvalue method for degree > 4, otherwise I use the algebraic solution.

[element90]

It looks like I'll need to use a combination of approaches as I want to allow complex coefficients. Am I right in thinking that the eigenvalue approach only works for real coefficients? Am I also right in thinking that Laguerre can't always find real roots when real coefficients are used?

[jdebord]

The eigenvalue codes used in both "Numerical Recipes" and my own libraries are translations of the Fortran codes in the EISPACK library (http://www.netlib.org/eispack/). This library also has routines for complex matrices which could be translated too, but this may take some time...

On the other hand, if the degree of the polynomial does not exceed 5, so that the degree of the derivative does not exceed 4, the standard formulas should be used with complex numbers I think.

Laguerre's method is an iterative method. There is no guarantee that it will find all the roots.

[element90]

Thanks for the link to the EISPACK library.  As I've already discovered that Laguerre can't guarantee to find all the roots, so a method of finding roots by eigenvalue of a complex matrix is what I want, translating from Fortran to C++ won't be a problem although it's a long time since I did any Fortran programming (Fortran IV and then Fortran 77).

I'm already on a derivative of degree 6.

EISPACK has been superseded by LAPACK of which there is a C++ implementation lapack++ which has in turn been superseded by the Template Numerical Toolkit (TNT). I'll see if I can find what I want in TNT.

[element90]

EISPACK has been superseded by LAPACK of which there is a C++ implementation lapack++ which has in turn been superseded by the Template Numerical Toolkit (TNT). I'll see if I can find what I want in TNT.

I can't find any sign of routines for finding the eigenvalues of a matrix of complex numbers in TNT or lapack++ and the documentation isn't much help. I resorted to EISPACK and found a routine for finding the eigenvalues of a matrix of complex numbers. Eek! It's Fortran spaghetti in Fortran IV, even Fortran 77 would've made it easier. I've disentangled the code to produce a C++ which probably has numerous bugs in int.

[jdebord]

There is a Fortran 90 version here:

http://people.sc.fsu.edu/~jburkardt%20/f_src/eispack/eispack.html

The same site has also versions in C and C++ but they don't seem to include the complex numbers.

I have started to extend my own routines, beginning with the classical formulas for degrees < 5. I will adapt the eigenvalue code later.

By the way, which is the general formula of a polynomial of degree n, for generating fractal graphics? If I understand correctly, the coefficient of z^n should be 1 and the coefficient of z^(n-1) should be 0, so:



Is this correct?

[element90]

By the way, which is the general formula of a polynomial of degree n, for generating fractal graphics? If I understand correctly, the coefficient of z^n should be 1 and the coefficient of z^(n-1) should be 0, so:

f(z) = z^n + a_{n-2} z^{n-2} + ... + a_2 z^2 + a_1 z + c

Is this correct?

I've been using polynomial formulae with parameters (coefficients) for each z term. The only thing that was lacking was a means of determining the critical points. My program, Saturn, only went up as far as quintics which reduces to a quartic derivative for determining the critical values. To find critical points easily I've been zeroing out some of the terms so that the derivative factors out to zero and an nth root of a number or zero and a quadratic. I'm now in the process of developing a program called Neptune which will automatically determines all the critical points and uses all the available coefficients.

There is a Fortran 90 version here:

http://people.sc.fsu.edu/~jburkardt%20/f_src/eispack/eispack.html

The same site has also versions in C and C++ but they don't seem to include the complex numbers.

I've also find it difficult finding code for C++ for matrices of complex numbers. I've also had trouble with the EISPACK code I downloaded, I built it using gfortran but for some reason the matrix gets corrupted when passed into a subroutine the first row was duplicated for all the remaining rows, using constant array dimensions fixed that problem so that I can step through the Fortran and C++ versions to find out where the C++ version goes astray. I suspect the problem may be because the code I downloaded looks like Fortran IV and the compiler doesn 't refer to any standards prior to Fortran 95. I'll have a look at the Fortran 90 code, that will be easier to deal with.

There haven't been any pictures for while:


https://copy.com/9RUM6ytyXcG6

This is for the polynomial for the above picture,

z = z^6 + z^5 + z^4 + z^3 + z^2 + z + c

So the critical points are the roots of

f'(z) = 6z^5 + 4z^4 + 3z^3 + 2z + 1 = 0

To make this general I'll need to add a parameter for each coefficient in the fractal formula and to multiply each of them by the associated power so that they can then be used as the coefficients of the polynomial used to find the critical points.

I've also included in Neptune

z = sin(z) + c

and will include formulae such as

z = alpha*z^2 + beta*z^(-2) + c

There will probably be a number of such formulae as it's easier to optimise specific formulae than it is for a general formula that allows the powers to changed with the proviso that they are always real and integers and of opposite signs. The critical points for these formulae are always the nth roots of a number and are consequently easy to find.

[element90]

I have now integrated a polynomial root solver into my new program called Neptune (currently under development) which allows parameters to added as coefficients of polynomials. As the coefficients are allowed to be complex the roots are found by finding the eigenvalues of a complex matrix, the code is a translation of Fortran 90 EISPACK routines into C++.

The critical points are determined once for each fractal to be calculated as only fixed coefficient values will be allowed (for version 1.0.0), the use of c (location in the complex plane) in the coefficients will be considered for a future version at which point the critical points will have to calculated for each location.

I've also been looking at formulae of this sort:

z = alpha*z^2 + beta*z^-2 + c

the critical points for this formula are solution of

f'(z) = 2*alpha*z - 2*beta*z^-3 = 0

z^4 = beta/alpha

Happily the solutions of this formulae is easy, they are just the 4th roots which are all equally spaced on a circle centred at the origin. Using the same values for alpha and beta the critical values are at 1, -1, i and -i which produces two pictures one for 1 and -1 and one for i and -i.

Critical points 1 & -1, alpha = beta = 1


https://copy.com/Q70jbCrO8QCm

Critical points i and -i, alpha = beta =1


https://copy.com/Hh9eIjXtF16N

Combined:


https://copy.com/HZ4FSRP83cPN

alpha = beta = 0.75


https://copy.com/KJ86MbRZzgnh

alpha = beta = 0.5


https://copy.com/L9Z8xTeycRQC

alpha = beta = 0.4


https://copy.com/LR8awazooOvm

alpha = beta = 0.34


https://copy.com/cVvs16ZPO82B

[jdebord]

I have tried the formula z^p + c / z^q (with only one critical point). Here is an example with p = q = 2 :

[element90]

The formula z = z^2 + c/z^2 has four critical points which are the roots of z^4 = c, all the pictures are the same so only one critical point is required.

Modifying the formula slightly to z = z^2 + c/z^2 + c does not alter the critical points but there are now two pictures one for +/-abs(z^0.25) and one for +/-abs(z^0.25)i.


https://copy.com/QYdpuokloBoP


https://copy.com/IDocLu0sa3cH

And as two layers:


https://copy.com/sf8ncxzC1Byx

[laser blaster]

I've also been playing around with fractals with multiple critical points lately. I think the "most official" way to display these fractals is to color a point black only if it's in the set for every critical point. That's because the resulting set of points is the connectedness map for the corresponding julia sets, which is a key property of the original mandelbrot set. I don't think there's a clear-cut best algorithm for coloring these things; the method I use is to iterate all of the critical points simultaneously and use the point with the greatest distance from the origin for coloring after bailout.

The problem with the above method is that the results can be kinda ugly, or at least nowhere near as graceful as the classic M-set. The resultant sets often don't seem to be truly fractal in some spots, although you can find interesting fractal detail in other spots.

So the other method I've experimented with is coloring each point black if it's in the set for ANY critical point. It comes a lot closer to the grace and symmetry of the M-set, but it it results in clearly overlapping and intersecting features, which make me feel it's somehow not ideal, either. Something I'm  reeally curious about: if the set of points on the parameter plane where all the critical points escape under iteration is the connectedness map of the corresponding Julia sets, then does the set of points where ANY critical points escape have some special topological meaning as well? I doubt anyone would know, but it's interesting to ponder.

Anyway, I read on Wikipedia that the behavior of cubic and higher degree polynomial maps is much more complicated than that of quadratic maps, so connectedness maps for cubic functions aren't really considered useful. I guess a simple, binary measure of behavior can't cut it for cubic fractals.

So, that brings me back to what you're doing. I think you're on the right track by combining the images from each critical point. But I wish I understood all this complex analysis and topology stuff so I could figure out what's really going on behind the scenes.

[xenodreambuie]

I too have been working with multiple critical points. The two main reasons for rendering Mandelbrot sets are to zoom into interesting/beautiful regions for their own sake, and as an atlas for Julia sets. For an atlas, individual sets are often incomplete, so a proper way to combine them could be useful if there is one. I'm not optimistic, because the situation is even worse. In addition to the sets between each critical point and infinity, there are sets between pairs of finite critical points and some may be different from the infinite ones. When infinity is not a critical point there are only sets between pairs of finite critical points. Any method using more than a subset would become unworkable for more complicated formulae. So I'm inclined to just flick through the possible sets and choose one at a time to work with. Any mathematicians' insights would be welcome.