矩阵分解是解线性方程组的一个方法,比较容易。
有LU,QR,SVD等等方法。在c#中用MathNet操作起来也是比较简单。
https://blog.csdn.net/c914620529/article/details/50393223/
https://www.cnblogs.com/asxinyu/p/4285245.html
https://www.cnblogs.com/asxinyu/p/Dotnet_Opensource_MathNet_Basic_Statistics_10.html
程序如下:
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using MathNet.Numerics.LinearAlgebra.Double;
namespace SimpleCalcModule
{
static public class MatrixCalc
{
//LU
public static bool calc1(int n, double[,] A, double[] b, ref double[] X)
{
if (n < 3)
return false;
if (b.Length != n)
return false;
if (A.GetLength(0) != n && A.GetLength(1) != n)
return false;
try
{
var matrixA = new DenseMatrix(n, n);
var vectorB = new DenseVector(n);
for (int i = 0; i < n; i++)
{
vectorB[i] = b[i];
X[i] = .0;
for (int j = 0; j < n; j++)
{
matrixA[i, j] = A[i, j];
}
}
//1.LU
var resultX = (DenseVector)matrixA.LU().Solve(vectorB);
X = resultX.ToArray();
if (X.Length == n)
return true;
else
return false;
}
catch
{
return false;
}
}
//QR
public static bool calc2(int n, double[,] A, double[] b, ref double[] X)
{
if (n < 3)
return false;
if (b.Length != n)
return false;
if (A.GetLength(0) != n && A.GetLength(1) != n)
return false;
try
{
var matrixA = new DenseMatrix(n, n);
var vectorB = new DenseVector(n);
for (int i = 0; i < n; i++)
{
vectorB[i] = b[i];
X[i] = .0;
for (int j = 0; j < n; j++)
{
matrixA[i, j] = A[i, j];
}
}
//1.LU
var resultX = (DenseVector)matrixA.QR().Solve(vectorB);
X = resultX.ToArray();
if (X.Length == n)
return true;
else
return false;
}
catch
{
return false;
}
}
//SVD
public static bool calc3(int n, double[,] A, double[] b, ref double[] X)
{
if (n < 3)
return false;
if (b.Length != n)
return false;
if (A.GetLength(0) != n && A.GetLength(1) != n)
return false;
try
{
var matrixA = new DenseMatrix(n, n);
var vectorB = new DenseVector(n);
for (int i = 0; i < n; i++)
{
vectorB[i] = b[i];
X[i] = .0;
for (int j = 0; j < n; j++)
{
matrixA[i, j] = A[i, j];
}
}
//1.LU
var resultX = (DenseVector)matrixA.Svd().Solve(vectorB);
X = resultX.ToArray();
if (X.Length == n)
return true;
else
return false;
}
catch
{
return false;
}
}
//Gram-Shmidt
public static bool calc4(int n, double[,] A, double[] b, ref double[] X)
{
if (n < 3)
return false;
if (b.Length != n)
return false;
if (A.GetLength(0) != n && A.GetLength(1) != n)
return false;
try
{
var matrixA = new DenseMatrix(n, n);
var vectorB = new DenseVector(n);
for (int i = 0; i < n; i++)
{
vectorB[i] = b[i];
X[i] = .0;
for (int j = 0; j < n; j++)
{
matrixA[i, j] = A[i, j];
}
}
//1.LU
var resultX = (DenseVector)matrixA.GramSchmidt().Solve(vectorB);
X = resultX.ToArray();
if (X.Length == n)
return true;
else
return false;
}
catch
{
return false;
}
}
}
}
测试函数
if (commOperate == 99)
{
//double[,] a = { { 8.1, 2.3, -1.5, 6.1 }, { 0.5, -6.23, 0.87, 2.3 }, { 2.5, 1.5, 10.2, 1.8 } };
//double[,] a = { { 2, -1, 3, 1 }, { 4, 2, 5, 4 }, { 1, 2, 0, 7 } };
double[,] A = { { .0, .0, 1.0 }, { 50.0, .0, 1.0 }, { .0, 50.0, 1.0 } };
int n = A.GetLength(0);//数组a的第一维长度,即行数,3
double[] b1 = new double[n];//存放解的数组,初始值为0
double[] b2 = new double[n];//存放解的数组,初始值为0
b1[0] = 37.5;
b1[1] = 62.5;
b1[2] = 37.5;
double[] X1 = new double[n];
double[] X2 = new double[n];
double[] X3 = new double[n];
double[] X4 = new double[n];
b2[0] = 38.75;
b2[1] = 38.75;
b2[2] = 63.75;
double[] X5 = new double[n];
double[] X6 = new double[n];
double[] X7 = new double[n];
double[] X8 = new double[n];
bool rst1 = SimpleCalcModule.MatrixCalc.calc1(n, A, b1, ref X1);
bool rst2 = SimpleCalcModule.MatrixCalc.calc2(n, A, b1, ref X2);
bool rst3 = SimpleCalcModule.MatrixCalc.calc3(n, A, b1, ref X3);
bool rst4 = SimpleCalcModule.MatrixCalc.calc4(n, A, b1, ref X4);
bool rst5 = SimpleCalcModule.MatrixCalc.calc1(n, A, b2, ref X5);
bool rst6 = SimpleCalcModule.MatrixCalc.calc2(n, A, b2, ref X6);
bool rst7 = SimpleCalcModule.MatrixCalc.calc3(n, A, b2, ref X7);
bool rst8 = SimpleCalcModule.MatrixCalc.calc4(n, A, b2, ref X8);
double x4 = 50;
double y4 = 50;
double Tx1 = X1[0] * x4 + X1[1] * y4 + X1[2];
double Tx2 = X2[0] * x4 + X2[1] * y4 + X2[2];
double Tx3 = X3[0] * x4 + X3[1] * y4 + X3[2];
double Tx4 = X4[0] * x4 + X4[1] * y4 + X4[2];
double Ty1 = X5[0] * x4 + X5[1] * y4 + X5[2];
double Ty2 = X6[0] * x4 + X6[1] * y4 + X6[2];
double Ty3 = X7[0] * x4 + X7[1] * y4 + X7[2];
double Ty4 = X8[0] * x4 + X8[1] * y4 + X8[2];
}
多谢,亲爱的美美。