博客
关于我
强烈建议你试试无所不能的chatGPT,快点击我
开源Math.NET基础数学类库使用(06)数值分析之线性方程组直接求解
阅读量:6243 次
发布时间:2019-06-22

本文共 6184 字,大约阅读时间需要 20 分钟。

原文:

开源Math.NET基础数学类库使用系列文章总目录:

  1. 

  2. 

  3.

  4.

  5.

  6.

  7. 

  8.

  9.

10.

11.

12.

13.

14.

后续继续更新中。。如文章链接打开有误,请关注博客,因为文章正在编辑修改中,所有已经列出的目录都将在1个月之内发表。

前言

  在前几篇关于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类型的线性方程的接口类型。该接口功能很多,看看下面的接口源代码和注释(本人进行了简单的翻译),就很清楚了。

1 using System; 2  3 namespace MathNet.Numerics.LinearAlgebra.Factorization 4 { 5     /// 形如AX = B的线性方程组求解的接口类型 6     /// 
泛型参数,支持类型有:double, single,
, 7 /// 以及
.
8 public interface ISolver
where T : struct, IEquatable
, IFormattable 9 {10 ///
求解AX=B的线性方程组11 ///
右矩阵
B
.12 ///
返回求解结果矩阵
X
.
13 Matrix
Solve(Matrix
input);14 15 ///
求解AX=B的线性方程组16 ///
右矩阵
B
.17 ///
求解结果矩阵,
X
.18 void Solve(Matrix
input, Matrix
result);19 20 ///
求解AX=b的线性方程组21 ///
等式的右边向量
b
.22 ///
返回求解结果的左边向量 ,
x
.
23 Vector
Solve(Vector
input);24 25 ///
求解AX=b的线性方程组26 ///
等式的右边向量
b
.27 ///
求解结果矩阵,
x
.28 void Solve(Vector
input, Vector
result);29 }30 }

 由于求解线性方程组主要用到了矩阵的分解,Math.NET实现了5种矩阵分解的算法:LU,QR,Svd,Evd,Cholesky。而GramSchmidt是继承QR的,每一个都是实现了ISolver<T>接口,因此就可以直接使用矩阵的分解功能,直接进行线性方程组的求解。为了方便,我们举一个LU的源码例子,简单的看看源码的基本情况:

1 public abstract class LU
: ISolver
where T : struct, IEquatable
, IFormattable 2 { 3 static readonly T One = BuilderInstance
.Matrix.One; 4 5 readonly Lazy
> _lazyL; 6 readonly Lazy
> _lazyU; 7 readonly Lazy
_lazyP; 8 9 protected readonly Matrix
Factors;10 protected readonly int[] Pivots;11 12 protected LU(Matrix
factors, int[] pivots)13 {14 Factors = factors;15 Pivots = pivots;16 17 _lazyL = new Lazy
>(ComputeL);18 _lazyU = new Lazy
>(Factors.UpperTriangle);19 _lazyP = new Lazy
(() => Permutation.FromInversions(Pivots));20 }21 22 Matrix
ComputeL()23 {24 var result = Factors.LowerTriangle();25 for (var i = 0; i < result.RowCount; i++)26 {27 result.At(i, i, One);28 }29 return result;30 }31 32 ///
Gets the lower triangular factor.33 public Matrix
L34 {35 get { return _lazyL.Value; }36 }37 ///
Gets the upper triangular factor.38 public Matrix
U39 {40 get { return _lazyU.Value; }41 }42 ///
Gets the permutation applied to LU factorization.43 public Permutation P44 {45 get { return _lazyP.Value; }46 }47 ///
Gets the determinant of the matrix for which the LU factorization was computed.48 public abstract T Determinant { get; }49 public virtual Matrix
Solve(Matrix
input)50 {51 var x = Matrix
.Build.SameAs(input, input.RowCount, input.ColumnCount);52 Solve(input, x);53 return x;54 } 55 public abstract void Solve(Matrix
input, Matrix
result);56 57 public virtual Vector
Solve(Vector
input)58 {59 var x = Vector
.Build.SameAs(input, input.Count);60 Solve(input, x);61 return x;62 } 63 public abstract void Solve(Vector
input, Vector
result); 64 public abstract Matrix
Inverse();65 }

  大家可能会注意到,上面是抽象类,这和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

 测试代码,由于求解的方法很多,只列举了几种,其他的以此类推:

1 var formatProvider = (CultureInfo) CultureInfo.InvariantCulture.Clone(); 2 formatProvider.TextInfo.ListSeparator = " "; 3  4 //先创建系数矩阵A  5 var matrixA = DenseMatrix.OfArray(new[,] {
{
5.00, 2.00, -4.00}, {
3.00, -7.00, 6.00}, {
4.00, 1.00, 5.00}}); 6 7 //创建向量b 8 var vectorB = new DenseVector(new[] {-7.0, 38.0, 43.0}); 9 10 // 1.使用LU分解方法求解11 var resultX = matrixA.LU().Solve(vectorB);12 Console.WriteLine(@"1. Solution using LU decomposition");13 Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));14 15 // 2.使用QR分解方法求解16 resultX = matrixA.QR().Solve(vectorB);17 Console.WriteLine(@"2. Solution using QR decomposition");18 Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));19 20 // 3. 使用SVD分解方法求解21 matrixA.Svd().Solve(vectorB, resultX);22 Console.WriteLine(@"3. Solution using SVD decomposition");23 Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));24 25 // 4.使用Gram-Shmidt分解方法求解26 matrixA.GramSchmidt().Solve(vectorB, resultX);27 Console.WriteLine(@"4. Solution using Gram-Shmidt decomposition");28 Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));29 30 // 5.验证结果,就是把结果和A相乘,看看和b是否相等31 var reconstructVecorB = matrixA*resultX;32 Console.WriteLine(@"5. Multiply coefficient matrix 'A' by result vector 'x'");33 Console.WriteLine(reconstructVecorB.ToString("#0.00\t", formatProvider));

结果如下: 

1. Solution using LU decompositionDenseVector 3-Double3.001.006.002. Solution using QR decompositionDenseVector 3-Double3.001.006.003. Solution using SVD decompositionDenseVector 3-Double3.001.006.004. Solution using Gram-Shmidt decompositionDenseVector 3-Double3.001.006.005. Multiply coefficient matrix 'A' by result vector 'x'DenseVector 3-Double-7.0038.0043.00

  今天的用法就到此为止,请关注博客,后续还有更多的分析和使用文章。

4.资源

  源码下载:

  如果本文资源或者显示有问题,请参考 : 

本博客还有大量的.NET开源技术文章,您可能感兴趣: 

1.:

2.:

3.

4.

5.

6. 

转载地址:http://zxsia.baihongyu.com/

你可能感兴趣的文章
Hexo Next底部powered by的logo栏更改以及注意事项(附官方文档,文末有福利链)
查看>>
我是如何进入阿里巴巴的-面向春招应届生Java面试指南(七)
查看>>
Android Studio 打包生成的 apk 安装包装到手机上闪退
查看>>
Mybatis技术内幕:初始化之加载 mybatis-config
查看>>
mysql与pymysql
查看>>
Fastlane(二):结构
查看>>
vue高阶组件
查看>>
Android消息机制Handler源码分析
查看>>
HashMap JDK1 8源码
查看>>
2018年互联网圈,程序员圈竟然脱单的这么多?
查看>>
数据结构:解读哈夫曼树
查看>>
重新学习web后端开发-003-了解http请求
查看>>
230. Kth Smallest Element in a BST
查看>>
关于Apt注解实践与总结【包含20篇博客】
查看>>
PAT A1004
查看>>
学习webpack4 - 第三方库的使用
查看>>
PAT A1052
查看>>
vue工程全局设置ajax的等待动效
查看>>
Sublime Text3插件安装及问题处理
查看>>
前端如何通过Nginx代理做到跨域访问API接口
查看>>