Coding4Fun - 高斯消去法解聯立方程式

用電腦解聯立方程式，一般會使用高斯消去法用矩陣求解：

在 .NET 做矩陣運算，Math.NET Numerics 程式庫是首選。

Math.Net Numerics 是個開源數學程式庫，囊括矩陣、線性代數求解、機率、線性回歸、積分、傅立葉轉換... 等各種你想得到的數學運算，且被學術論文和期刊廣泛引用，其專業性毋庸置疑。總之，遇到 .NET 算數學需求，先找 Math.Net Numerics 就對了，別急著自己造輪子。

練習題為假設有聯立方程式組如下，求解 x, y, z 為何？

7x - 2y + 3z = 12
2x - y + 4z = 12
-x + 3y - z = 2

先用 MathNet.Numerics.LinearAlgebra 函式庫求解。在專案 dotnet add package MathNet.Numerics 加入參照，只需幾行程式建個 Matrix 物件呼叫 .Solve() 便能完成：參考

using MathNet.Numerics.LinearAlgebra;
public class MathNetDemo
{
    public static void Solve()
    {
        var matrix = Matrix<double>.Build.DenseOfArray(new double[,]
        {
            {7, -2, 3},
            {2, -1, 4},
            {-1, 3, -1}
        });
        var vector = Vector<double>.Build.Dense(new double[] {12, 12, 2});
        var result = matrix.Solve(vector);
        Console.WriteLine($"x = {result[0]}, y = {result[1]}, z = {result[2]}");
    }
}

小缺點是 MathNet 的 Matrix/Vector 的數字只支援 Double、Single、Complex 及 Complex32 四種型別(Matrices and vectors of type 'Frac' are not supported. Only Double, Single, Complex or Complex32 are supported at this point.)，使用 double 浮點數計算會包含浮點數誤差(Floating-Point Error)，得到 x = 1.0000000000000002, y = 2, z = 2.9999999999999996，有時遇到分數無法整除只能以循環小數表示，雖然實務應用不致造成問題，但在做數學習題時，這答案就是不完美。

前幾天寫的分數型別剛好可以避開上述問題，而高斯消去法邏輯不難，我打算自己實作玩看看，順便把解題過程印出來幫助了解原理，以上就是本次的計劃。

用矩陣運算求解聯立方程式的概念是將常數值轉成矩陣形式(我用 Frac[,] 二維陣列模擬)，先排序將每一行最大值移到對角線，接著將對角線係數調成 1，對角線以下的係數調成 0 形成列階梯形矩陣，最後再將非對角線係數也調成 0，即可得解。

實作程式邏輯如下：(註：我有為 Frac 增加 Abs() 方法 public Frac Abs() => this > 0 ? this : -this; 計算絕對值)

public class FracMatrixDemo 
{
    public static void GaussianEliminationDemo() 
    {
/*
聯立方程組
7x - 2y + 3z = 12
2x - y + 4z = 12
-x + 3y - z = 2
*/        
        Frac[,] matrix = new Frac[3, 4] 
        {
            {7, -2, 3, 12},
            {2, -1, 4, 12},
            {-1, 3, -1, 2}
        };
        PrintMatrix(matrix, "原矩陣");
        // 排序將絕對值最大的列放在最上面
        SortRows(matrix);
        PrintMatrix(matrix, "排序後");
        // 將矩陣化為上三角矩陣
        Console.WriteLine("列階梯形矩陣");
        ToRowEchelonForm(matrix);
        // 將對角線以外係數調成 0
        Console.WriteLine("簡化列階梯形矩陣");
        ToDiagonal(matrix);
    }

    static void ToRowEchelonForm(Frac[,] matrix)
    {
        var rowCount = matrix.GetLength(0);
        var colCount = matrix.GetLength(1);
        for (int i = 0; i < rowCount; i++)
        {
            Console.WriteLine($"** 第 {i + 1} 列 **");
            // 將對角線係數調成 1
            var diag = matrix[i, i];
            for (int j = i; j < colCount; j++)
            {
                matrix[i, j] /= diag;
            }
            PrintMatrix(matrix, $"({i}, {i}) 改為 1");
            // 將對角線以下係數調成 0
            for (int k = i + 1; k < rowCount; k++)
            {
                var factor = matrix[k, i];
                for (int l = i; l < colCount; l++)
                {
                    matrix[k, l] -= matrix[i, l] * factor;
                    if (k == l && matrix[k, l] == 0) 
                    {
                        throw new ApplicationException("無限多組解");
                    }                    
                }
            }
            if (i < rowCount - 1)
                PrintMatrix(matrix, $"({i}, {i}) 以下改為 0");            
        }
    }

    static void ToDiagonal(Frac[,] matrix)
    {
        var rowCount = matrix.GetLength(0);
        var colCount = matrix.GetLength(1);
        for (int i = rowCount - 1; i > 0; i--)
        {
            Console.WriteLine($"** 第 {i + 1} 列 **");
            for (int j = i - 1; j >= 0; j--)
            {
                var factor = matrix[j, i];
                for (int k = i; k < colCount; k++)
                {
                    matrix[j, k] -= matrix[i, k] * factor;
                }
            }
            PrintMatrix(matrix, $"({i}, {i}) 以上改為 0");
        }
    }

    static void SortRows(Frac[,] matrix)
    {
        // 排序將對角線元素絕對值較大者放在最上面
        var rowCount = matrix.GetLength(0);
        var colCount = matrix.GetLength(1);
        for (int i = 0; i < rowCount; i++)
        {
            int maxRow = i;
            for (int j = i + 1; j < rowCount; j++)
            {
                if (matrix[j, i].Abs() > matrix[maxRow, i].Abs())
                {
                    maxRow = j;
                }
            }
            // 與最大列交換
            if (maxRow != i)
            {
                for (int k = 0; k < colCount; k++)
                {
                    var temp = matrix[i, k];
                    matrix[i, k] = matrix[maxRow, k];
                    matrix[maxRow, k] = temp;
                }
            }
        }

    }

    public static void PrintMatrix(Frac[,] matrix, string title)
    {
        var len = Math.Max(5, matrix.Cast<Frac>().Max(o => o.ToString().Length) + 2);
        Console.ForegroundColor = ConsoleColor.Cyan;
        Console.WriteLine(title);
        Console.ForegroundColor = ConsoleColor.Yellow;
        for (int i = 0; i < matrix.GetLength(0); i++)
        {
            for (int j = 0; j < matrix.GetLength(1); j++)
            {
                Console.Write(matrix[i, j].ToString().PadLeft(len));
            }
            Console.WriteLine();
        }
        Console.WriteLine();
        Console.ResetColor();
    }
}

執行結果如下，成功。

註：範例程式已更新至 Github

This article introduces how to solve a system of simultaneous equations using the MathNet library and implement the Gaussian elimination method with a custom fraction type to obtain accurate solutions