我正在尝试使用 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
另请注意，他们已经消除了一些在计算过程中应该自然消失的行和列。这些行和列仍然存在于我正在使用的矩阵中