Mat*_* VS 3 .net math matrix linear-algebra mathnet-numerics
我正在尝试使用Math.NET来解决以下系统:
系数矩阵A:
var matrixA = DenseMatrix.OfArray(new[,] {
{ 20000, 0, 0, -20000, 0, 0, 0, 0, 0 },
{ 0, 666.66666666666663, 2000, 0, -666.66666666666663, 2000, 0, 0, 0 },
{ 0, 2000, 8000, 0, -2000, 4000, 0, 0, 0 },
{ -20000, 0, 0, 20666.66666666666663, 0, 2000, -666.66666666666663, 0, 2000 },
{ 0, -666.66666666666663, -2000, 0, 20666.66666666666663, -2000, 0, -20000, 0 },
{ 0, 2000, 4000, 2000, -2000, 16000, -2000, 0, 4000 },
{ 0, 0, 0, -666.66666666666663, 0, -2000, 666.66666666666663, 0, -2000 },
{ 0, 0, 0, 0, -20000, 0, 0, 20000, 0 },
{ 0, 0, 0, 2000, 0, 4000, -2000, 0, 7999.9999999999991 }});
结果矢量b:
double[] loadVector = { 0, 0, 0, 5, 0, 0, 0, 0, 0 };
var vectorB = MathNet.Numerics.LinearAlgebra.Vector<double>.Build.Dense(loadVector);
我从有限元分析示例问题中提取了这些数字,因此基于该示例我期望的答案是:
[0.01316, 0, 0.0009199, 0.01316, -0.00009355, -0.00188, 0, 0, 0]
但是,Math.NET和我发现的在线矩阵计算器主要给我零(来自迭代求解器),NaN或大的不正确数字(来自直接求解器)作为解决方案.
在Math.NET中,我尝试将我的矩阵插入到提供的示例中,包括:
迭代的例子:
namespace Examples.LinearAlgebra.IterativeSolversExamples
{
/// <summary>
/// Composite matrix solver
/// </summary>
public class CompositeSolverExample : IExample
{
public void Run()
{
// Format matrix output to console
var formatProvider = (CultureInfo)CultureInfo.InvariantCulture.Clone();
formatProvider.TextInfo.ListSeparator = " ";
// Solve next system of linear equations (Ax=b):
// 5*x + 2*y - 4*z = -7
// 3*x - 7*y + 6*z = 38
// 4*x + 1*y + 5*z = 43
// Create matrix "A" with coefficients
var matrixA = DenseMatrix.OfArray(new[,] { { 20000, 0, 0, -20000, 0, 0, 0, 0, 0 }, { 0, 666.66666666666663, 2000, 0, -666.66666666666663, 2000, 0, 0, 0 },
{ 0, 2000, 8000, 0, -2000, 4000, 0, 0, 0 }, { -20000, 0, 0, 20666.66666666666663, 0, 2000, -666.66666666666663, 0, 2000 },
{0, -666.66666666666663, -2000, 0, 20666.66666666666663, -2000, 0, -20000, 0 }, { 0, 2000, 4000, 2000, -2000, 16000, -2000, 0, 4000 },
{ 0, 0, 0, -666.66666666666663, 0, -2000, 666.66666666666663, 0, -2000 }, { 0, 0, 0, 0, -20000, 0, 0, 20000, 0 },
{0, 0, 0, 2000, 0, 4000, -2000, 0, 7999.9999999999991 }});
Console.WriteLine(@"Matrix 'A' with coefficients");
Console.WriteLine(matrixA.ToString("#0.00\t", formatProvider));
Console.WriteLine();
// Create vector "b" with the constant terms.
double[] loadVector = {0,0,0,5,0,0,0,0,0};
var vectorB = MathNet.Numerics.LinearAlgebra.Vector<double>.Build.Dense(loadVector);
Console.WriteLine(@"Vector 'b' with the constant terms");
Console.WriteLine(vectorB.ToString("#0.00\t", formatProvider));
Console.WriteLine();
// Create stop criteria to monitor an iterative calculation. There are next available stop criteria:
// - DivergenceStopCriterion: monitors an iterative calculation for signs of divergence;
// - FailureStopCriterion: monitors residuals for NaN's;
// - IterationCountStopCriterion: monitors the numbers of iteration steps;
// - ResidualStopCriterion: monitors residuals if calculation is considered converged;
// Stop calculation if 1000 iterations reached during calculation
var iterationCountStopCriterion = new IterationCountStopCriterion<double>(500000);
// Stop calculation if residuals are below 1E-10 --> the calculation is considered converged
var residualStopCriterion = new ResidualStopCriterion<double>(1e-10);
// Create monitor with defined stop criteria
var monitor = new Iterator<double>(iterationCountStopCriterion, residualStopCriterion);
// Load all suitable solvers from current assembly. Below in this example, there is user-defined solver
// "class UserBiCgStab : IIterativeSolverSetup<double>" which uses regular BiCgStab solver. But user may create any other solver
// and solver setup classes which implement IIterativeSolverSetup<T> and pass assembly to next function:
var solver = new CompositeSolver(SolverSetup<double>.LoadFromAssembly(Assembly.GetExecutingAssembly()));
// 1. Solve the matrix equation
var resultX = matrixA.SolveIterative(vectorB, solver, monitor);
Console.WriteLine(@"1. Solve the matrix equation");
Console.WriteLine();
// 2. Check solver status of the iterations.
// Solver has property IterationResult which contains the status of the iteration once the calculation is finished.
// Possible values are:
// - CalculationCancelled: calculation was cancelled by the user;
// - CalculationConverged: calculation has converged to the desired convergence levels;
// - CalculationDiverged: calculation diverged;
// - CalculationFailure: calculation has failed for some reason;
// - CalculationIndetermined: calculation is indetermined, not started or stopped;
// - CalculationRunning: calculation is running and no results are yet known;
// - CalculationStoppedWithoutConvergence: calculation has been stopped due to reaching the stopping limits, but that convergence was not achieved;
Console.WriteLine(@"2. Solver status of the iterations");
Console.WriteLine(monitor.Status);
Console.WriteLine();
// 3. Solution result vector of the matrix equation
Console.WriteLine(@"3. Solution result vector of the matrix equation");
Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));
Console.WriteLine();
// 4. Verify result. Multiply coefficient matrix "A" by result vector "x"
var reconstructVecorB = matrixA*resultX;
Console.WriteLine(@"4. Multiply coefficient matrix 'A' by result vector 'x'");
Console.WriteLine(reconstructVecorB.ToString("#0.00\t", formatProvider));
Console.WriteLine();
Console.Read();
}
}
}
直接示例:
namespace Examples.LinearAlgebraExamples
{
/// <summary>
/// Direct solvers (using matrix decompositions)
/// </summary>
/// <seealso cref="http://en.wikipedia.org/wiki/Numerical_analysis#Direct_and_iterative_methods"/>
public class DirectSolvers : IExample
{
/// <summary>
/// Gets the name of this example
/// </summary>
public string Name
{
get
{
return "Direct solvers";
}
}
/// <summary>
/// Gets the description of this example
/// </summary>
public string Description
{
get
{
return "Solve linear equations using matrix decompositions";
}
}
/// <summary>
/// Run example
/// </summary>
public void Run()
{
// Format matrix output to console
var formatProvider = (CultureInfo) CultureInfo.InvariantCulture.Clone();
formatProvider.TextInfo.ListSeparator = " ";
// Solve next system of linear equations (Ax=b):
// 5*x + 2*y - 4*z = -7
// 3*x - 7*y + 6*z = 38
// 4*x + 1*y + 5*z = 43
matrixA = DenseMatrix.OfArray(new[,] { { 20000, 0, 0, -20000, 0, 0, 0, 0, 0 }, { 0, 666.66666666666663, 2000, 0, -666.66666666666663, 2000, 0, 0, 0 },
{ 0, 2000, 8000, 0, -2000, 4000, 0, 0, 0 }, { -20000, 0, 0, 20666.66666666666663, 0, 2000, -666.66666666666663, 0, 2000 },
{0, -666.66666666666663, -2000, 0, 20666.66666666666663, -2000, 0, -20000, 0 }, { 0, 2000, 4000, 2000, -2000, 16000, -2000, 0, 4000 },
{ 0, 0, 0, -666.66666666666663, 0, -2000, 666.66666666666663, 0, -2000 }, { 0, 0, 0, 0, -20000, 0, 0, 20000, 0 },
{0, 0, 0, 2000, 0, 4000, -2000, 0, 7999.9999999999991 }});
Console.WriteLine(@"Matrix 'A' with coefficients");
Console.WriteLine(matrixA.ToString("#0.00\t", formatProvider));
Console.WriteLine();
// Create vector "b" with the constant terms.
double[] loadVector = { 0, 0, 0, 5000, 0, 0, 0, 0, 0 };
var vectorB = MathNet.Numerics.LinearAlgebra.Vector<double>.Build.Dense(loadVector);
Console.WriteLine(@"Vector 'b' with the constant terms");
Console.WriteLine(vectorB.ToString("#0.00\t", formatProvider));
Console.WriteLine();
// 1. Solve linear equations using LU decomposition
var resultX = matrixA.LU().Solve(vectorB);
Console.WriteLine(@"1. Solution using LU decomposition");
Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));
Console.WriteLine();
// 2. Solve linear equations using QR decomposition
resultX = matrixA.QR().Solve(vectorB);
Console.WriteLine(@"2. Solution using QR decomposition");
Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));
Console.WriteLine();
// 3. Solve linear equations using SVD decomposition
matrixA.Svd().Solve(vectorB, resultX);
Console.WriteLine(@"3. Solution using SVD decomposition");
Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));
Console.WriteLine();
// 4. Solve linear equations using Gram-Shmidt decomposition
matrixA.GramSchmidt().Solve(vectorB, resultX);
Console.WriteLine(@"4. Solution using Gram-Shmidt decomposition");
Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));
Console.WriteLine();
// 5. Verify result. Multiply coefficient matrix "A" by result vector "x"
var reconstructVecorB = matrixA*resultX;
Console.WriteLine(@"5. Multiply coefficient matrix 'A' by result vector 'x'");
Console.WriteLine(reconstructVecorB.ToString("#0.00\t", formatProvider));
Console.WriteLine();
// To use Cholesky or Eigenvalue decomposition coefficient matrix must be
// symmetric (for Evd and Cholesky) and positive definite (for Cholesky)
// Multipy matrix "A" by its transpose - the result will be symmetric and positive definite matrix
var newMatrixA = matrixA.TransposeAndMultiply(matrixA);
Console.WriteLine(@"Symmetric positive definite matrix");
Console.WriteLine(newMatrixA.ToString("#0.00\t", formatProvider));
Console.WriteLine();
// 6. Solve linear equations using Cholesky decomposition
newMatrixA.Cholesky().Solve(vectorB, resultX);
Console.WriteLine(@"6. Solution using Cholesky decomposition");
Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));
Console.WriteLine();
// 7. Solve linear equations using eigen value decomposition
newMatrixA.Evd().Solve(vectorB, resultX);
Console.WriteLine(@"7. Solution using eigen value decomposition");
Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));
Console.WriteLine();
// 8. Verify result. Multiply new coefficient matrix "A" by result vector "x"
reconstructVecorB = newMatrixA*resultX;
Console.WriteLine(@"8. Multiply new coefficient matrix 'A' by result vector 'x'");
Console.WriteLine(reconstructVecorB.ToString("#0.00\t", formatProvider));
Console.WriteLine();
Console.Read();
}
}
}
示例问题中的数字可能是错误的,但我需要确保我在继续之前正确使用Math.NET.我是否按照它们的使用方式使用这些求解器示例？还有什么我可以尝试这些例子不包含的内容吗？
他们似乎搞砸了某个地方的单位,所以为了让我的矩阵匹配,我们不得不使用以下输入:
Member A (mm^2) E (N/mm^2) I (mm^4) L (mm)
AB 600000000 0.0002 60000000 6
BC 600000000 0.0002 60000000 6
还要注意,它们已经消除了一些在计算过程中自然会消失的行和列.这些行和列仍然存在于我正在使用的矩阵中
Math.NET可以解决任何矩阵吗？
不,它不能.具体来说,它不能解决没有解的方程组,也不能解决任何其他求解器.
在这种情况下,您的矩阵A是单数的,即它没有反转.这意味着您的方程组没有解决方案,即它是不一致的,或者它具有无限的解决方案(参见数值方法简介中的 6.5节).奇异矩阵具有零行列式.你可以使用以下Determinant方法在mathnet中看到这个:
Console.WriteLine("Determinant {0}", matrixA.Determinant());
这给了
Determinant 0
A是奇异的条件是当其行(或列)的线性组合为零时.例如,这里第2行,第5行和第8行的总和为零.这些不是唯一一起添加到零的行.(稍后你会看到另一个例子.实际上有三种不同的方法,技术上意味着这个9x9矩阵是"等级6"而不是"等级9".).
请记住,当你试图解决时,你所做的Ax=b就是解决一组联立方程.在二维中,您可能有一个系统,如
A = [1 1 b = [1
2 2], 2]
而解决这等同于发现x0和x1这样
x0 + x1 = 1
2*x0 + 2*x1 = 2
这里有无限的解决方案x1 = 1 - x0,即沿着这条线x0 + x1 = 1.或者替代
A = [1 1 b = [1
1 1], 2]
这相当于
x0 + x1 = 1
x0 + x1 = 2
显然没有解决方案,因为我们可以从第二个方程中减去第一个方程式0 = 1!
在你的情况下,第1,第4和第7个方程是
20000*x0 -20000 *x3 = 0
-20000*x0 +20666.66666666666663*x3 +2000*x5 -666.66666666666663*x6 +2000*x8 = 5
-666.66666666666663*x3 -2000*x5 +666.66666666666663*x6 -2000*x8 = 0
将这些添加在一起会产生0=5,因此您的系统没有解决方案.
最简单的方法是在Matlab或R等交互式环境中探索矩阵.由于Python在Visual Studio中可用,它通过numpy提供类似Matlab的环境,我已经用Python中的一些代码演示了上述内容.我建议使用Visual Studio的Python工具,我在Visual Studio 2012和2013中都成功使用过.
# numpy is a Matlab-like environment for linear algebra in Python
import numpy as np
# matrix A
A = np.matrix ([
[ 20000, 0, 0, -20000, 0, 0, 0, 0, 0 ],
[ 0, 666.66666666666663, 2000, 0, -666.66666666666663, 2000, 0, 0, 0 ],
[ 0, 2000, 8000, 0, -2000, 4000, 0, 0, 0 ],
[ -20000, 0, 0, 20666.66666666666663, 0, 2000, -666.66666666666663, 0, 2000 ],
[ 0, -666.66666666666663, -2000, 0, 20666.66666666666663, -2000, 0, -20000, 0 ],
[ 0, 2000, 4000, 2000, -2000, 16000, -2000, 0, 4000 ],
[ 0, 0, 0, -666.66666666666663, 0, -2000, 666.66666666666663, 0, -2000 ],
[ 0, 0, 0, 0, -20000, 0, 0, 20000, 0 ],
[ 0, 0, 0, 2000, 0, 4000, -2000, 0, 7999.9999999999991 ]])
# vector b
b = np.array([0, 0, 0, 5, 0, 0, 0, 0, 0])
b.shape = (9,1)
# attempt to solve Ax=b
np.linalg.solve(A,b)
这失败了,提供了一条信息性错误消息:LinAlgError: Singular matrix.你可以看到它A是单数的,例如,显示第2行,第5行和第8行的总和为零
A[1,]+A[4,]+A[7,]
注意到行是零索引的.
为了证明第1,第4和第7个方程0=5通过将柱状向量附加b到其上A,然后将相应的(0索引)行添加到一起来形成增强矩阵
Aaug = np.append(A,b,1)
Aaug[0,] + Aaug[3,] + Aaug[6,]
最后,即使你的矩阵不是单数,你仍然可以有一个数字上不稳定的问题:在这种情况下,问题被称为病态.检查矩阵的条件数,看看如何做到这一点(维基百科,np.linalg.cond(A),matrixA.ConditionNumber()).
小智 4
您问题中的最后两句话是问题的根源：
另请注意，他们已经消除了一些在计算过程中自然消失的行和列。这些行和列仍然存在于我正在使用的矩阵中。
在您的示例问题中，您的关节是固定的，以防止在某些方向上的运动（称为边界条件）。有时在进行有限元分析时，如果不根据这些边界条件从刚度矩阵和载荷矩阵中删除适当的行和列，最终将得到一个无法求解的系统，这里就是这种情况。
使用以下命令再次尝试 DirectSolver：
var matrixA = DenseMatrix.OfArray(new[,] { {20000, 0, -20000, 0, 0}, {0, 8000 ,0, -2000 ,4000},
{-20000, 0, 20666.667 ,0, 2000}, {0, -2000 ,0, 20666.67, -2000},
{0, 4000 ,2000 ,-2000, 16000}});
和
double[] loadVector = { 0, 0, 5, 0, 0 };
var vectorB = MathNet.Numerics.LinearAlgebra.Vector<double>.Build.Dense(loadVector);
要回答你的问题，是的，你使用的方法是正确的，但你正在解决错误的系统。修正你的输入，你应该得到你正在寻找的输出。
我还应该指出，我建议使用直接求解器示例的原因是因为您似乎正在寻找精确的解决方案。迭代求解器仅通过近似解来节省计算时间。
| 归档时间： |
|
| 查看次数： |
4210 次 |
| 最近记录： |