How to estimate multiplicity of the polynomial root?

Question

I want to estimate multiplicity of polynomial roots.

I have found some info about it, choosed the test example and made c program

Here should be 4 roots. One simple root and one with multiplicity 3.

#include <complex.h>
#include <math.h>
#include <stdio.h>

complex long double z0 = +1.5; // exact period = 1 stability = 3.000000000000000000 multiplicity = ?
complex long double z1 = -0.5; // exact period = 2 stability = 0.999999999999900080 multiplicity = ?
complex long double c  = -0.75; // parameter of the f function




/*
https://en.wikibooks.org/wiki/Fractals/Mathematics/Newton_method

*/
int GiveMultiplicity(const complex long double c, const complex long double z0 ,  const int pMax){

    complex long double z = z0;
    complex long double d = 1.0; /* d = first derivative with respect to z */

    complex long double e = 0.0; // second derivative with respect to z
    complex long double m;
    int multiplicity;

    int p;

    for (p=0; p < pMax; p++){

        d = 2*z*d; // f' = first derivative with respect to z */
        e = 2*(d*d +z*e); // f'' = second derivative with respect to z
        z = z*z +c ; // f = complex quadratic polynomial 
    }   

    m = (d*d)/(d*d -z*e);
    multiplicity = (int) round(cabs(m));


    return multiplicity;     
}






int main(){

    int m;


    m = GiveMultiplicity(c, z0, 1);
    printf("m = %d \n", m);

    m = GiveMultiplicity(c, z1,  1);
    printf("m = %d \n", m);

    m = GiveMultiplicity(c, z1, 2);
    printf("m = %d \n", m);




    return 0;
}

The result is :

m=1
m=1
m=1

Is it good ? Maybe I should simply add the results ?

Good results using symbolic computations are roots: [ 3/2, -1/2] and its multiplicities : [1,3]

Here is a graph of the function f(z)= (z^2-0.75)^2-z-0.75 = z^4-1.5*z^2-z-3/16

Is it possibly to compute the similar values numerically ?

I don't know how numericians solve this problem, but what about shifting the variable to bring the zero at the origin, and count the tiny low-order coefficients ? — , Mar 23 '20 at 12:42
@YvesDaoust : like here ? https://math.stackexchange.com/questions/2716268/what-is-the-intuition-for-the-multiplicity-of-a-root-of-a-polynomial-equation — Adam, Mar 23 '20 at 19:30
@Claude I have changed z after the loop : z = z - z0; Should I also change d and e ? — Adam, Apr 04 '20 at 14:38
@Adam if I update `e` before `d`, and after the loop `d -= 1` to correspond to `z -= z0`, and *also* perturb the input `z0` by a (not too big, not too small) value (e.g. `h = 1e-5`) then I get the right answers -- need to perturb to avoid 0/0 and NaN — Claude, Apr 05 '20 at 00:24
https://www.sciencedirect.com/science/article/abs/pii/S0096300315005548 — Adam, May 01 '20 at 11:42

score 2 · Answer 1 · answered Mar 22 '20 at 15:56

2

You do this with contour integration, see here. Software is available.

answered Mar 22 '20 at 15:56

user14717

4,757
2
44
68

https://www.degruyter.com/view/j/math.2018.16.issue-1/math-2018-0131/math-2018-0131.xml?lang=en – Adam Mar 22 '20 at 16:46
Nice reference. Also deals with the poles. – user14717 Mar 22 '20 at 17:06

score 2 · Answer 2 · answered Apr 15 '20 at 01:35

Summary of changes:

evaluate e before evaluating d inside the loop;
when subtracting z0 from z after the loop, you also need to subtract 1 from d to match;
perturb input a small amount from true root location to avoid 0/0 = NaN result: h must be small enough, but not too small...

Complete program:

#include <complex.h>
#include <math.h>
#include <stdio.h>

complex long double h = 1.0e-6; // perturb a little; not too big, not too small
complex long double z0 = +1.5; // exact period = 1 stability = 3.000000000000000000 multiplicity = ?
complex long double z1 = -0.5; // exact period = 2 stability = 0.999999999999900080 multiplicity = ?
complex long double c  = -0.75; // parameter of the f function

/*
https://en.wikibooks.org/wiki/Fractals/Mathematics/Newton_method
*/
int GiveMultiplicity(const complex long double c, const complex long double z0, const int pMax){
    complex long double z = z0;
    complex long double d = 1.0; /* d = first derivative with respect to z */
    complex long double e = 0.0; // second derivative with respect to z
    complex long double m;
    int multiplicity;
    int p;
    for (p=0; p < pMax; p++){
        e = 2*(d*d +z*e); // f'' = second derivative with respect to z
        d = 2*z*d; // f' = first derivative with respect to z */
        z = z*z +c ; // f = complex quadratic polynomial
    }
    d = d - 1;
    z = z - z0;

    m = (d*d)/(d*d -z*e);
    multiplicity = (int) round(cabs(m));

    return multiplicity;
}

int main(){
    int m;

    m = GiveMultiplicity(c, z0 + h, 1);
    printf("m = %d\n", m);

    m = GiveMultiplicity(c, z1 + h, 1);
    printf("m = %d\n", m);

    m = GiveMultiplicity(c, z1 + h, 2);
    printf("m = %d\n", m);

    return 0;
}

Output:

m = 1
m = 1
m = 3

Adam · Answer 3 · 2020-03-30T19:54:58.030

I have choosed the method based on the geometrical notation of the root It is described in The Fundamental Theorem of Algebra: A Visual Approach by Daniel J. Velleman

I count how many times color chages along a circle around root. I use carg function which returns the phase angle of z in the interval [−π; π]. So count the sign change of the argument and divide it by 2. This estimates the multiplicity of the root. It is probly the same method as above, but easier to understand and implement for me.

Here is the image of dynamical plane

before transformation:

and after f(z):

and the code:

// gcc p.c -Wall -lm
// ./a.out

#include <complex.h>
#include <math.h>
#include <stdio.h>


// parameter c of the function fc(z) = z^2+c is c = -0.7500000000000000 ; 0.0000000000000000 

const long double pi = 3.1415926535897932384626433832795029L;
long double EPS2 = 1e-18L*1e-18L; // 
complex double c = -0.75;
complex double z = 1.5; //-0.5;






//https://stackoverflow.com/questions/1903954/is-there-a-standard-sign-function-signum-sgn-in-c-c
int sign(long double x){
    if (x > 0.0) return 1;
    if (x < 0.0) return -1;
    return 0;

}

int DifferentSign(long double x, long double y){
    if (sign(x)!=sign(y)) return 1;
    return 0;


}




long double complex Give_z0(long double InternalAngleInTurns, long double radius )
{

  //0 <= InternalAngleInTurns <=1
  long double a = InternalAngleInTurns *2.0*pi; // from turns to radians
  long double Cx, Cy; /* C = Cx+Cy*i */


      Cx = radius*cosl(a); 
      Cy = radius*sinl(a); 


  return Cx + Cy*I;
}


int GiveMultiplicity(complex long double zr, int pMax){

    int s; // number of starting point z0
    int sMax = 5*pMax; // it should be greater then 2*pMax
    long double t= 0.0; // angle of circle around zr, measured in turns 
    long double dt = 1.0 / sMax; // t step 
    long double radius = 0.001; // radius should be smaller then minimal distance between roots 

    int p; 

    long double arg_old = 0.0;
    long double  arg_new = 0.0;

    int change = 0;
    complex long double z;
    complex long double z0;
    //complex long double zp;



    //
    for (s=0; s<sMax; ++s){
        z0 = zr + Give_z0(t, radius); // z =  point on the circle around root zr

        // compute zp = f^p(z)
        z = z0;
        for (p=0; p < pMax; ++p){z = z*z + c ;} /* complex quadratic polynomial */
        // turn (zp-z0) 
        z = z - z0; // equation for periodic_points of f for period p 

        arg_new = carg(z);  
        if (DifferentSign(arg_new, arg_old)) {change+=1;}
        arg_old = arg_new;



        //printf("z0 = %.16f %.16f  zp = %.16f %.16f\n", creal(z0), cimag(z0), creal(zp), cimag(zp));
        t += dt; // next angle using globl variable dt
    }

    return change/2;

}



int main(){

    printf("multiplicity = %d\n", GiveMultiplicity(z,2));

    return 0;
}

And here is the image of argument of z around root ( it uses carg )

Can you do this faster than `O(d^2*M(n))`, d the degree of the polynomial, n the working precision of the number type? Else expanding the polynomial around the root and arguing about coefficient sizes will be faster. Or the `O(d*log(d)*M(n))` method in the first answer that computes the winding number directly using FFT. — Lutz Lehmann, Apr 15 '20 at 06:42

score 0 · Answer 4 · answered Mar 30 '20 at 17:13

0

I have found one error im my initial program. Function for finding periodic points should be

f^n(z) - z

so

for (p=0; p < pMax; p++){

        d = 2*z*d; // f' = first derivative with respect to z */
        e = 2*(d*d +z*e); // f'' = second derivative with respect to z
        z = z*z +c ; // f = complex quadratic polynomial 
    }
z = z - z0; // new line

answered Mar 30 '20 at 17:13

Adam

1,254
12
25

1

then you should also do `d = d - 1;` I suppose – Claude Apr 04 '20 at 23:56

How to estimate multiplicity of the polynomial root?

4 Answers4