How o use singular value decomposition function to find the natural frequency and mode shapes

I have a matrix R which is having a dependency on ω, I am interested in finding the natural frequency and mode shapes of this matrix using SVD how to do this?. The regular way is to take the determinant and find the cross over of ω at zero, which is cumbersome and takes more time.

R={{-634148.+78.5 ω^2,-9.77188*10^-10-1.71581*10^-14 ω^2,2.02044*10^-8+2.92574*10^-13 ω^2,494545. -61.2197 ω^2,-168584.+20.8681 ω^2,-4.5629+0.0000151949 ω^2,0,0,0,0,0,0,0,0,0,0,0,0,0.5,0,0.392699},{-9.77188*10^-10-1.71581*10^-14 ω^2,-1.01464*10^7+78.5 ω^2,1.32366*10^-7+3.44729*10^-13 ω^2,5.13413*10^6-39.7225 ω^2,5.0401*10^6-38.9937 ω^2,45.1142 -0.0000202285 ω^2,0,0,0,0,0,0,0,0,0,0,0,0,-0.707107,0,-1.74472*10^-8},{2.02044*10^-8+2.92574*10^-13 ω^2,1.32366*10^-7+3.44729*10^-13 ω^2,-5.1366*10^7+78.5 ω^2,-2.51028*10^6+3.83722 ω^2,4.60004*10^7-70.3001 ω^2,-54.5166-0.0000184728 ω^2,0,0,0,0,0,0,0,0,0,0,0,0,0.5,0,-1.1781},{494545. -61.2197 ω^2,5.13413*10^6-39.7225 ω^2,-2.51028*10^6+3.83722 ω^2,-3.14241*10^6+68.0354 ω^2,-44.667+0.000130837 ω^2,-167.801+0.00112853 ω^2,0,0,0,0,0,0,0,0,0,0,0,0,0.0000454151,0,-0.355056},{-168584.+20.8681 ω^2,5.0401*10^6-38.9937 ω^2,4.60004*10^7-70.3001 ω^2,-44.667+0.000130837 ω^2,-4.45509*10^7+87.9732 ω^2,73.5031 -0.0000932266 ω^2,0,0,0,0,0,0,0,0,0,0,0,0,-4.51929*10^-6,0,1.11558},{-4.5629+0.0000151949 ω^2,45.1142 -0.0000202285 ω^2,-54.5166-0.0000184728 ω^2,-167.801+0.00112853 ω^2,73.5031 -0.0000932266 ω^2,-1.62342*10^8+78.4998 ω^2,0,0,0,0,0,0,0,0,0,0,0,0,-1.83318*10^-6,0,-2.22145},{0,0,0,0,0,0,-1.2337*10^9+78.5 ω^2,0,0,-1.34797*10^9+85.7707 ω^2,-3.43845*10^8+21.8787 ω^2,-2.65699*10^8+16.9064 ω^2,0,0,0,0,0,0,-0.00222144,0,0},{0,0,0,0,0,0,0,-4.9348*10^9+78.5 ω^2,0,2.65789*10^9-42.2802 ω^2,-3.82643*10^9+60.8687 ω^2,-1.88146*10^9+29.9291 ω^2,0,0,0,0,0,0,-0.00444285,0,0},{0,0,0,0,0,0,0,0,-1.11033*10^10+78.5 ω^2,-1.30209*10^9+9.20574 ω^2,-9.94639*10^9+70.3207 ω^2,5.15779*10^9-36.4654 ω^2,0,0,0,0,0,0,0.00666423,0,0},{0,0,0,0,0,0,-1.34797*10^9+85.7707 ω^2,2.65789*10^9-42.2802 ω^2,-1.30209*10^9+9.20574 ω^2,-3.28329*10^9+117.828 ω^2,8.7716*10^6,1.65686*10^7,0,0,0,0,0,0,-0.0000101084,0,0},{0,0,0,0,0,0,-3.43845*10^8+21.8787 ω^2,-3.82643*10^9+60.8687 ω^2,-9.94639*10^9+70.3207 ω^2,8.7716*10^6,-1.31329*10^10+117.827 ω^2,3.32568*10^7,0,0,0,0,0,0,-0.0000405303,0,0},{0,0,0,0,0,0,-2.65699*10^8+16.9064 ω^2,-1.88146*10^9+29.9291 ω^2,5.15779*10^9-36.4654 ω^2,1.65686*10^7,3.32568*10^7,-1.17495*10^10+47.2269 ω^2,0,0,0,0,0,0,-0.0118894,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,-3.21934*10^14+0.313998 ω^2,2.54828*10^9+2.39864*10^-7 ω^2,-3.22467*10^9+4.19847*10^-7 ω^2,0,0,0,0,1.99999,688.251},{0,0,0,0,0,0,0,0,0,0,0,0,2.54828*10^9+2.39864*10^-7 ω^2,-1.26437*10^16+0.314 ω^2,3.28974*10^10-5.13839*10^-7 ω^2,0,0,0,0,-2.,-2390.39},{0,0,0,0,0,0,0,0,0,0,0,0,-3.22467*10^9+4.19847*10^-7 ω^2,3.28974*10^10-5.13839*10^-7 ω^2,-9.91289*10^16+0.314 ω^2,0,0,0,0,2.00001,3924.33},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6.1685*10^14+157. ω^2,-0.104144+1.51402*10^-16 ω^2,-0.0315765+6.35558*10^-15 ω^2,0,4.01771*10^-14,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.104144+1.51402*10^-16 ω^2,-5.55165*10^15+157. ω^2,-0.547232+2.3335*10^-14 ω^2,0,6.36878*10^-14,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.0315765+6.35558*10^-15 ω^2,-0.547232+2.3335*10^-14 ω^2,-1.54213*10^16+157. ω^2,0,3.03551*10^-14,0},{0.5,-0.707107,0.5,0.0000454151,-4.51929*10^-6,-1.83318*10^-6,-0.00222144,-0.00444285,0.00666423,-0.0000101084,-0.0000405303,-0.0118894,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,1.99999,-2.,2.00001,4.01771*10^-14,6.36878*10^-14,3.03551*10^-14,0,0,0},{0.392699,-1.74472*10^-8,-1.1781,-0.355056,1.11558,-2.22145,0,0,0,0,0,0,688.251,-2390.39,3924.33,0,0,0,0,0,0}}

• {u,v,w}=SingularValueDecomposition[R]; such that u.v.ConjugateTranspose[w] == R, and DiagonalMatrixQ[v]. Jul 23, 2020 at 12:03
• But there is an unknown in the matrix R that is ω, which I am interested in finding Jul 23, 2020 at 12:11
• Why is that a problem? The SVD works on symbolic matrices too. Jul 23, 2020 at 12:12
• It's not working, It is still running. Jul 23, 2020 at 12:36
• I was trying this line of thinking: You can write R as $R = M+W \omega^2$ using W = Map[Coefficient[#, ω^2] &, R, {2}]; M = R - ω^2 W;. We have both W and M poorly conditioned matrices, but W is definitely a singular matrix. We can then write (R - M).PseudoInverse[W] which should give us a mostly diagonal matrix - e.g have a look at (R - M).PseudoInverse[W] /.ω-> 200 // MatrixPlot - however I got a bit lost here. Jul 23, 2020 at 17:24

I think this is a numerical problem and I only have a partial solution. The solution to your problem is to find eigenvalues that make the determinate 0. I looked at the determinate using Table and then plotting the results

d = Table[{ω, Det[R]}, {ω, 0, 400, 10}];
ListPlot[d]


So there is a root around 270. This is the first eigenvalue but note the values of the determinate they are huge and this is the problem.

In an attempt to stabilise the determinant I divided by a large number and then I can plot.

Plot[Det[R/10^7], {ω, 0, 10000},
PlotRange -> {All, 10^9 {-1, 1}}]


The zero crossings are where the eigenvalues are located. The range is still huge. Trying to find the zero crossing using FindRoot fails due to the number size

FindRoot[Det[R/10^7], {ω, 267}]

General::munfl: 1.41268*10^-183 (-9.0127*10^-179) is too small to represent as a normalized machine number; precision may be lost.


So we need an expert to tell us how to normalise your numbers. Can you see a way of normalising your problem to bring the numbers into a narrower range? I note that the matrix R has values of 10^16.

Once you have got the number range down I suggest you look at getting the M and K matrices separately. Then you can use the generalised version of Eigensystem.

• Actually I am directly constructing the dynamic stiffness matrix using the variational energy approach where I have constructed Lagrangian. The minimization of the lagrangian directly results in the dynamic stiffness matrix. So I don't get the mass and stiffness matrix separately. Is there anything that I am doing wrong? Jul 23, 2020 at 17:16
• I don't think you are doing anything wrong. Looking at your R matrix it appears you can seperate the mass and stiffness matrix. If I have time I will show how. However, your real problem is the dynamic range. You can normalise by mass and time. You now know that the first eigenvalue is at about 270 (radians per second) so normalise time to make that about 1. You should be able to normalise the mass as well to make that about 1. There may be standard methods to normalise an eigen problem so that may be worth investigating. Perhaps post another question.
– Hugh
Jul 23, 2020 at 18:16
• Thanks I will take a look at how to normalize eigen value problem Jul 23, 2020 at 18:20