How to work with Matrices that have a great difference in element size (e.g. 1 to 10e+15)?

Dear Mathworks-Community,

I have been trying for the past days to create a script that calculates the displacement of a suspension bridge using the Finite-Element-Method.

All has been well so far, but as you might already suspect I had to deal with matrices getting quite huge in dimensions, but also with huge numbers within said matrix. Especially since my main units are millimeters for a bridge, this yields great differences in magnitude when calculating the moment of interia (height^3) compared to the cross-sectional area (height^1) alone.

I have noticed that my results are sometimes all over the place, not making sense at all, for a given set of parameters. Then, if I change the input parameters only slightly, the results miraculously make sense. This happened for example when I input the starting angle of the main cable of the bridge. If I input 50.0001° its all fine, although if I input 50.001°, the results are going completely awry, no patterns to be seen.

The same thing happens when I change the length of my bridge from 200000mm to 20000mm. First, no sensibe results, but with the smaller length it just makes sense.

This has led me to believe that I might have a problem during the inversion of my matrix. As soon as specific elements in the matrix become big enough to matter, the results are all over the place.

In the beginning I started to calculate the displacements using the backslash-operator, however this method was out of the question quite early because it never yielded sensible results.

GrandDisplacementMatrix = GrandStiffnessMatrix \ GrandForceMatrix

Next up was the pinv-operator, which is less accurate but I can live with that. With this all was fine, but as I said before, it is only fine for some parameters.

GrandDisplacementMatrix = pinv(GrandStiffnessMatrix) * GrandForceMatrix

I am marking results as sensible when the y-coordinate of my displacement is in descending order and symmatric. For example:

118.9413
207.2895
259.9209
277.4011
259.9212
207.2900
118.9418

However when the script is not working as intended with only miniscule changes to the parameters, I get something like this:

-0.4346
0.3620
-0.5611
0.5137
-0.4465
0.9111
-0.4104

After playing around with the parameter I concluded that I have an issue with the handling of the matrix, perhaps the number sizes vary too much. Some entries in the stiffness matrix have a magnitude of 1e5, others have 1e15. As far as I know this should not be too much of an issue, but apparently it is.

My question is, what can I do to make it easier for MATLAB to handle this matrix? Can I do some pre-treatment of the matrix to make it more easily invertible? (I believe my RCOND is always around 1e-30, which is very very low.) I know this value is quite close to singularity, but I kept getting good results... For some parameters, that is.

Any tips would be fantastic! :)

Best Regards,

Lennart

[Sindar]

One option would be to transform to a more natural unit system, so that most values are close to 1. This especially helps with powers (using mm as base unit, mm^3 is 1 in unit system. Using meters, mm^3 is 1e-9)

Answer 1

[Thiago Henrique Gomes Lobato]

I don't believe you're facing problems due to the difference of magnitudes, but rather because your problem may be very ill-conditioned. You see, the pseudo-inverse uses SVD and then invert only the single values above a specific threshold. For a square well-conditioned matrix the result should actually be the same as using "\" or "inv". This means that your matrix is either rank-deficient or ill-conditioned with very small single values.

A way to improve this is using regularization. The regularization topic is very wide but one thing that you may try right aways is the Tikhonov regularization. This line should do it:

lambda = 0.1;
GrandDisplacementMatrix = (GrandStiffnessMatrix'*GrandStiffnessMatrix+lambda*eye(size(GrandStiffnessMatrix,2)))...
    \ (GrandStiffnessMatrix'*GrandForceMatrix);

You may have to play a little with lambda. Normally an initial approach would be some percentage of the maximum eigenvalue of your matrix, but it really depends of the problem. I'm sure that, using this regularization, you will get more consistent results, although the higher the lambda the more bias you will have in your result, so it should be as small as possible.