Can Math.NET solve any matrix?
I am trying to use Math.NET to solve the following system:
Coefficient Matrix 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 }});
Result Vector b:
double[] loadVector = { 0, 0, 0, 5, 0, 0, 0, 0, 0 };
var vectorB = MathNet.Numerics.LinearAlgebra.Vector<double>.Build.Dense(loadVector);
I pulled these numbers from the leaf element example example, so the answer I'm expecting is based on this example:
[0.01316, 0, 0.0009199, 0.01316, -0.00009355, -0.00188, 0, 0, 0]
However, Math.NET and the online Matrix calculator , I found mostly to give me zeros (from iterative solvers), NaNs, or large wrong numbers (from direct solvers) as a solution.
In Math.NET, I tried to connect my matrices to the examples given, including:
Iterative example:
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();
}
}
}
Direct example:
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();
}
}
}
The numbers from the example problem may be wrong, but I need to be sure I am using Math.NET correctly before proceeding. Am I using these solver examples the way they are intended to be used? Is there anything else I can try so that these examples don't cover?
Example leaf element example (page 8, example 1):
They seem to confuse something with units, so we had to use the following input to match my matrix:
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
Also note that they have excluded some rows and columns, which should naturally disappear during the calculation. These rows and columns are still present in the matrix I am using
source to share
The last two sentences in your question are the source of your problem:
Also note that they have excluded some rows and columns, which should naturally disappear during the calculation. These rows and columns are still present in the matrix I am using.
In your example problem, you have joints that are locked against movement in certain directions (called boundary conditions). Sometimes in finite element analysis, if you do not remove the corresponding rows and columns from your stiffness matrix and load matrix according to these boundary conditions, you will end up with a system that cannot be resolved, which is the case here.
Try DirectSolver again:
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}});
and
double[] loadVector = { 0, 0, 5, 0, 0 };
var vectorB = MathNet.Numerics.LinearAlgebra.Vector<double>.Build.Dense(loadVector);
To answer your question, yes, you are using the methods correctly, but you are solving the wrong system. Correct your input and you should get the result you are looking for.
I should also point out that the reason I suggest using the direct solution example is because it looks like you are looking for the exact solution. Iterative solvers save computation time by simply approximating the solution.
source to share
Can Math.NET solve any matrix?
No, he can not. In particular, it cannot solve a system of equations that has no solution, and no other solver.
In this case, your matrix A
is singular, meaning it doesn't have the opposite. This means that your system of equations either has no solution, i.e. Is inconsistent or has infinite solutions (see Section 6.5 in Introduction to numeric methods, for example). The singular matrix has a zero determinant. You can see it in mathnet as follows using the method Determinant
:
Console.WriteLine("Determinant {0}", matrixA.Determinant());
This gives
Determinant 0
The condition for the singularity of A when the linear combination of its rows (or columns) is zero. For example, here the sum of the 2nd, 5th and 8th rows is zero. These are not the only strings that are concatenated to give zero. (You will see another example later. There are actually three different ways to do this, which technically means that this 9x9 matrix is "rank 6" and not "rank 9".).
Remember that all you do when you try to solve Ax=b
is solve a set of simultaneous equations. In two dimensions, you can have a system like
A = [1 1 b = [1 2 2], 2]
and solving this is tantamount to finding x0
and x1
such that
x0 + x1 = 1
2*x0 + 2*x1 = 2
There are infinite solutions here that satisfy x1 = 1 - x0
, i.e. along the line x0 + x1 = 1
. Alternative for
A = [1 1 b = [1 1 1], 2]
which is equivalent
x0 + x1 = 1
x0 + x1 = 2
there is clearly no solution, because we can subtract the first equation from the second to get 0 = 1
!
In your case 1st, 4th and 7th equations
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
Adding them together yields 0=5
, and therefore your system has no solution.
The easiest way to explore matrices is in an interactive environment like Matlab or R. Since Python is available in Visual Studio and it provides a Matlab-like environment using numpy. I've demonstrated this above with some Python code. I would recommend the visual studio Python tools that I have used successfully in Visual Studio 2012 and 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)
This is not an informative error message: LinAlgError: Singular matrix
. You can see what A
is singular for example by showing that the sum of 2nd, 5th and 8th rows is zero
A[1,]+A[4,]+A[7,]
noting that the lines are zero-indexed.
To demonstrate what the 1st, 4th, and 7th equations lead to 0=5
, form an expanded matrix by adding a columnar vector b
on A
and then adding the corresponding (0-indexed) rows together
Aaug = np.append(A,b,1) Aaug[0,] + Aaug[3,] + Aaug[6,]
Finally, even if your matrix is not singular, you may still have a numerically unstable problem: in this case, the problem is known as ill-conditioned. Check the condition number of the matrix to see how to do it ( the wikipedia , np.linalg.cond(A)
, matrixA.ConditionNumber()
).
source to share
No, it cannot solve singular matrices. But there is no other code there either, since there is no solution.
In your specific case, the posted matrix A
is the only one. Size 9 × 9, but grade 6. MATLAB
report condition 1.9e17
. So please check your stiffness matrix composition before you can expect a reasonable response. Perhaps you need to normalize your matrix i.e. Extract coefficients E I
to reduce numbers from 1e5
to 1
, which are more numerically acceptable.
FYI
If you don't like Math.NET
it or want to test the code, use pure c#
. Read this MSDN Magazine article by James McCaffrey and use the code in the list.
var A = new [,] { ... };
var b = new [] { ... };
var x = LU.SustemSolve(A,b);
source to share