2

プロジェクトの一環として、C++ で閉じた形式で 4 次多項式を解く必要があります。

A*x 4 + B*x 3 + C*x 2 + D*x + E = 0

この目的へのリンクがいくつか見つかりました。それらの 1 つがここにあります。しかし、すべての根を計算しますが、実際の根だけが必要です。アルゴリズムは、主にフェラーリの方法を使用して次数を減らします。

bool solveQuartic(double a, double b, double c, double d, double e, double &root)
{
// I switched to this method, and it seems to be more numerically stable.
// http://www.gamedev.n...topic_id=451048 

// When a or (a and b) are magnitudes of order smaller than C,D,E
// just ignore them entirely. This seems to happen because of numerical
// inaccuracies of the line-circle algorithm. I wanted a robust solver,
// so I put the fix here instead of there.
if(a == 0.0 || abs(a/b) < 1.0e-5 || abs(a/c) < 1.0e-5 || abs(a/d) < 1.0e-5)
    return solveCubic(b, c, d, e, root);

double B = b/a, C = c/a, D = d/a, E = e/a;
double BB = B*B;
double I = -3.0*BB*0.125 + C;
double J = BB*B*0.125 - B*C*0.5 + D;
double K = -3*BB*BB/256.0 + C*BB/16.0 - B*D*0.25 + E;

double z;
bool foundRoot2 = false, foundRoot3 = false, foundRoot4 = false, foundRoot5 = false;
if(solveCubic(1.0, I+I, I*I - 4*K, -(J*J), z))
{
    double value = z*z*z + z*z*(I+I) + z*(I*I - 4*K) - J*J;

    double p = sqrt(z);
    double r = -p;
    double q = (I + z - J/p)*0.5;
    double s = (I + z + J/p)*0.5;

    bool foundRoot = false, foundARoot;
    double aRoot;
    foundRoot = solveQuadratic(1.0, p, q, root);
    root -= B/4.0;

    foundARoot = solveQuadratic(1.0, r, s, aRoot);
    aRoot -= B/4.0;
    if((foundRoot && foundARoot && ((aRoot < root && aRoot >= 0.0) 
        || root < 0.0)) || (!foundRoot && foundARoot)) 
    {
        root = aRoot;
        foundRoot = true;
    }

    foundARoot = solveQuadraticOther(1.0, p, q, aRoot);
    aRoot -= B/4.0;
    if((foundRoot && foundARoot && ((aRoot < root && aRoot >= 0.0) 
        || root < 0.0)) || (!foundRoot && foundARoot)) 
    {
        root = aRoot;
        foundRoot = true;
    }

    foundARoot = solveQuadraticOther(1.0, r, s, aRoot);
    aRoot -= B/4.0;
    if((foundRoot && foundARoot && ((aRoot < root && aRoot >= 0.0) 
        || root < 0.0)) || (!foundRoot && foundARoot)) 
    {
        root = aRoot;
        foundRoot = true;
    }
    return foundRoot;
}
return false;
}

これは、実数と虚数の両方の解を与える solveCubic() を使用します。

bool solveCubic(double &a, double &b, double &c, double &d, double &root)
{
if(a == 0.0 || abs(a/b) < 1.0e-6)
    return solveQuadratic(b, c, d, root);

double B = b/a, C = c/a, D = d/a;

double Q = (B*B - C*3.0)/9.0, QQQ = Q*Q*Q;
double R = (2.0*B*B*B - 9.0*B*C + 27.0*D)/54.0, RR = R*R;

// 3 real roots
if(RR<QQQ)
{
    /* This sqrt and division is safe, since RR >= 0, so QQQ > RR,    */
    /* so QQQ > 0.  The acos is also safe, since RR/QQQ < 1, and      */
    /* thus R/sqrt(QQQ) < 1.                                     */
    double theta = acos(R/sqrt(QQQ));
    /* This sqrt is safe, since QQQ >= 0, and thus Q >= 0             */
    double r1, r2, r3;
    r1 = r2 = r3 = -2.0*sqrt(Q);
    r1 *= cos(theta/3.0);
    r2 *= cos((theta+2*PI)/3.0);
    r3 *= cos((theta-2*PI)/3.0);

    r1 -= B/3.0;
    r2 -= B/3.0;
    r3 -= B/3.0; 

    root = 1000000.0;

    if(r1 >= 0.0) root = r1;
    if(r2 >= 0.0 && r2 < root) root = r2;
    if(r3 >= 0.0 && r3 < root) root = r3;

    return true;
}
// 1 real root
else
{
    double A2 = -pow(fabs®+sqrt(RR-QQQ),1.0/3.0);
    if (A2!=0.0) {
        if (R<0.0) A2 = -A2; 
        root = A2 + Q/A2; 
    }
    root -= B/3.0;
    return true;
}
}

コードを説明するリンクをいくつか示します。solveCubicとsolveQuartic

実根の 4 次多項式を解くためにコードを変更できる人はいますか?

できるだけ効率よく導入したい。ところで、誰かがこの目的のために LAPACK のような便利なライブラリを導入してくれたら幸いです (4 次多項式の根を直接計算することはできないようです)。

4

1 に答える 1

1

おそらく、この方程式を実根の閉形式で解く最も効率的な方法は、すべての根の閉形式で解いてから、虚数の根を破棄することです。

try/catch ペアを使用して虚数が発生しているかどうかを判断できると思うかもしれませんが、実根の計算で生成する中間値の一部が虚数である可能性があるため、これはあまり良い戦略ではありません。

したがって、C++ 複素数ライブラリ (こちらまたはこちらを参照) を使用して計算を行うことができます。

その後、数値の虚数部がゼロでないかどうかを確認し、ゼロである場合は破棄します。ただし、浮動小数点演算は不正確であるため、「ゼロ」にはゼロに非常に近い数値の範囲が含まれることを覚えておいてください。

于 2012-10-20T19:15:56.037 に答える