【原創】開源Math.NET基礎數學類庫使用(06)直接求解線性方程組

               本博客所有文章分類的總目錄：【總目錄】本博客博文總目錄-實時更新 

開源Math.NET基礎數學類庫使用總目錄：【目錄】開源Math.NET基礎數學類庫使用總目錄

前言

　　在前幾篇關於Math.NET的博客中(見上面鏈接)，主要是介紹了Math.NET中主要的數值功能，並進行了簡單的矩陣向量計算例子，接着使用Math.NET的矩陣等對象，對3種常用的矩陣數據交換格式的讀寫。一方面可以了解Math.NET的使用，另一方面以后也可以直接讀取和保存數據為這兩種格式，給大家的科研或者工作帶來便利。接下來的文章將繼續對Math.NET的功能進行講解和演示，並附帶一些數學方面的基礎知識。畢竟很多人沒有精力去研究Math.NET，那我就把我的研究心得一一寫出來，方便后來人。

1.數值分析與線性方程

　　數值分析的基本含義與特點：引用

　　數值分析(numerical analysis)是研究分析用計算機求解數學計算問題的數值計算方法及其理論的學科，是數學的一個分支，它以數字計算機求解數學問題的理論和方法為研究對象。為計算數學的主體部分。數值分析有如下特點：
1·面向計算機
2·有可靠的理論分析
3·要有好的計算復雜性
4·要有數值實驗
5.要對算法進行誤差分析　　

　　數值分析的主要內容：插值法，函數逼近，曲線擬和，數值積分，數值微分，解線性方程組的直接方法，解線性方程組的迭代法，非線性方程求根，常微分方程的數值解法。

　　所以我們今天要解決的就是數值分析的一個很小但又最常接觸到的部分，方程(組)的求解。方程(組)的求解有2種方法，一種是直接求解，一種是迭代求解。直接方法比較好理解，相當與使用公式進行直接計算，結果也比較精確。而另外一種是迭代方法。從最初的猜測,迭代方法逐次近似形式收斂只在極限的精確解。

　　正如上面所說，方程（組）的形式很多，不同的形式以及實際要求的精度都可以使用不同的方法，而本文主要介紹最簡單，也是最常見的線性方程的直接求解方法。

2.Math.NET解線性方程源碼分析

　　本文首先對Math.NET中求解線性方程的相關源碼進行分析，這樣大家碰到了問題，也可以更好的查找源碼解決，或者進行擴展，實現自己的一些特殊需求。Math.NET中MathNet.Numerics.LinearAlgebra.Factorization命名空間下有一個泛型接口：ISolver<T>，就是解AX = B類型的線性方程的接口類型。該接口功能很多，看看下面的接口源代碼和注釋(本人進行了簡單的翻譯)，就很清楚了。

using System;

namespace MathNet.Numerics.LinearAlgebra.Factorization
{
    /// <summary>形如AX = B的線性方程組求解的接口類型</summary>
    /// <typeparam name="T">泛型參數，支持類型有：double, single, <see cref="Complex"/>, 
    /// 以及 <see cref="Complex32"/>.</typeparam>
    public interface ISolver<T> where T : struct, IEquatable<T>, IFormattable
    {
        /// <summary>求解AX=B的線性方程組</summary>
        /// <param name="input">右矩陣<c>B</c>.</param>
        /// <returns>返回求解結果矩陣 <c>X</c>.</returns>
        Matrix<T> Solve(Matrix<T> input);

        /// <summary>求解AX=B的線性方程組</summary>
        /// <param name="input">右矩陣 <c>B</c>.</param>
        /// <param name="result">求解結果矩陣, <c>X</c>.</param>
        void Solve(Matrix<T> input, Matrix<T> result);

        /// <summary>求解AX=b的線性方程組</summary>
        /// <param name="input">等式的右邊向量 <c>b</c>.</param>
        /// <returns>返回求解結果的左邊向量 , <c>x</c>.</returns>
        Vector<T> Solve(Vector<T> input);

        /// <summary>求解AX=b的線性方程組</summary>
        /// <param name="input">等式的右邊向量 <c>b</c>.</param>
        /// <param name="result">求解結果矩陣, <c>x</c>.</param>
        void Solve(Vector<T> input, Vector<T> result);
    }
}

 由於求解線性方程組主要用到了矩陣的分解，Math.NET實現了5種矩陣分解的算法：LU，QR，Svd，Evd，Cholesky。而GramSchmidt是繼承QR的，每一個都是實現了ISolver<T>接口，因此就可以直接使用矩陣的分解功能，直接進行線性方程組的求解。為了方便，我們舉一個LU的源碼例子，簡單的看看源碼的基本情況：

public abstract class LU<T> : ISolver<T> where T : struct, IEquatable<T>, IFormattable
{
    static readonly T One = BuilderInstance<T>.Matrix.One;

    readonly Lazy<Matrix<T>> _lazyL;
    readonly Lazy<Matrix<T>> _lazyU;
    readonly Lazy<Permutation> _lazyP;

    protected readonly Matrix<T> Factors;
    protected readonly int[] Pivots;

    protected LU(Matrix<T> factors, int[] pivots)
    {
        Factors = factors;
        Pivots = pivots;

        _lazyL = new Lazy<Matrix<T>>(ComputeL);
        _lazyU = new Lazy<Matrix<T>>(Factors.UpperTriangle);
        _lazyP = new Lazy<Permutation>(() => Permutation.FromInversions(Pivots));
    }

    Matrix<T> ComputeL()
    {
        var result = Factors.LowerTriangle();
        for (var i = 0; i < result.RowCount; i++)
        {
            result.At(i, i, One);
        }
        return result;
    }

    /// <summary>Gets the lower triangular factor.</summary>
    public Matrix<T> L
    {
        get { return _lazyL.Value; }
    }
    /// <summary>Gets the upper triangular factor.</summary>
    public Matrix<T> U
    {
        get { return _lazyU.Value; }
    }
    /// <summary>Gets the permutation applied to LU factorization.</summary>
    public Permutation P
    {
        get { return _lazyP.Value; }
    }
    /// <summary>Gets the determinant of the matrix for which the LU factorization was computed.</summary>
    public abstract T Determinant { get; }
    public virtual Matrix<T> Solve(Matrix<T> input)
    {
        var x = Matrix<T>.Build.SameAs(input, input.RowCount, input.ColumnCount);
        Solve(input, x);
        return x;
    }    
    public abstract void Solve(Matrix<T> input, Matrix<T> result);
    
    public virtual Vector<T> Solve(Vector<T> input)
    {
        var x = Vector<T>.Build.SameAs(input, input.Count);
        Solve(input, x);
        return x;
    }    
    public abstract void Solve(Vector<T> input, Vector<T> result);    
    public abstract Matrix<T> Inverse();
}

　　大家可能會注意到，上面是抽象類，這和Math.NET的實現是有關的。最終都是實現相應版本的Matrix矩陣，然后實現對應版本的類型的分解方法。下面例子會介紹具體使用，大家有疑問，可以拿着源碼和例子，調試一番，知道上面的2個實現過程，就比較簡單了

3.Math.NET求解線性方程的實例

 　　假設下面是要求解的線性方程組(Ax=b):

5*x + 2*y  - 4*z = -7
3*x -  7*y + 6*z = 38
4*x + 1*y + 5*z = 43

 測試代碼，由於求解的方法很多，只列舉了幾種，其他的以此類推：

var formatProvider = (CultureInfo) CultureInfo.InvariantCulture.Clone();
formatProvider.TextInfo.ListSeparator = " ";

//先創建系數矩陣A 
var matrixA = DenseMatrix.OfArray(new[,] {{5.00, 2.00, -4.00}, {3.00, -7.00, 6.00}, {4.00, 1.00, 5.00}});

//創建向量b
var vectorB = new DenseVector(new[] {-7.0, 38.0, 43.0});

// 1.使用LU分解方法求解
var resultX = matrixA.LU().Solve(vectorB);
Console.WriteLine(@"1. Solution using LU decomposition");
Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));

// 2.使用QR分解方法求解
resultX = matrixA.QR().Solve(vectorB);
Console.WriteLine(@"2. Solution using QR decomposition");
Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));

// 3. 使用SVD分解方法求解
matrixA.Svd().Solve(vectorB, resultX);
Console.WriteLine(@"3. Solution using SVD decomposition");
Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));

// 4.使用Gram-Shmidt分解方法求解
matrixA.GramSchmidt().Solve(vectorB, resultX);
Console.WriteLine(@"4. Solution using Gram-Shmidt decomposition");
Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));

// 5.驗證結果，就是把結果和A相乘，看看和b是否相等
var reconstructVecorB = matrixA*resultX;
Console.WriteLine(@"5. Multiply coefficient matrix 'A' by result vector 'x'");
Console.WriteLine(reconstructVecorB.ToString("#0.00\t", formatProvider));

結果如下： 

1. Solution using LU decomposition
DenseVector 3-Double
3.00
1.00
6.00

2. Solution using QR decomposition
DenseVector 3-Double
3.00
1.00
6.00

3. Solution using SVD decomposition
DenseVector 3-Double
3.00
1.00
6.00

4. Solution using Gram-Shmidt decomposition
DenseVector 3-Double
3.00
1.00
6.00

5. Multiply coefficient matrix 'A' by result vector 'x'
DenseVector 3-Double
-7.00
38.00
43.00

　　今天的用法就到此為止，請關注博客，后續還有更多的分析和使用文章。

4.資源

　　如果本文資源或者顯示有問題，請參考 本文原文地址：http://www.cnblogs.com/asxinyu/p/4285245.html