 • June 19, 2021, 10:46:40 PM

Login with username, password and session length

### Author Topic:  Bairstow's method, Durand-Kerner method and other root-finding methods  (Read 6665 times)

0 Members and 1 Guest are viewing this topic.

#### FractalAlex ##### Bairstow's method, Durand-Kerner method and other root-finding methods
« on: May 19, 2020, 10:12:45 PM »
For several weeks, on ManpWin, I tried to make Bairstow fractals, without success, although I had similar and different results.
Of course, here's the link: https://en.wikipedia.org/wiki/Bairstow%27s_method

Can someone help me? Some help would be appreciated. ManpWin uses .frm files, also used by FractInt.

Hopefully I'm the first one to start a topic about Bairstow's method, a root-finding method.

Here are some images below as reference, from Wikipedia.

The first image uses the function f(x) = x^5 - 1, the second one is f(x) = x^6 - x, and the third one is f(x) = 6x^5 + 11x^4 - 33x^3 - 33x^2 + 11x + 6.

Linkback: https://fractalforums.org/index.php?topic=3499.0
« Last Edit: May 20, 2020, 06:11:18 PM by FractalAlex »

#### FractalAlex ##### Articles and test images: Bairstow's method
« Reply #1 on: May 20, 2020, 02:02:10 PM »
I found some interesting links about Bairstow's method several weeks back.

https://pdfs.semanticscholar.org/d577/0c83a488f6e71d188892b91a6eded1f13a18.pdf

https://core.ac.uk/download/pdf/33474975.pdf

To note, this fractal is implemented in hrkalona's Fractal Zoomer Java-based application. I try to make variations of it.

And the link for Fractal Zoomer's branch on Bairstow's method is present at:
https://github.com/hrkalona/Fractal-Zoomer/tree/master/src/fractalzoomer/functions/root_finding_methods/bairstow
where you can find several .java files to help me make these fractals.

Edit: here is the attempted code that is used on my copy of ManpWIN, version 4.00e:
Code: [Select]
Bairstow3 {; Julia set of Bairstow's method applied to z^3 - 1 = 0  z = pixel  i = sqrt(-1):  u = real(z)  v = imag(z)  f = z^3 - 1  qz = f/(z^2 - u*z - v)  a3 = 1  a2 = 0  a1 = 0  a0 = -1  b3 = a3  b2 = a2 - (u*b3)  b1 = a1 - (u*b2) - (v*b3)  b0 = a0 - (u*b1) - (v*b2)  xqz = (z^2 - u*z - v)*tz + f2*(z - u) + f1  tz = 1 - z^(-3)  f3 = b3  f2 = b2 - (u*f3) - (v*f4)  f1 = b1 - (u*f2) - (v*f3)  c = a1 - (u*b0) - (v*b1)  d = a0 - (v*b0)  g = b1 - (u*f0) - (v*f1)  h = b0 - (v*f0)  ug = u*g  vg = v*g  deltav = (-b1 + (f2/f1) * b0)/((-f2^2/f1)+f3)  deltau = (-b0 - f2*deltav)/f1  u = deltau + u  v = deltav + v  z = -2*u + (u^2 + v*abs(v))*i  |z| >= 10^(-21)}
And two images for comparison
« Last Edit: May 20, 2020, 02:20:23 PM by FractalAlex »

#### superheal

• Strange Attractor
•      • Posts: 97 ##### Re: Bairstow's method
« Reply #2 on: May 20, 2020, 03:27:52 PM »
Hi, I am the developer of Fractal Zoomer. I think I added this method in the last release.
I can help you out if you need any help.
Just to point out though, Bairstow's method will not find the roots of the input polynomial.
This method only finds the coefficient u and v for a quadratic polynomial. The roots of that quadratic polynomial are a subset of the original polynomial's roots.
When those coefficients are determined, the original polynomial should be divided by the quadratic polynomial, and then the same procedure should be repeated to the quotient polynomial.

My implementation only does the first step, so it matches the images posted on that wiki.

#### FractalAlex ##### Re: Bairstow's method
« Reply #3 on: May 20, 2020, 04:02:02 PM »
You are correct : it was introduced in version 1.0.7.5. I see, at least on this point. I see that the wikibooks page hasn't been updated to include this formula and Durand-Kerner. Also, I'd like to have the Aberth-Ehrlich method as a feature request in Fractal Zoomer, despite it's complexity. I'll keep trying and keep you up to date for further development of my efforts on ManpWIN. Also, thank you for the amazing software, you did an excellent job right there! I'd love to see how these methods work. Also I discovered a page with a lot of new root-finding methods on a German page made by Jürgen Meier, I think it would be a good idea to check it out:
http://www.3d-meier.de/tut20/Seite1.html
These methods are variations of Newton's method, in section 5.
You can also find downloable code files for each of these methods, at: http://www.3d-meier.de/tut20/Download/Seite1.html

Can you explain the following code here, it would be helpful to me to have a detailed explanation:

Code: [Select]
        double u = z.getRe();        double v = z.getIm();
Code: [Select]
return z.sub_mutable(new Complex((-h * c + g * d) * (factor), (-vg * c + (ug - h) * d) * (factor)))
« Last Edit: May 20, 2020, 04:39:35 PM by FractalAlex »

#### FractalAlex ##### Re: Bairstow's method
« Reply #4 on: May 20, 2020, 04:33:01 PM »
One of my attempts yields an interesting result with this :
Code: [Select]
Bairstow3 {; Julia set of Bairstow's method applied to z^3 - 1 = 0  z = pixel  i = sqrt(-1):  u = real(z)  v = imag(z)  a3 = 1  a2 = 0  a1 = 0  a0 = -1  b3 = a3  b2 = a2 - (u*b3)  b1 = a1 - (u*b2) - (v*b3)  b0 = a0 - (u*b1) - (v*b2)  f3 = b3  f2 = b2 - (u*f3) - (v*f4)  f1 = b1 - (u*f2) - (v*f3)  c = a1 - (u*b0) - (v*b1)  d = a0 - (v*b0)  g = b1 - (u*f0) - (v*f1)  h = b0 - (v*f0)  ug = u*g  vg = v*g  factor = 1 / ((v*g) * g + h * (h - (u*g)))  u = u - (-h * c + g * d) * (factor)  v = v - (-vg * c + (u*g - h) * d) * (factor)/f1  z = -2*u + (u^2 + v*abs(v))*i  |u| >= 10^(-21)}
Just a small question: Is the Bairstow fractal a Julia set, and is it using a two-dimensional version of Newton's method?

#### superheal

• Strange Attractor
•      • Posts: 97 ##### Re: Bairstow's method
« Reply #5 on: May 20, 2020, 05:16:54 PM »
Thanks for sharing more methods and your kind words! As for the wiki is a bit outdated, its hard to find a balance between work, relaxation and working my fractal project Anyway, here is all you have to know for bairstow's method.

First of all, the initialization step:
pixel = C
z(0) = -2 * Re(pixel) + (Re(pixel)^2 + Im(pixel) * Abs(Im(pixel))) * i
In other words, z = -2s + (s^2 + t|t|)i , where s is the real part of pixel and t is the imaginary part.

You have to define 3 helper vectors, in my case a,b,f
Their lengths are n + 1, where n the the degree of the polynomial.
For instance the polynomial of z^3 - 1, has a degree of 3.
and here is the initialization for that polynomial (https://github.com/hrkalona/Fractal-Zoomer/blob/master/src/fractalzoomer/functions/root_finding_methods/bairstow/Bairstow3.java#L56)

Code: [Select]
        n = 3;        a = new double[n + 1]; //z^3 - 1        a = -1;        a = 0;        a = 0;        a = 1;        b = new double[n + 1];        b[n] = b[n - 1] = 0;        f = new double[n + 1];        f[n] = f[n - 1] = 0;
The initialization part needs to be done for every pixel.

Then you need the function that is going to be iterated (https://github.com/hrkalona/Fractal-Zoomer/blob/master/src/fractalzoomer/functions/root_finding_methods/bairstow/BairstowRootFindingMethod.java#L48):

Code: [Select]
    public Complex bairstowMethod(Complex z) {        double u = z.getRe();        double v = z.getIm();        for (int i = n - 2; i >= 0; i--) {            b[i] = a[i + 2] - u * b[i + 1] - v * b[i + 2];            f[i] = b[i + 2] - u * f[i + 1] - v * f[i + 2];        }        double c = a - u * b - v * b;        double d = a - v * b;        double g = b - u * f - v * f;        double h = b - v * f;        double ug = u * g;        double vg = v * g;        double factor = 1 / (vg * g + h * (h - ug));        return z.sub_mutable(new Complex((-h * c + g * d) * (factor), (-vg * c + (ug - h) * d) * (factor)));    }
When calling this function, the input is z(n) and the output z(n+1)
so,
initially u = Re(z_n) and v = Im(z_n)
The rest of the operations are just plain doubles, using the vectors that you already have set up before.
The last line means: z = z - ((-h * c + g * d) * factor + (-vg * c + (ug - h) * d) * factor * i)
or in general z = z - complex_number
and the real part of the complex number is = (-h * c + g * d) * factor
and the imaginary = (-vg * c + (ug - h) * d) * factor

The bailout technique is the same for all root finding methods:

norm(z - z_prev) <= epsilon

I believe this is all you need.

#### superheal

• Strange Attractor
•      • Posts: 97 ##### Re: Bairstow's method
« Reply #6 on: May 20, 2020, 05:21:26 PM »
I will add 3 more methods on the next update, but I guess there are so many of them so it might take some years to add them all #### FractalAlex ##### Re: Bairstow's method
« Reply #7 on: May 20, 2020, 05:45:51 PM »
Just to note that an .frm file works very differently than that of a .java file. It might take some time for me figure it out. An .frm file is used by FractInt and ManpWIN.

#### FractalAlex ##### Durand-Kerner method
« Reply #8 on: May 20, 2020, 06:12:44 PM »
For now, can you explain me how Durand-Kerner method works? I have also tried to make it work, but with no success on ManpWIN.

#### superheal

• Strange Attractor
•      • Posts: 97 ##### Re: Bairstow's method, Durand-Kerner method and other root-finding methods
« Reply #9 on: May 20, 2020, 06:53:32 PM »
Durand Kerner is based on this wiki page: https://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method
https://de.wikipedia.org/wiki/Datei:Weierstrass_1_m3_3_m5.png

First of all the initialization step (https://github.com/hrkalona/Fractal-Zoomer/blob/master/src/fractalzoomer/functions/root_finding_methods/durand_kerner/DurandKernerRootFindingMethod.java#L88):

Code: [Select]
       Complex[] complex = new Complex[degree];        complex = new Complex(pixel);//z        for (int i = 1; i < degree; i++) {            complex[i] = complex.times(a.pow(i));//(a^i) * z        }degree stands for the polynomial degree. You need to create an array of complex numbers an initialize them like that,
complex = pixel
complex = pixel * a
complex = pixel * a^2
...
complex = pixel * a^i

where a is an arbitrary complex constant, for example 0.4+0.9i as used in that wiki.
This array will hold all the roots of the polynomial, which means that in the end, you will have all the solutions stored.

Two more helper complex vectors are needed, to be allocated,
new_val = new Complex[degree];
fz = new Complex[degree];

now comes the iteration function (https://github.com/hrkalona/Fractal-Zoomer/blob/master/src/fractalzoomer/functions/root_finding_methods/durand_kerner/DurandKernerRootFindingMethod.java#L88):

Code: [Select]
    public Complex[] durandKernerMethod(Complex[] roots, Complex[] fz) {        for (int i = 0; i < degree; i++) {            Complex denominator = new Complex(1, 0); // denominator = 1 + 0i            for (int j = 0; j < degree; j++) {                if (j == i) { // if i equals j then the later multiplication will make the denominator equal to zero                    continue;                }                denominator.times_mutable(roots[i].sub(roots[j])); // denominator = denominator * (roots[i] - roots[j])            }            new_val[i] = roots[i].sub(fz[i].divide_mutable(denominator)); // new_val[i] = roots[i] - (fz[i] / denominator)        }        for (int i = 0; i < degree; i++) {            roots[i] = new_val[i];        }        return roots;    }
which is called like that for z^3-1 (https://github.com/hrkalona/Fractal-Zoomer/blob/master/src/fractalzoomer/functions/root_finding_methods/durand_kerner/DurandKerner3.java#L65)

Code: [Select]
@Override    protected void function(Complex[] complex) {                for(int i = 0; i < degree; i++) {            fz[i] = complex[i].cube().sub_mutable(1); // fz[i] = complex[i]^3 - 1        }                durandKernerMethod(complex, fz);            }
As a bailout technique, it is mentioned that all roots (p,q,r,s... wiki's terminology), should not differ much from their previous value, but in my current implementation I only check only the first root's convergenece.
All the root values are stored in the array which I call complex[degree],
and the check bailout check is made right here: https://github.com/hrkalona/Fractal-Zoomer/blob/master/src/fractalzoomer/functions/root_finding_methods/durand_kerner/DurandKernerRootFindingMethod.java#L106

at some point I will have to fix it so it matches the wiki's bailout.

#### superheal

• Strange Attractor
•      • Posts: 97 ##### Re: Bairstow's method, Durand-Kerner method and other root-finding methods
« Reply #10 on: May 20, 2020, 06:55:43 PM »
both of those implementations do not look like the usual root finding methods, so they are a little bit trickier to implement.

#### FractalAlex ##### Re: Bairstow's method, Durand-Kerner method and other root-finding methods
« Reply #11 on: May 20, 2020, 07:06:13 PM »
But the greatest advantage, numerically, is that it can find multiple or all the roots simultaneously. Easier when written and calculated on paper than rendered as a fractal, as you said. On the wiki page, there is a parameter $z_j$ used. What does it mean and how is it calculated? #### superheal

• Strange Attractor
•      • Posts: 97 ##### Re: Bairstow's method, Durand-Kerner method and other root-finding methods
« Reply #12 on: May 20, 2020, 07:17:11 PM »
a is constant and it can have any value you want. In the wiki they used a = 0.4+0.9i.
In my implementation I used 0.4+0.9i as well, to match the wiki, but in the custom polynomial function I let the user to pick his desired value.

The drawback of those two methods, is that they only work with polynomials. All the other root finding methods can work on any function that has roots.

#### FractalAlex ##### Re: Bairstow's method, Durand-Kerner method and other root-finding methods
« Reply #13 on: May 20, 2020, 07:32:58 PM »
Oof, it does turn out that these methods are much harder to replicate in a .frm file for ManpWIN. I'm new to these complex methods. Especially if the letter j is involved. Speaking of which, have you heard of Bisection method and Regula falsi? These are relatively simple methods. Regula falsi is a combination of the Secant and bisection methods.
« Last Edit: May 20, 2020, 07:43:53 PM by FractalAlex »

#### superheal

• Strange Attractor
•      • Posts: 97 ##### Re: Bairstow's method, Durand-Kerner method and other root-finding methods
« Reply #14 on: May 21, 2020, 09:24:42 AM »
Yes I am familiar with bisection method, I think in high school I learned the principle, that is method is based uppon (https://en.wikipedia.org/wiki/Intermediate_value_theorem)
I dont know Regula falsi, I will do a little bit of searching.
Now that I read the last comments again, you said something about Zj Can you give me the link, so I can see what it is all about?

### Similar Topics ###### Root-finding via eigenvalue/matrix methods

Started by FractalAlex on Fractal Mathematics And New Theories

4 Replies
521 Views April 25, 2021, 06:50:31 PM
by FractalAlex ###### Root finding via subdivision/IA

Started by marcm200 on Fractal Mathematics And New Theories

13 Replies
853 Views March 25, 2021, 03:03:23 PM
by marcm200 ###### Is it possible to speed up rendering with this method?
3 Replies
418 Views April 19, 2019, 12:59:29 AM
by gerrit ###### Continuous Newton's method

Started by FractalAlex on Programming

0 Replies
218 Views December 04, 2020, 01:19:40 AM
by FractalAlex ###### Method of Early Checking Whether a Point Is In the Set

Started by Byte11 on Fractal Mathematics And New Theories

12 Replies
1066 Views December 24, 2017, 12:33:10 AM
by Byte11