クラメルの公式

たしか動いたはずのもの。
C#で書かれたもの。

以下の書物に記述がある。
ゲームプログラミングのためのリアルタイム衝突判定p29

理工系の数理 線形代数
永井 敏隆 (著), 永井 敦 (著)p103

連立方程式が以下の形式で与えられている時

$$
\begin{equation}
\left\{
\begin{aligned}
& a_{11}x_1+a_{12}x_2+…a_{1n}x_n=b_1 \\
& a_{21}x_1+a_{22}x_2+…a_{2n}x_n=b_2 \\
& \vdots\\
& a_{n1}x_1+a_{n2}x_2+…a_{nn}x_n=b_n \\
\end{aligned}
\right.
\end{equation}
$$

ここで係数aを行列Aにまとめて,xとbもまとめる。

$$
A=\begin{bmatrix}
a_{11} & a_{12} & \cdots & a_{1n} \\
a_{21} & a_{22} & \cdots & a_{2n} \\
\vdots & \vdots & \ddots & \vdots \\
a_{n1} & a_{n2} & \cdots & a_{nn} \\
\end{bmatrix},
\quad
x=\begin{bmatrix}
x_1 \\
x_2 \\
\vdots \\
x_n  \\
\end{bmatrix},
\quad
b=\begin{bmatrix}
b_1 \\
b_2 \\
\vdots \\
b_n  \\
\end{bmatrix},

$$

この時、求めるべき解の成分は

$$
x_i=\frac{det(A_i)}{det(A)}=\frac{1}{det(A)}det(A_i)
$$

ただし、$${A_i}$$はA行列のi列をbに置き換えたもの。
以下のvar Ai = PasteCol(A, x_index, b);部分がそれにあたる。
CramersRule関数は解であるxベクトルの成分を一つ求める。

       //クラメルの公式要素１つ分
       public double CramersRule(double[,] A, int x_index, double[] b)
       {
           var Ai = PasteCol(A, x_index, b);
           return Det(Ai) / Det(A);
       }

        public static double[,] PasteCol(double[,] matrix, int col_index, double[] new_col)
        {
            int row_count = matrix.GetLength(0);//rowの個数
            int col_count = matrix.GetLength(1);//colの個数
            double[,] ret = new double[row_count, col_count];
            for (int row_index = 0; row_index < row_count; row_index++)
            {
                //入力配列が短い場合は0で埋める
                //必然的に、長い場合は切り落とされる
                if (new_col.Length - 1 < row_index)
                {
                    ret[row_index, col_index] = 0;
                }
                ret[row_index, col_index] = new_col[row_index];
            }
            return ret;
        }

解のxベクトルのすべての成分を求めると以下の形式となる。
ただし以下の場合、xは適切な長さの配列で初期化してなければならず、そもそもインデックスしか使わないので引数はCramersRules(double[,] A,  double[] b)の形に直すことができる。

この関数はクラメルの公式を用いて連立方程式を解いたのと同じである。

       public double[] CramersRules(double[,] A, double[] x, double[] b)
       {
           double[] ret = new double[x.Length];
           for (int x_index = 0; x_index < x.Length; x_index++)
           {
               ret[x_index] = CramersRule(A, x_index, b);
           }
           return ret;
       }

↓直した

       public double[] CramersRules(double[,] A, double[] b)
       {
           double[] ret = new double[b.Length];
           for (int x_index = 0; x_index < b.Length; x_index++)
           {
               ret[x_index] = CramersRule(A, x_index, b);
           }
           return ret;
       }

応用し、余因子展開を用いると逆行列を求めることもできる。

$$
A^{-1}=\frac{1}{det(A)}adj(A)
$$

ここでdetは行列式。

行列式

クラメルの公式を用いるには行列式を求める必要がある。

N×N次行列の行列式を力づくで求めるには余因子展開をする必要がある。
余因子展開するためには余因子を求める必要がある。余因子を求めるためには小行列を求めてその行列式である小行列式を求める必要がある。

このように、行列式をまともに相手するとコンピューターが死ぬので、プログラム的に行列式を求めたい時はガウスの消去法ないし、各種反復法を用いる。

行列式をそのままコードで求めるもの

       //行列式略記
       public static double Det(double[,] matrix)
       {
           return DeterminantN(matrix);
       }
       
       public static double DeterminantN(double[,] matrix)
       {
           int m = matrix.GetLength(0);
           int n = matrix.GetLength(1);

           if (m != n) { throw new Exception(); }
           else if (m == 1) { return matrix[0, 0]; }
           else if (m == 2) { return Determinant2(matrix); }
           else if (m == 3) { return Determinant3(matrix); }
           else { return CofactorExpansion(matrix); }
       }
       
       //行列式
       public static double Determinant2(double[,] m)
       {
           //二次の場合
           //|a b|
           //|c d|
           var a = m[0, 0];
           var b = m[0, 1];
           var c = m[1, 0];
           var d = m[1, 1];
           return a * d - b * c;
       }
       
       public static double Determinant2(
           double a, double b,
           double c, double d)
       {
           //二次の場合
           //|a b|
           //|c d|
           return a * d - b * c;
       }
       
       public static double Determinant3(double[,] m)
       {
           //三次の場合
           //|a b c|
           //|d e f|
           //|g h i|

           var a = m[0, 0];
           var b = m[0, 1];
           var c = m[0, 2];

           var d = m[1, 0];
           var e = m[1, 1];
           var f = m[1, 2];

           var g = m[2, 0];
           var h = m[2, 1];
           var i = m[2, 2];

           return a * Determinant2(e, f, h, i) - d * Determinant2(b, c, h, i) + g * Determinant2(b, c, e, f);
       }
       
       public static double Determinant3(
           double a, double b, double c,
           double d, double e, double f,
           double g, double h, double i)
       {
           //三次の場合
           //|a b c|
           //|d e f|
           //|g h i|
           return a * Determinant2(e, f, h, i) - d * Determinant2(b, c, h, i) + g * Determinant2(b, c, e, f);
       }

N×N次の行列の行列式を求めるには余因子展開（上のコードではCofactorExpansion）をして、再帰的に(N-1)×(N-1)次の行列式を求めることになる。２×２や３×３の行列式になったら普通に計算する。

余因子展開

余因子展開するためには小行列式を求める必要があるため小行列を取得する必要がある。

       //余因子展開
       public static double CofactorExpansion(double[,] matrix)
       {
           //ループはどこかの１列、あるいは１行のみ行えばよい
           //どの行、列を選択しても結果は同じとなる            
           int col_index = 0;
           int row_count = matrix.GetLength(0);
           //int col_count = matrix.GetLength(1);
           double ret = 0;
           for (int row_index = 0; row_index < row_count; row_index++)
           {
               ret += matrix[row_index, col_index] * CalcuCofactor(matrix, row_index, col_index);
           }
           return ret;
       }
       
       //余因子を取得する
       public static double CalcuCofactor(double[,] matrix, int input_row_index, int input_col_index)
       {
           return _sign(input_row_index + input_col_index) * Det(GetMinorMatrix(matrix, input_row_index, input_col_index));
       }

小行列式

小行列を取得する

       //小行列を取得する(略記）
       public static double[,] M(double[,] matrix, int i, int j)
       {
           return GetMinorMatrix(matrix, i, j);
       }
       
       //小行列を取得する
       public static double[,] GetMinorMatrix(double[,] matrix, int input_row_index, int input_col_index)
       {
           int row_count = matrix.GetLength(0);//row_count, 行数
           int col_count = matrix.GetLength(1);//col_count, 列数
           if (row_count != col_count) { return null; }
           int ret_row_index = 0;
           int ret_col_index = 0;
           double[,] ret = new double[row_count - 1, col_count - 1];
           for (int r = 0; r < row_count; r++)
           {
               if (r == input_row_index) { continue; }
               for (int c = 0; c < col_count; c++)
               {
                   if (c == input_col_index) { continue; }
                   ret[ret_row_index, ret_col_index] = matrix[r, c];
                   ret_col_index++;
               }
               ret_col_index = 0;
               ret_row_index++;
           }
           return ret;
       }

       //符号
       public static int _sign(int value)
       {
           if (value % 2 == 0)
           {
               //偶数
               return 1;
           }
           else
           {
               //奇数
               return -1;
           }
       }

余因子行列adj(A)を用いた逆行列

$$
A^{-1}=\frac{1}{det(A)}adj(A)
$$

余因子行列は行列Aの各成分の余因子を成分とする行列の転置をとったもの。

        //随伴行列、余因子行列
        //要素を余因子とした行列の転置
        public double[,] AdjugateMatrix(double[,] matrix)
        {
            //転置
            return Transpose(CofactorMatrix(matrix));
        }

        //余因子行列を用いた逆行列
        public double[,] InverseMatrixFromCofactor(double[,] matrix)
        {
            return Divide(AdjugateMatrix(matrix), Det(matrix));
        }

転置

        public static double[,] Transpose(double[,] matrix)
        {
            int row_count = matrix.GetLength(0);
            int col_count = matrix.GetLength(1);

            //反転
            double[,] ret = new double[col_count, row_count];

            for (int row_index = 0; row_index < row_count; row_index++)
            {
                for (int col_index = 0; col_index < col_count; col_index++)
                {
                    ret[col_index, row_index] = matrix[row_index, col_index];
                }
            }
            return ret;
        }

Divideは行列の要素ごとに割り算する演算

double[,] Divide(double[,] matrix, double val)
{
    double[,] ret = new double[matrix.GetLength(0), matrix.GetLength(1)];
    for(int i = 0; i<matrix.GetLength(0); i++)
    {
        for(int j = 0; j<matrix.GetLength(1); j++)
        {
            ret[i,j]=matrix[i,j]/val;
        }
    }
    return ret;
}