Bairstow's method, Durand-Kerner method and other root-finding methods

  • 103 Replies
  • 2039 Views

0 Members and 1 Guest are viewing this topic.

Offline FractalAlex

  • *
  • Fractal Flamingo
  • ****
  • Posts: 344
  • Experienced root-finding method expert
« 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/fractal-mathematics-and-new-theories/28/bairstows-method-durand-kerner-method-and-other-root-finding-methods/3499/
« Last Edit: May 20, 2020, 06:11:18 PM by FractalAlex »
"I am lightning, the rain transformed."
- Raiden, Metal Gear Solid 4: Guns of the Patriots

Offline FractalAlex

  • *
  • Fractal Flamingo
  • ****
  • Posts: 344
  • Experienced root-finding method expert
« 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 »

Offline superheal

  • *
  • Strange Attractor
  • ******
  • Posts: 92
« 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.

Offline FractalAlex

  • *
  • Fractal Flamingo
  • ****
  • Posts: 344
  • Experienced root-finding method expert
« 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! :thumbs: 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 »

Offline FractalAlex

  • *
  • Fractal Flamingo
  • ****
  • Posts: 344
  • Experienced root-finding method expert
« 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?

Offline superheal

  • *
  • Strange Attractor
  • ******
  • Posts: 92
« 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 :P
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[0] = -1;
        a[1] = 0;
        a[2] = 0;
        a[3] = 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[1] - u * b[0] - v * b[1];
        double d = a[0] - v * b[0];

        double g = b[1] - u * f[0] - v * f[1];
        double h = b[0] - v * f[0];

        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.

Offline superheal

  • *
  • Strange Attractor
  • ******
  • Posts: 92
« 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 :P

Offline FractalAlex

  • *
  • Fractal Flamingo
  • ****
  • Posts: 344
  • Experienced root-finding method expert
« 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.

Offline FractalAlex

  • *
  • Fractal Flamingo
  • ****
  • Posts: 344
  • Experienced root-finding method expert
« 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.

Offline superheal

  • *
  • Strange Attractor
  • ******
  • Posts: 92
« 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[0] = new Complex(pixel);//z

        for (int i = 1; i < degree; i++) {
            complex[i] = complex[0].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[0] = pixel
complex[1] = pixel * a
complex[2] = 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.

Offline superheal

  • *
  • Strange Attractor
  • ******
  • Posts: 92
« 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.

Offline FractalAlex

  • *
  • Fractal Flamingo
  • ****
  • Posts: 344
  • Experienced root-finding method expert
« 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? ???

Offline superheal

  • *
  • Strange Attractor
  • ******
  • Posts: 92
« 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.

Offline FractalAlex

  • *
  • Fractal Flamingo
  • ****
  • Posts: 344
  • Experienced root-finding method expert
« 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. :head_banging: 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 »

Offline superheal

  • *
  • Strange Attractor
  • ******
  • Posts: 92
« 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?