鸟语天空

追赶法求解方程组

作者：追风剑情 发布于：2018-10-26 22:27 分类：Algorithms

示例

using System;

namespace ConsoleApp2
{
    class Program
    {
        static void Main(string[] args)
        {
            //测试数据
            Console.WriteLine(@"求解方程组:
            17* x + 6.9*y + 0  *z = 200
            4 * x + 18 *y + 9.1*z = 300
            0 * x + 9  *y + 26 *z = 400
            ");
            float[] matrixA =
            {
                17, 6.9f, 0,
                4, 18, 9.1f,
                0, 9, 26
            };
            int DIM = 3;
            float[] dValue ={ 200, 300, 400 };

            float[] xValue;
            bool success = ThomasEquation.Resolve(matrixA, DIM, dValue, out xValue);
            if (success)
            {
                Console.WriteLine("方程组的解为: x={0}, y={1}, z={2}", 
                    xValue[0], xValue[1], xValue[2]);

                Console.WriteLine("验证解:");
                float x0 = xValue[0];
                float x1 = xValue[1];
                float x2 = xValue[2];
                for (int i = 0; i < DIM; i++)
                {
                    int j = i * DIM;
                    float m0 = matrixA[j];
                    float m1 = matrixA[j + 1];
                    float m2 = matrixA[j + 2];

                    float d = m0 * x0 + m1 * x1 + m2 * x2;

                    Console.WriteLine("{0}x{1} + {2}x{3} + {4}x{5} = {6}",
                        m0,x0, m1,x1, m2,x2, d);
                }
            }
            Console.Read();
        }
    }

    /**
     * 追赶法(托马斯法)求解多元一次方程组
     */
    public class ThomasEquation
    {
        /**
         * 追赶法求解方程
         * 算法原理: http://www.devacg.com/?post=590
         * @param   matrixA         系数矩阵
         * @param   DIM             矩阵维数
         * @param   dValue          Ax=d
         * @param   xValue          返回的解
         */
        public static bool Resolve(float[] matrixA, int DIM, float[] dValue, out float[] xValue)
        {
            xValue = new float[DIM];

            if (!CheckValidity(matrixA, DIM))
            {
                Console.WriteLine("系数矩阵不满足追赶法条件");
                return false;
            }

            float[] L = new float[DIM];
            float[] M = new float[DIM];
            float[] U = new float[DIM];
            float[] Y = new float[DIM];

            /* 消元，追的过程 */
            L[0] = matrixA[Index(0, 0, DIM)];
            U[0] = matrixA[Index(0, 1, DIM)] / L[0];
            Y[0] = dValue[0] / L[0];
            for(int i = 1; i < DIM; i++)
            {
                M[i] = matrixA[Index(i, i - 1, DIM)];
                L[i] = matrixA[Index(i, i, DIM)] - M[i] * U[i-1];
                if (i + 1 < DIM)
                {
                    U[i] = matrixA[Index(i, i + 1, DIM)] / L[i];
                }
                Y[i] = (dValue[i] - M[i] * Y[i - 1]) / L[i];
            }

            /* 回代求解，赶的过程 */
            xValue[DIM - 1] = Y[DIM - 1];
            for(int i = DIM - 2; i >= 0; i--)
            {
                xValue[i] = Y[i] - U[i] * xValue[i+1];
            }

            return true;
        }

        // 通过行列下标计算所在一维数组的索引
        private static int Index(int i, int j, int DIM)
        {
            return i * DIM + j;
        }

        // 检查系数矩阵是否满足追赶法求解条件
        private static bool CheckValidity(float[] matrixA, int DIM)
        {
            float a = 0;
            float b = 0;
            float c = 0;

            //检查：条件二
            a = matrixA[Index(0, 0, DIM)];
            c = matrixA[Index(0, 1, DIM)];
            if (Math.Abs(a) <= Math.Abs(c))
            {
                Console.WriteLine("(1)不满足求解条件二");
                return false;
            }

            a = matrixA[Index(DIM - 1, DIM - 1, DIM)];
            b = matrixA[Index(DIM - 1, DIM - 2, DIM)];
            if (Math.Abs(a) <= Math.Abs(b))
            {
                Console.WriteLine("(2)不满足求解条件二");
                return false;
            }

            for (int i = 1; i < DIM; i++)
            {
                a = matrixA[Index(i, i, DIM)];
                
                //检查: 条件一
                if (a == 0)
                {
                    Console.WriteLine("不满足求解条件一");
                    return false;
                }

                //检查: 条件三
                if (i < DIM - 1)
                {
                    b = matrixA[Index(i, i - 1, DIM)];
                    c = matrixA[Index(i, i + 1, DIM)];

                    if (Math.Abs(a) <= Math.Abs(b) + Math.Abs(c))
                    {
                        Console.WriteLine("不满足求解条件三");
                        return false;
                    }
                }
            }
            return true;
        }
    }
}

运行测试

参见 克洛脱(Crout)矩阵分解——LU分解

标签: Algorithms