高斯消元法求解方程组(C#实现)

作者：追风剑情发布于：2017-11-5 12:30 分类：Algorithms

高斯消元法的主要思想是通过对系数矩阵进行行变换，将方程组的系数矩阵由对称矩阵变为三角矩阵，从而达到消元的目的，最后通过回代逐个获得方程组的解。

高斯消元法的实现简单，主要由两个步骤组成：第一个步骤就是通过选择主元，逐行消元，最终形成方程组系数矩阵的三角矩阵形式；第二步骤就是逐步回代的过程，最终矩阵的对角线上的元素就是方程组的解。

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