C#--Gauss消元之直接三角形分解法

Gauss消元之直接三角形分解法

基本介绍

若能通过 正交变换，将系数矩阵A分解为A=LU，其中L是单位下三角矩阵( 主对角线元素为1的下三角矩阵)，而U是上三角矩阵，则线性方程组Ax=b变为LUx=b，若令Ux=y，则线性方程组Ax=b的求解分为两个三角方程组的求解：

(1)求解Ly=b，得y；

(2)再求解Ux=y，即得方程组的解x；

因此三角分解法的关键问题在于 系数矩阵 A的LU分解。

矩阵能LU分解的充分条件

一般地，任给一个矩阵不一定有LU分解，下面给出一个 矩阵能LU分解的充分条件。

定理1 对任意矩阵

，若A的各阶顺序主子式均不为零，则A有唯一的Doolittle分解(或Crout分解)。

定理2 若矩阵

非奇异，且其LU分解存在，则A的LU分解是唯一的。

求解Ly=b与Ux=y的计算公式为

 核心代码实现，计算公式：

 public void Huanyuan4()
        {
            for (int r = 0; r < n; r++)
            {
                for (int j = r; j < n + 1; j++)
                {
                    double sum = 0;
                    for (int k = 0; k <= r - 1; k++)
                    {
                        sum += a[r, k] * a[k, j];
                    }
                    a[r, j] = a[r, j] - sum;//计算Uj,k
                }
                if (r < n)
                {
                    for (int i = r + 1; i < n; i++)
                    {
                        double sum = 0;
                        for (int k = 0; k <= r - 1; k++)
                        {
                            sum += a[i, k] * a[k, r];
                        }
                        a[i, r] = (a[i, r] - sum) / a[r, r];//计算Uir
                    }
                }
            }
            //调格式
            for (int r = 0; r < n; r++)
            {
                for (int j =0; j < n; j++)
                {
                    if (r == j && j > 0 && r > 0)
                    {
                        a[r, j - 1] = 0;
                        for (j = 1; j < r; j++)
                            a[r, j - 1] = 0;
                    }
                }
            }
        }

完整部分：

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace Guass4
{
    class Guass4
    {
        int n;

        public int N
        {
            get { return n; }
            set { n = value; }
        }  
        double[,] a;
        public double[,] A
        {
            get{  return a; }
            set {a = value; }
        }
        double[] x;
        public double[] X
        {
            get {return x; }
            set{x = value;}
        }
        public void Input()
        {
            Console.WriteLine("请输入阶数：");
            n = Convert.ToInt32(Console.ReadLine());
            a = new double[n, n + 1];
            x = new double[n];
            Console.WriteLine("请输入各行系数(','或' '隔开)：");
            for (int i = 0; i < n; i++)
            {
                string s = Console.ReadLine();
                string[] ss = s.Split(' ', ',');
                for (int j = 0; j < n + 1; j++)
                {
                    a[i, j] = Convert.ToDouble(ss[j]);
                }
            }
        }
        public void Huanyuan4()
        {
            for (int r = 0; r < n; r++)
            {
                for (int j = r; j < n + 1; j++)
                {
                    double sum = 0;
                    for (int k = 0; k <= r - 1; k++)
                    {
                        sum += a[r, k] * a[k, j];
                    }
                    a[r, j] = a[r, j] - sum;//计算Uj,k
                }
                if (r < n)
                {
                    for (int i = r + 1; i < n; i++)
                    {
                        double sum = 0;
                        for (int k = 0; k <= r - 1; k++)
                        {
                            sum += a[i, k] * a[k, r];
                        }
                        a[i, r] = (a[i, r] - sum) / a[r, r];//计算Uir
                    }
                }
            }
            //调格式
            for (int r = 0; r < n; r++)
            {
                for (int j =0; j < n; j++)
                {
                    if (r == j && j > 0 && r > 0)
                    {
                        a[r, j - 1] = 0;
                        for (j = 1; j < r; j++)
                            a[r, j - 1] = 0;
                    }
                }
            }
        }
        public void Output()
        {
            Console.WriteLine("方程系数为：");
            for (int i = 0; i < n; i++)
            {
                string s = null;
                for (int j = 0; j < n + 1; j++)
                {
                    s += string.Format("{0,8:f2}", a[i, j]);
                }
                Console.WriteLine(s);
            }
        }

        public void Huidai4()
        {
            for (int i=n-1;i>=0;i--)
            {
                double sum = 0;
                for(int j=i+1;j<n;j++)
                {
                    sum+= a[i, j] * x[j];
                }
                x[i] = (a[i, n] - sum) / a[i, i];
            }
            Console.WriteLine("\n方程组的解是：");
            for(int i=0;i<n;i++)
            {
                Console.WriteLine("x{0}={1}", i + 1, x[i]);
            }
        }
    }
    class Program
    {
        static void Main(string[] args)
        {
            Guass4 abc = new Guass4();
            abc.Input();
            abc.Output();
            abc.Huanyuan4();
            abc.Output();
            abc.Huidai4();
        }
    }
}

结果：

谢谢！！！！