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:
However when the script is not working as intended with only miniscule changes to the parameters, I get something like this:
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! :)