高斯消元法的主要思想是通过对系数矩阵进行行变换,将方程组的系数矩阵由对称矩阵变为三角矩阵,从而达到消元的目的,最后通过回代逐个获得方程组的解。
高斯消元法的实现简单,主要由两个步骤组成:第一个步骤就是通过选择主元,逐行消元,最终形成方程组系数矩阵的三角矩阵形式;第二步骤就是逐步回代的过程,最终矩阵的对角线上的元素就是方程组的解。
using System;
using System.Collections.Generic;
namespace Test4
{
class Program
{
static void Main(string[] args)
{
//测试数据
/*
* 2x + 3y + z = 6
* x - y + 2z = -1
* x + 2y - z = 5
*/
float[] matrix = new float[] { 2, 3, 1,
1, -1, 2,
1, 2, -1};
float[] value = new float[] { 6, -1, 5 };
float[] result = GuassEquation.Resolve(matrix, 3, value);
if (null != result)
{
Console.WriteLine("解: x={0}, y={1}, z={2}", result[0], result[1], result[2]);
}
Console.Read();
}
}
/// <summary>
/// 高斯消元法求解方程组
/// </summary>
public class GuassEquation
{
//矩阵维数
private static int m_DIM;
//一维数组形式存放的系数矩阵(必须是方阵)
private static float[] m_Matrix;
//方程组的值
private static float[] m_Value;
public static float[] Resolve(float[] matrix, int dim, float[] value)
{
if (matrix.Length / dim != dim || value.Length != dim)
{
throw new ArgumentException("Error, not square matrix.");
}
m_Matrix = matrix;
m_DIM = dim;
m_Value = value;
/*消元,得到上三角阵*/
for(int i = 0; i < m_DIM - 1; i++)
{
/*按列选主元*/
int pivotRow = SelectPivotalElement(i);
if (pivotRow != i)
{
SwapRow(i, pivotRow);
}
//主元为0时,不存在唯一解
if (IsPrecisionZero(i))
{
Console.WriteLine("不存在唯一解");
return null;
}
/*对系数归一化处理,使每行的第一个系数是1.0*/
SimplePivotalRow(i, i);
/*逐行进行消元*/
for (int j = i + 1; j < m_DIM; j++)
{
RowElimination(i, j);
}
}
/*回代求解*/
//对最后一个变量求解
m_Matrix[(m_DIM - 1) * m_DIM + m_DIM - 1] = m_Value[m_DIM - 1] / m_Matrix[(m_DIM - 1) * m_DIM + m_DIM - 1];
//对剩余变量求解
for(int i = m_DIM - 2; i >= 0; i--)
{
double totalCof = 0.0;
//将第i行对角线项之后的项相加
for (int j = i + 1; j < m_DIM; j++)
{
totalCof += m_Matrix[i * m_DIM + j] * m_Matrix[j * m_DIM + j];
}
//求出对角线上变量的解
m_Matrix[i * m_DIM + i] = (float)((m_Value[i] - totalCof) / m_Matrix[i * m_DIM + i]);
}
/*将对角线元素的解逐个存入解向量*/
float[] xValue = new float[m_DIM];
for (int i = 0; i < m_DIM; i++)
{
xValue[i] = m_Matrix[i * m_DIM + i];
}
return xValue;
}
// 判断主元是否为0
private static bool IsPrecisionZero(int i)
{
return m_Matrix[i * m_DIM + i] == 0;
}
/// <summary>
/// 选择主元
/// 在消元的过程中,如果某一行的对角线元素的值太小,
/// 在计算过程中就会出现很大的数除以很小的数的情况,
/// 有除法溢出的可能,因此在消元的过程中,通常都会增加
/// 一个主元选择的步骤。
/// 再通过行交换操作,将当前列绝对值最大的行交换到当前位置,
/// 避免了除法溢出的问题,增加了算法的稳定性。
/// </summary>
/// <param name="i">第i行</param>
/// <returns></returns>
private static int SelectPivotalElement(int i)
{
int col = i;
int pivotalRow = i;
float max_val = 0;
for (; i < m_DIM; i++)
{
float val = Math.Abs(m_Matrix[i * m_DIM + col]);
if (val > max_val)
{
max_val = val;
pivotalRow = i;
}
}
return pivotalRow;
}
// 交换行
private static void SwapRow(int i, int pivotRow)
{
for (int j = 0; j < m_DIM; j++)
{
float tmp = m_Matrix[i * m_DIM + j];
m_Matrix[i * m_DIM + j] = m_Matrix[pivotRow * m_DIM + j];
m_Matrix[pivotRow * m_DIM + j] = tmp;
}
float tmp_val = m_Value[i];
m_Value[i] = m_Value[pivotRow];
m_Value[pivotRow] = tmp_val;
}
// 对系数归一化处理
private static void SimplePivotalRow(int row, int col)
{
for (; row < m_DIM; row++)
{
float a = m_Matrix[row * m_DIM + col];
for (int j = col; j < m_DIM; j++)
{
m_Matrix[row * m_DIM + j] /= a;
}
m_Value[row] /= a;
}
}
// 消元
private static void RowElimination(int i, int j)
{
int i_start = i * m_DIM;
int j_start = j * m_DIM;
//第j行与第i行各项相减
for (int m = 0; m < m_DIM; m++)
{
m_Matrix[j_start + m] -= m_Matrix[i_start + m];
}
m_Value[j] -= m_Value[i];
}
public static void Print()
{
for (int i = 0; i < m_DIM; i++)
{
for (int j = 0; j < m_DIM; j++)
{
Console.Write(m_Matrix[i * m_DIM + j] + ",");
}
Console.WriteLine();
}
}
}
}
测试