Why does matrix inversion cause extremely large values to appear in the inverted matrix?

I'm working on making an fea program where I have a stiffness matrix thats 158 x 158.

When I do \ to find the displacements, I get a sparse vector with several values around 1E16.

I inspected my sitffness matrix and found that most values were in the range of 0.001 to .500

So I inverted the matrix to see what kind of values would appear. I found that my last 4 rows and columns of the inverted stiffness matrix

are all around 1E16.

What could be causing this behavior?

Answer 1

[Chunru]

Consider ax = b where all variables are scaler, what will happen if a is very samll? Then x wil be large

In linear system Ax = b, if A is not well conditioned, in a sense that its inverse can be very large, x = A\b may have some very large elements. You can use "cond(A)" to check if A is well conditioned or not in matlab.

[raymond bryant]

So since my cond() came back as 2.4E32 this matrix poorly conditioned.

Thank you.

[Walter Roberson]

Yes, that is very poorly conditioned. Did you check its rank()?

[raymond bryant]

If I remember correctly, the rank was 152.

[Walter Roberson]

Rank 152 on a matrix 158 is a singular matrix which has 0 as the mathematical determinant. Calculation of the matrix inverse would involve division by that 0.

When you have a singular stiffness matrix if I recall correctly that signifies that you have too many degrees of freedom, that something is not anchored properly. Like if you build a bridge but forget to fasten it to the ends.

[raymond bryant]

Yes, that's what I've been looking into, it may be an issue with my rotation degrees of freedom. What should the rank be for a non-singular matrix? Is it 158?

[Walter Roberson]

Yes, a non-singular matrix must be square and full rank.

When you use the \ operator if you have more rows than columns then the system is over determined and the operator will return the best-fit coefficients with lowest norm if I recall correctly.

[Chunru]

A non-singulare square matrix is full rank. The rank should be size of the matrix.

Answer 2

[John D'Errico]

As others have said, a singular stiffness matrix typically means one of several things in a Finite Element problem.

It might be that you made a mistake in the creation of your global stiffness matrix. It happens. Since we cannot possibly know what you did wrong there, then only you can go through your code to fix that.

But perhaps the most common error made in FEA is to insufficiently constrain the problem. That is, if your entire system can freely rotate or translate in some direction as a RIGID body? Then the stiffness matrix will be singular. Remember that FEA works by minimizing the potential energy of the object, so if an enrgy free rotation of translatino exists, then there will be infinitely many solutions. You will then have a singularity. The clue is often a simple one, where the rank is typically one less than the size of the matrix.

Other possibilities are less simple, sometimes reducing to a divide by zero. For example, suppose I designed a very simple truss that connnects nodes a-b, b-c, c-d, but where nodes b and c are coincident? Then the truss member between b amd c will have zero length. Since the stiffness of that truss is found by dividing by the length of that truss member, you have a divide by zero. Other cases might be a 2-d FEA problem composed of triangular elements. But if some of those triangles are constrained to have all nodes in a straight line, then they have zero ara. And again, the stiffness matrix will divide by that area.

There are other ways I could surely come up with to create a numerically singular stiffness matrix, even though it is still mathematically well-posed. But they might involve systems where some elements were made of titanium, and others of polystyrene, so extremely varying stiffnesses.

Regardless, you CANNOT invert a singular stiffness matrix.

[raymond bryant]

Okay, thank you very much. I think part of my problem is I've only work on truss FEA systems and this is a beam. So I need to add in the rotation dof's for each node. I might also use a different method for adding the boundary conditions to my K matrix.