L'algorithme hybride de Newton Raphson ne parvient pas à une solution

Question

Brève explication du problème : j'utilise l'algorithme de Newton Raphson pour la recherche de racines dans les polynômes et il ne fonctionne pas dans certains cas. pourquoi ?

J'ai pris dans "numerical recipes in c++" un algorithme hybride Newton Raphson, qui bissecte dans le cas où New-Raph ne converge pas correctement (avec une faible valeur de dérivée ou si la vitesse de convergence n'est pas rapide).

J'ai vérifié l'algorithme avec plusieurs polynômes et il a fonctionné. Maintenant, je teste dans le logiciel que j'ai et j'ai toujours obtenu une erreur avec un polynôme spécifique. Mon problème est que je ne sais pas pourquoi ce polynôme n'arrive pas au résultat, alors que beaucoup d'autres y arrivent. Comme je veux améliorer l'algorithme pour n'importe quel polynôme, j'ai besoin de savoir lequel est la raison de la non convergence afin de pouvoir le traiter correctement.

Je publierai ci-dessous toutes les informations que je peux fournir sur l'algorithme et le polynôme dans lequel j'ai l'erreur.

Le polynôme :

f(t)= t^4 + 0,557257315256597*t^3 - 3,68254086033178*t^2 +
+ 0,139389107255627*t + 1,75823776590795

Il s'agit d'une dérivée première :

 f'(t)= 4*t^3 + 1.671771945769790*t^2 - 7.365081720663563*t + 0.139389107255627

Plot :

Racines (par Matlab) :

  -2.133112008595826          1.371976341295347          0.883715461977390 
  -0.679837109933505

Algorithme :

double rtsafe(double* coeffs, int degree, double x1, double x2,double xacc,double xacc2)
    {
    int j;
    double df,dx,dxold,f,fh,fl;
    double temp,xh,xl,rts;
    double* dcoeffs=dvector(0,degree);
    for(int i=0;i<=degree;i++)
        dcoeffs[i]=0.0;
    PolyDeriv(coeffs,dcoeffs,degree);
    evalPoly(x1,coeffs,degree,&fl);
    evalPoly(x2,coeffs,degree,&fh);
    evalPoly(x2,dcoeffs,degree-1,&df);
if ((fl > 0.0 && fh > 0.0) || (fl < 0.0 && fh < 0.0))
    nrerror("Root must be bracketed in rtsafe");

if (fl == 0.0) return x1;
if (fh == 0.0) return x2;

if (fl < 0.0) { // Orient the search so that f(xl) < 0.
    xl=x1;
    xh=x2;
} else {
    xh=x1;
    xl=x2;
}
rts=0.5*(x1+x2);    //Initialize the guess for root,
dxold=fabs(x2-x1);  //the "stepsize before last,"
dx=dxold;           //and the last step

evalPoly(rts,coeffs,degree,&f);
evalPoly(rts,dcoeffs,degree-1,&dx);

for (j=1;j<=MAXIT;j++) { //Loop over allowed iterations

    if ((((rts-xh)*df-f)*((rts-xl)*df-f) > 0.0) //Bisect if Newton out of range,
            || (fabs(2.0*f) > fabs(dxold*df))) { //or not decreasing fast enough.
        dxold=dx;
        dx=0.5*(xh-xl);
        rts=xl+dx;
        if (xl == rts) 
            return rts; //Change in root is negligible.
    } else {// Newton step acceptable. Take it.
        dxold=dx;
        dx=f/df;
        temp=rts;
        rts -= dx;
        if (temp == rts)
            return rts;
    }
    if (fabs(dx) < xacc) 
        return rts;// Convergence criterion
    evalPoly(rts,coeffs,degree,&f);
    evalPoly(rts,dcoeffs,degree-1,&dx);
    //The one new function evaluation per iteration.
    if (f < 0.0) //Maintain the bracket on the root.
        xl=rts;
    else
        xh=rts;

}
//As the Accuracy asked to the algorithm is really high (but usually easily reached)
//the results precission is checked again, but with a less exigent result
dx=f/df;
if(fabs(dx)<xacc2)
    return rts;
nrerror("Maximum number of iterations exceeded in rtsafe");
return 0.0;// Never get here.
}

L'algorithme est appelé avec les variables suivantes :

x1=0.019
x2=1.05
xacc=1e-10
xacc2=0.1
degree=4
MAXIT=1000
coeffs[0]=1.75823776590795;
coeffs[1]=0.139389107255627;
coeffs[2]=-3.68254086033178;
coeffs[3]=0.557257315256597;
coeffs[4]=1.0;

Le problème est que l'algorithme dépasse le nombre maximum d'itérations et qu'il y a un arrondi à la racine d'approximativement 0.15 .

Ma question directe et abrégée est donc la suivante : pourquoi ce polynôme n'atteint-il pas une erreur précise alors que de nombreux (au moins 1000) autres polynômes très similaires y parviennent (avec une précision de 1e-10 et peu d'itérations !).

Je sais que la question est difficile et qu'il n'y a peut-être pas de réponse vraiment directe, mais je suis bloqué depuis quelques jours et je ne sais pas comment la résoudre. Merci beaucoup d'avoir pris le temps de lire ma question.

Answer 1

Sans avoir exécuté votre code, je pense que vous comparez des valeurs en virgule flottante pour l'égalité afin de déterminer si votre solution a convergé.

   if (xl == rts) 
        return rts; //Change in root is negligible.

Peut-être devriez-vous le calculer comme un ratio :

   diff = fabs(xl - rts);
   if (diff/xl <= 1.0e-8)  // pick your own accuracy value here
      return rts;

Answer 2

Je ne sais pas exactement pourquoi, mais la vérification de la vitesse de décroissance de la fonction ne semble pas fonctionner dans ce cas.

Cela fonctionne si je procède de la manière suivante :

  double old_f = f;

  .
  .
  .

    if ((((rts-xh)*df-f)*((rts-xl)*df-f) > 0.0) //Bisect if Newton out of range,
        || (fabs(2.0*f) > old_f)) { //or not decreasing fast enough.
  .
  .
  .

    if (fabs(dx) < xacc)
      return rts;// Convergence criterion
    old_f = f;

MISE À JOUR

Il semble qu'il y ait un problème dans votre code :

evalPoly(rts,dcoeffs,degree-1,&dx);

devrait être

evalPoly(rts,dcoeffs,degree-1,&df);