三次方程式の解

今日はアルゴリズムのお話を。

FEMのポスト処理ソフトを作っていると三次方程式の解き方が必要な時が多々でてきます。

三次方程式の出現元は2階のテンソルの主値を求める処理。

例えば
- 応力テンソルから最大主応力を求めたり
- 剛体の慣性テンソル（らしきもの）からバウンディングボックスの各軸方向を求めたり。

バウンディングボックスについてはバウンディングボックス - malibu-bulldogの日記を。

さてさて2階のテンソルが出現元なので，三次方程式を解くことはようは行列の固有値を求める事と結果が同じになります。

行列の固有値を求めるのは

2008-06-15 - malibu-bulldogの日記

に書いたヤコビ法があります。

でも，ヤコビ法は反復法なんで直接法があればそっちの方が計算はやいかも…。

つまり直接テンソルの主値に関する三次方程式が解けるならそっちの方がいいかもしれないってことです。

さてさて，応力テンソルも慣性テンソルも実対称テンソルなので複素数という鬼門は通らずに進めます。つまりカルダノの方法をチョイスできます。

三次方程式 - Wikipedia

カルダノの方法を使った三次方程式の解法をCのプログラムに起こすと，

対称のテンソルを3x3の配列として引数にとる関数は次のようになります。

int solve( double m[3][3], double x[3] )

{

    double a = -1.0 * m[0][0] -  m[1][1] - m[2][2];

    double b = m[0][0]*m[2][2] + m[1][1]*m[2][2] + m[0][0]*m[1][1]

             - m[1][2]*m[2][1] - m[0][1]*m[1][0] - m[0][2]*m[2][0];

    double c = m[0][1]*m[1][0]*m[2][2] - m[0][1]*m[1][2]*m[2][0]

             + m[0][2]*m[1][1]*m[2][0] - m[0][2]*m[1][0]*m[2][1]

             + m[0][0]*m[1][2]*m[2][1] - m[0][0]*m[1][1]*m[2][2];

    double p = -1.0 * a * a / 9.0 + b / 3.0;

    double q = a * a * a /27.0 - a * b / 6.0 + c / 2.0;



    double D = p * p * p + q * q;


    if( D > 0.0 ){
        double r = cbrt( -1.0 * q + sqrt( D ) );

        double s = cbrt( -1.0 * q - sqrt( D ) );


        x[0] = r + s -1.0 * a / 3.0;

        return 1;

    }



    if( D == 0.0 ){
        double r = cbrt( -1.0 * q + sqrt( D ) );

        x[0] = 2.0 * r - 1.0 * a / 3.0;

        x[1] = x[2] = -1.0 * r - 1.0 * a / 3.0;

        return 2;

    }

    if( D  < 0 ){
        double theta =acos( -1.0 * q / sqrt( -1.0 * p * p * p ) ) / 3.0;

        x[0] = 2.0 * sqrt( -1.0 * p ) * cos( theta ) - 1.0 * a / 3.0;

        x[1] = 2.0 * sqrt( -1.0 * p ) * cos( theta + 2.0 * M_PI / 3.0) - 1.0 * a / 3.0;

        x[2] = 2.0 * sqrt( -1.0 * p ) * cos( theta - 2.0 * M_PI / 3.0) - 1.0 * a / 3.0;

    }

    return 3;

}

上記で求まった，x[0],x[1],x[2]がテンソルの主値になります。

主値に対応する主方向（バウンディングボックスを求めるときはこっちの方が重要かな）はハミルトンケーリーの公式（だったかな）より行列の演算だけで簡単にもとまります（解が三つの実数の時に）。

そこら辺のテクニックは

http://www004.upp.so-net.ne.jp/s_honma/urawaza/eigenvector.htm

にまとまってます。

解が二つとか一つだった場合は…特別な処理が必要になります。