3

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

 7 months ago
source link: https://blog.darkthread.net/blog/guassian-elimination/
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.
neoserver,ios ssh client

Coding4Fun - 高斯消去法解聯立方程式-黑暗執行緒

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

Fig1_638435121923942262.png

在 .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();
    }
}

執行結果如下，成功。

Fig2_638435121926223442.png

註：範例程式已更新至 Github


report

Recommend

About Joyk


Aggregate valuable and interesting links.
Joyk means Joy of geeK