如何估计多项式根的多重性？

Question

我要estimate multiplicity of polynomial roots.

我找到了一些info about it，选择了测试示例并制作了c程序

这里应该是 4 个根。一个单根和一个重数 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;
}

结果是：

m=1
m=1
m=1

好吃吗？也许我应该简单地添加结果？

使用符号计算的良好结果是根：[3/2，-1/2] 及其多重性：[1,3]

这是函数 f(z)= (z^2-0.75)^2-z-0.75 = z^4-1.5*z^2-z-3/16 的图

是否可以用数值计算相似的值？

Answer 1

您可以通过轮廓积分来执行此操作，请参阅here。软件是available。

Answer 2

变更摘要：

• 在循环内评估 d 之前先评估 e；
• 循环后z减去z0时，还需要d减去1才能匹配；
• 从真正的根位置扰动输入少量以避免0/0 = NaN 结果：h 必须足够小，但不能太小...

完整的程序：

#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;
}

输出：

m = 1
m = 1
m = 3

Answer 3

我在我的初始程序中发现了一个错误。寻找周期点的函数应该是

f^n(z) - z 

所以

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 

Answer 4

我选择了基于the geometrical notation of the root的方法 在The Fundamental Theorem of Algebra: A Visual Approach by Daniel J. Velleman中有描述

我数了数围绕根的圆圈有多少次颜色变化。 我使用 carg 函数返回区间 [−π; 中 z 的相位角； π]。所以计算参数的符号变化并将其除以 2。这估计了根的多重性。 它可能与上面的方法相同，但对我来说更容易理解和实施。

这是动力平面的图像

改造前：

f(z) 之后：

和代码：

// 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;
}

这是z围绕根的参数图像（它使用carg）