Coding4Fun - 高斯消去法解聯立方程式
0 | 1,079 |
在 .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
Comments
Be the first to post a comment