多項式の実根を計算するアルゴリズム

このアルゴリズムの根底にある考え方は非常に単純で、2つの文で述べることができます。 多項式の実際の根は、常に多項式の単調変化の部位にあります。 多項式の導関数の根の間。 ただし、多項式の導関数も多項式ですが、その程度は小さく、その実根を見つけたので、半分の方法で導関数の根の間で元の多項式の根を探す必要があります。



そして今、順番に。



代数多項式の根を見つける問題は、少なくとも中世以来、長い間知られています。 学校では、二次方程式を解くことを教えます。 ウィキペディアでは、3次方程式を解くためのカルダノの式と、ラジカルの4次方程式を解くためのフェラーリ法の説明を見つけることができます。 また、任意の次数の代数方程式を解くためのロバチェフスキー法も説明されています。 ロバチェフスキー法の本質は次のように簡潔に要約されます。



元の多項式と同じ係数の根を持ち、反対の符号を持つ多項式2を構築することは、初期多項式を使用して難しくありません。 元の多項式と多項式2を乗算すると、元の多項式の根の平方に等しい根をもつ多項式が得られます。



この変換(二乗)は、数回繰り返すのに便利です。 その結果、元の多項式の根が互いに等しくなかった場合、それらの複数の2乗値の大きさは大きく離れ、それらの近似値は対応する2次多項式の係数によって非常に単純に表現されます。



特に、多項式の引数の最大次数の係数が1に等しい場合、次に高い係数は(反対の符号で)方程式の根の合計であり、これらの根の値は強く間隔が空いているので、この合計は絶対値の最大根に等しいと見なすことができます。



具体的には、ルート1、2、3、4の次数4の多項式の場合、4番目の四分割後のロバチェフスキー法により、小数点第2位まで正しいルートの値が得られることをお知らせします。 同時に、多項式の係数を表すにはlong double形式で十分です。



間違いなく、この方法は知性に恵まれた研究者の手にある貴重なツールです。 しかし、現代のコンピューターテクノロジー向けのプログラミングは、根の場所のあらゆる種類の特殊なケースで結果の信頼性を厳密に保証する必要がある場合、深刻な問題を引き起こします。



次に、別の方法について説明します。 公開報道では、それへの参照は[1]で始まります。 そのような方法の適用に関する独立した出版物は知りません。 このアルゴリズムは、初期多項式の単調変化の間隔の逐次研究に還元されます。 多項式値がこの単調間隔の境界で異なる符号を持っている場合、セグメントを半分に分割する手順は、次のルートの正確な値の計算を開始します。 単調間隔の境界は、多項式の導関数が消滅するポイントです。 これらは導出された多項式の根です。 導出された多項式は、元の多項式よりも次数が1つ少ないため、多項式の導関数の係数を計算するプロセスは、セグメントを半分に分割する手順を必要とせずに、その根が直接配置されている1次の多項式まで継続する必要があります。 このステップの結果として、派生した2次多項式の単調変動の2つの区間を取得します。



これで、2次の微分多項式の2つの実根(存在する場合)を見つけ、多項式の微分から元の多項式の根まで階段を上ることができます。 元の多項式と導出された多項式の単調性区間のプラス、マイナスの無限境界がどのように技術的に実現されるかは、まだ説明されていません。



引数の次数が最も高い係数が1に等しくなるように、多項式を正規化します。 Mを他の係数の中で最大のモジュロ値とします。 次に、M + 1より大きい引数のすべての値について、多項式の値が1より大きくなります。



これを証明するために、ホーナースキームによる多項式p(x)= x ^ n + k [n-1] * x ^(n-1)+ ... + k [1] * x + k [0]の計算を検討します。



最初のステップで、 p [1] = k [n-1] + xが計算され、 p [1]> 1であることは明らかです。

2番目のステップでは、 p [2] = k [n-2] + x * p [1]が計算され、再びp [2]> 1であることが明らかです。

同様のことが次のステップで行われます。

最後のステップで、 p(x)= k [0] + x * p [n-1]が計算され、最終的にp(x)> 1が得られます。



したがって、引数の無限値を使用して多項式の符号を決定する必要がある場合、M + 1に等しい引数を取る必要があります。



対応するプログラムの添付テキストは、ここで説明するアルゴリズムの個々の技術的詳細の退屈なプレゼンテーションを完全に置き換えます。



最後に、セグメントを半分に分割するアルゴリズムの実装のあまり明らかではない機能についてコメントします。



セグメントの現在の端ngとvgの中間にあるテストポイントptは、演算子pt = 0.5 *(ng + vg)によって計算されます。 (pt <= ng || pt> = vg)breakの場合 、半分の除算のサイクルはオペレーターによって中断されます。



マシンで実数を表現する際の精度は有限であるため、遅かれ早かれ、新しい数値ではなく半分に分割する操作が初期境界の1つの値を与える状態が発生します。 半分のサイクルを停止するのはこの状態です。 この状態は、結果の最大達成可能精度に対応します。



最近、このアルゴリズムを使用して、実数根を持たない多項式の複素根を計算する問題を解決することができました。 しかし、次の記事でHabréでそれについて話す予定です。



以下に、アプリケーションとして、説明したアルゴリズムを実装するpolynomRealRoots.cppファイルの全文を示します。



ファイルpolynomRealRoots.cpp
//************************************************************************* double polinom(int n,double x,double *k) //  x^n+k[n-1]*x^(n-1)+k[n-2]*x^(n-2)+...+k[1]*x+k[0] //  ,   { double s=1; for(int i=n-1;i>=0;i--) s=s*x+k[i]; return s; }//polinom //      ,    double dihot(int degree,double edgeNegativ,double edgePositiv,double *kf) { for(;;) {//    double x=0.5*(edgeNegativ+edgePositiv); if(x==edgeNegativ||x==edgePositiv)return x; if(polinom(degree,x,kf)<0)edgeNegativ=x; else edgePositiv=x; }//    }//dihot //        void stepUp( int level, //    double **A, //  double **B, //   int *currentRootsCount //    ) { // ,      double major=0; for(int i=0;i<level;i++) {// major double s=fabs(A[level][i]); if(s>major)major=s; }// major major+=1.0; currentRootsCount[level]=0; //  //     //A[level][0]+A[level][1]*x+...+A[level][level-1]*x^(level-1)+x^level=0 //    for(int i=0;i<=currentRootsCount[level-1];i++) {//   //signLeft signRight -   A-        int signLeft,signRight; //       double edgeLeft,edgeRight; //  ,        double edgeNegativ,edgePositiv; //    if(i==0)edgeLeft=-major; else edgeLeft=B[level-1][i-1]; //  A-    double rb=polinom(level,edgeLeft,A[level]); if(rb==0) {//     B[level][currentRootsCount[level]]=edgeLeft; currentRootsCount[level]++; continue; }//     //   A-    if(rb>0)signLeft=1;else signLeft=-1; //    if(i==currentRootsCount[level-1])edgeRight=major; else edgeRight=B[level-1][i]; //  A-    rb=polinom(level,edgeRight,A[level]); if(rb==0) {//     B[level][currentRootsCount[level]]=edgeRight; currentRootsCount[level]++; continue; }//     //   A-    if(rb>0)signRight=1;else signRight=-1; //       , //   if(signLeft==signRight)continue; //          if(signLeft<0){edgeNegativ=edgeLeft;edgePositiv=edgeRight;} else {edgeNegativ=edgeRight;edgePositiv=edgeLeft;} //          B[level][currentRootsCount[level]]=dihot(level,edgeNegativ,edgePositiv,A[level]); currentRootsCount[level]++; }//   return; }//stepUp //        void polynomRealRoots( int n, //   double *kf, //    double *rootsArray, //    int &rootsCount //   ) { //   A  B,    //A    B    // A-  , //        //A[n] -       //B[n] -      //A[n-1] -       //B[n-1] -      //  //A[n-2]  B[n-2] -        // A[1] -      //    //    B[1] -    , //    //     double **A=new double *[n+1]; double **B=new double *[n+1]; //      A- int *currentRootsCount=new int[n+1]; for(int i=1;i<=n;i++) { A[i]=new double[i]; B[i]=new double[i]; } //   for(int i=0;i<n;i++)A[n][i]=kf[i]/kf[n]; //  A- for(int i1=n,i=n-1;i>0;i1=i,i--) { for(int j1=i,j=i-1;j>=0;j1=j,j--) { A[i][j]=A[i1][j1]*j1/i1; } } //      currentRootsCount[1]=1; B[1][0]=-A[1][0]; //     for(int i=2;i<=n;i++)stepUp(i,A,B,currentRootsCount); //  rootsCount=currentRootsCount[n]; for(int i=0;i<rootsCount;i++)rootsArray[i]=B[n][i]; //  for(int i=1;i<=n;i++) { delete[]A[i]; delete[]B[i]; } delete[]A; delete[]B; delete[]currentRootsCount; return; }//polynomRealRoots //*************************************************************************
      
      









ヘッダーファイルpolynomRealRoots.hのテキストも受け入れます。これにより、上記のプログラムモジュールへのリンクを簡単に整理できます。



 //************************************************************************* //     void polynomRealRoots(int n,double *kf,double *rootsArray,int &rootsCount); //**
      
      



***************************************************** **********************



文学



1.コスティンI.K. 軌道のパラメータを決定する際の縦方向の誤差を考慮に入れて、円形の地上物体上の宇宙船の通過間隔を計算するための一連のアルゴリズム



ラジオエレクトロニクスの質問、ser。 RLT、2004年、発行。 1



このリンクは、引用された語句「計算アルゴリズムのファミリー」を検索することでYandexで見つけることができますが、電子形式のこの記事のテキストは利用できないようです。 したがって、この記事の2つの文から引用します。

多項式の実際の根は、常に多項式の単調変化の部位にあります。 多項式の導関数の根の間。



ただし、多項式の導関数も多項式ですが、その程度は小さく、その実根を見つけたので、半分の方法で導関数の根の間で元の多項式の根を探す必要があります。



All Articles