我正在尝试使用数学。. 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]
数学。. NET和一个在线矩阵计算器，我发现大多数给我零(从迭代求解器)，NaN，或大的不正确的数字(从直接求解器)作为解决方案。
在数学上。. 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.00t", 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.00t", 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.00t", 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.00t", 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.00t", 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.00t", 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.00t", 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.00t", 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.00t", 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.00t", 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.00t", 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.00t", 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.00t", 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.00t", 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.00t", formatProvider));
Console.WriteLine();
Console.Read();
}
}
}
示例问题中的数字很可能是错误的，但我需要确保我正在使用数学。. NET正确，然后再继续。我用的这些解算器的例子是正确的吗?还有什么我可以尝试的，这些例子没有涵盖?
有限元分析例题(p.8，例1):
他们似乎在某个地方搞砸了单位，所以为了让我的矩阵匹配那里，我们必须使用以下输入:
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
还要注意，它们已经消除了一些在计算过程中应该自然消失的行和列。这些行和列仍然存在于我使用的矩阵
可以数学。。NET解任何矩阵?
不，它不能。具体地说，它不能解一个没有解的方程组，其他任何求解器也不能。
在这种情况下，你的矩阵A是奇异的，即它没有逆。这意味着你的方程组要么没有解，即它是不一致的，要么它有无限的解(参见数值方法介绍中的6.5节的例子)。一个奇异矩阵的行列式为零。您可以在mathnet中使用Determinant方法看到如下内容:
Console.WriteLine("Determinant {0}", matrixA.Determinant());
这给
Determinant 0
A是奇异的条件是它的行(或列)的线性组合为零。比如这里第二，第五，第八行的和是零。并不是只有这几行加起来等于零。(稍后您将看到另一个示例。实际上有三种不同的方法，这在技术上意味着这个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,]
注意行是零索引的。
通过将列向量b附加到A上，然后将相应的(0索引的)行相加，证明第1、第4和第7方程导致0=5形成增宽矩阵
Aaug = np.append(A,b,1)
Aaug[0,] + Aaug[3,] + Aaug[6,]
最后，即使你的矩阵不是奇异的，你仍然可以有一个数值不稳定的问题:在这种情况下，这个问题被称为病态的。检查矩阵的条件数，看看如何做到这一点(维基百科，np.linalg.cond(A), matrixA.ConditionNumber())。
你问题的最后两句话是你问题的根源:
还要注意，它们已经消除了一些在计算过程中应该自然消失的行和列。这些行和列仍然存在于我使用的矩阵中。
在你的例子问题中，你有固定的关节，反对在某些方向上的运动(称为边界条件)。有时在做有限元分析时，如果不根据这些边界条件从刚度矩阵和载荷矩阵中去掉适当的行和列，就会得到一个无法求解的系统，这里就是这种情况。
再试一次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);
回答你的问题，是的，你使用的方法是正确的，但是你解决的系统是错误的。修复你的输入，你应该得到你想要的输出。
我还应该指出，我建议使用直接求解器示例的原因是因为看起来您正在寻找精确的解决方案。迭代求解器仅通过近似解来节省计算时间。
不，它不能解奇异矩阵。但是其他代码也没有，因为这里没有解决方案。
对于您的特殊情况，A发布的矩阵是奇异的。大小是9×9，但排名是6。MATLAB报告1.9e17的条件。所以在你得到一个合理的答案之前，请检查你的刚度矩阵组成。也许你需要规范化你的矩阵，即提取E I系数，使数字从1e5降至1附近，这在数字上更容易接受。
供参考
如果你不喜欢Math.NET，或者你想验证代码使用纯c#。阅读这篇由James McCaffrey撰写的MSDN杂志文章，并使用列出的代码。
var A = new [,] { ... };
var b = new [] { ... };
var x = LU.SustemSolve(A,b);