Usage of Trigonometric Decomposition

三角分解解线性方程组

相比于高斯消元法，应用三角分解来解线性方程组在某些情况下能提高计算的效率。

三角分解解线性方程的公式推导

即，我们要解 LUx=b, 我们把这个方程拆为两个方程来解

$$\left\{ \begin{array}\ Ly=b\\ Ux=y \end{array} \right.$$

解方程组Ly=b的第一个方程为$y_{1}=b_{1}$，而第i个方程为

$$\sum_{j=1}^{i-1}(l_{ij}y_{ij}+y_{i}=b_{i})\Longrightarrow\\ \left\{ \begin{array}\ y_{1}=b_{1},\\ y_{i}=b_{i}-\sum_{j=1}^{i-1}l_{ij}y_{j},\quad i=1,2,3,\dots,n \end{array} \right.$$

将解得的y代入Ux=y，由高斯消元法回代公式可得

$$\left\{ \begin{array}\ x_{n}=\frac{y_{n}}{u_{nn}},\\ x_{i}=\frac{y_{i}-\sum_{j=i+1}^{n}u_{ij}x_{j}}{u_{ii}},\quad i=n-1,n-2,\dots,2,1 \end{array} \right.$$

---

一个LU分解的例子

$$A= \begin{pmatrix} 9&18&9&-27\\ 18&45&0&-45\\ 9&0&126&9\\ -27&-45&9&135 \end{pmatrix} \quad\quad\quad\quad\quad\quad b= \begin{pmatrix} 1\\ 2\\ 16\\ 8 \end{pmatrix}$$

现在对A进行分解

即可得

$$L= \begin{pmatrix} 1&\\ 2&1&\\ 1&-2&1\\ -3&1&\frac{2}{3}&1 \end{pmatrix} \quad\quad U= \begin{pmatrix} 9&18&9&-27\\ 0&9&-18&9\\ 0&0&81&54\\ 0&0&0&9 \end{pmatrix}\\ \quad\\ \quad\\ y= \begin{pmatrix} 1\\0\\15\\1 \end{pmatrix}$$

LU分解解线性方程组的算法分析

由矩阵LU分解公式可知，计算$u_{ij}$和$l_{ij}$各有

​ $\sum_{j=2}^{n-1}(j(n-j)+n-1=\sum_{1}^{n-1}j(n-j)=1/6n^{3}-1/6n$

个乘除法，所以矩阵分解共有$1/3n^3-1/3n$个乘除法，而求解下三角方程Ly=b有$1/2n(n-1)$个乘除法，求解上三角方程Ux=y有$1/2n(n-1)+n$个乘除法，综上，用LU分解来解线性方程组共有$\frac{1}{3}n^{2}+n^{2}-\frac{1}{3}n$个乘除法，与高斯消元法的运算次数一样。

但是从LU分解的形式上看，不难发现如果我们要解一系列系数矩阵相同但右端不同的方程组时，LU分解可以大大减少计算量。

编程实现LU分解

1. 程序如下

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

namespace Trigonometric_Decomposition
{
    class Program
    {
        static void Main(string[] args)
        {
           Matrix m= new Matrix(6, 6);
            m.LuDecomposition();
            m.DisPlay();
            Console.Read();
        }

        class Matrix
        {
            public int row = 0;
            public int column = 0;
            private float[,] matrix;
            \\矩阵的输入
            public Matrix(int row, int column)
            {
                this.row = row;
                this.column = column;
                Console.WriteLine($"Please enter a matrix with {this.row}rows,{this.column}colums");
                matrix = new float[row, column];
                try
                {
                    for (int i = 0; i < row; i++)
                    {
                        for (int j = 0; j < column; j++)
                        {
                            string mid = Console.ReadLine();
                            var rows = mid.Split(' ');
                            foreach (var item in rows)
                            {
                                matrix[i, j++] = float.Parse(item);
                            }

                        }

                    }
                    Console.WriteLine("\r\n");
                }
                catch (Exception e)
                {
                    Console.WriteLine($"something wrong!\n{e.Message}");
                }
            }
            \\输出矩阵
            public void DisPlay()
            {
                for (int i = 0; i < this.row; i++)
                {
                    for (int j = 0; j < this.column; j++)
                    {
                        Console.Write(matrix[i, j] + "   ");
                    }
                    Console.WriteLine("\r");
                }
            }
            \\矩阵的LU分解
            public void LuDecomposition()
            {
            //初始化第一行、第一列
                for (int i = 1; i < this.column; i++)
                {
                    matrix[i, 0] = matrix[i,0] / matrix[0, 0];
                }
                
                for (int i = 1; i < this.row; i++)
                {
                //先计算对角元，避免重复计算（放在第二个循环内会计算两次）
                    matrix[i, i] = matrix[i, i] - FigureTemp(i, i, i );
                    for (int j = i+1; j < this.column; j++)
                    {
                    //第i行i列的更新
                        matrix[i, j] = matrix[i, j] - FigureTemp(i, j, i);
                        matrix[j, i] = (matrix[j, i] - FigureTemp(i, j, i))/matrix[i,i];
                    }
                }
                
            }
            private float FigureTemp(int i,int j,int k)
            {
            //用于计算求和中间量
                float result = 0;
                for (int s = 0; s < k; s++)
                {
                    result += matrix[i, s] * matrix[s, j];
                }
                return result;
            }
        }
    }
}

1. 程序的验证

当我们输入矩阵

$$A= \begin{pmatrix} 9&18&9&-27\\ 18&45&0&-45\\ 9&0&126&9\\ -27&-45&9&135 \end{pmatrix}$$

将得到输出如下

经对比，与上面计算的一致，程序无误。