How to find the characteristic vibration frequencies of a system with any number of masses and springs using eigendecomposition?

Question

I have created a program to solve any system of three springs and two masses using Mathematica's Eigendecomposition functions. My goal is to solve a generalization of the spring system presented in the book Mathematical Methods in the Physical Sciences, 3rd Edition by Mary L. Boas on page 185 on the section on applications of diagonalization.

Here is what I have so far.

GetCoefficients[k1_, k2_, k3_] := 
 Coefficient[Expand[k1*x^2 + k2*(x - y)^2 + k3*y^2], #] & /@ {x^2, 
   x y, y^2}
SpringFunction[k1_, k2_, k3_, m1_, m2_] := 
 Module[{PotentialEnergyMatrix = {{GetCoefficients[k1, k2, k3][[1]], 
      GetCoefficients[k1, k2, k3][[2]]/
      2}, {GetCoefficients[k1, k2, k3][[2]]/2, 
      GetCoefficients[k1, k2, k3][[3]]}}, 
   KineticEnergyMatrix = {{m1, 0}, {0, m2}}}, <|
   "Potential Energy V" -> Expand[k1*x^2 + k2*(x - y)^2 + k3*y^2], 
   "Kinetic Energy Matrix T" -> MatrixForm@KineticEnergyMatrix, 
   "V matrix" -> MatrixForm@PotentialEnergyMatrix, 
   "\!\(\*SuperscriptBox[\(T\), \(-1\)]\)V" -> 
    MatrixForm[Inverse@KineticEnergyMatrix . PotentialEnergyMatrix], 
   "Eigenvalues" -> 
    Eigenvalues[Inverse@KineticEnergyMatrix . PotentialEnergyMatrix], 
   "Eigenvectors" -> 
    MatrixForm@
     Eigenvectors[
      Inverse@KineticEnergyMatrix . PotentialEnergyMatrix], 
   "\[Lambda]" -> 
    "m*\!\(\*FractionBox[SuperscriptBox[\(\[Omega]\), \(2\)], \
\(k\)]\)", 
   "\!\(\*SubscriptBox[\(\[Omega]\), \(1\)]\)" -> \[Sqrt](1/m k*
       Eigenvalues[
         Inverse@KineticEnergyMatrix . PotentialEnergyMatrix][[1]]), 
   "\!\(\*SubscriptBox[\(\[Omega]\), \(2\)]\)" -> \[Sqrt](1/m k*
       Eigenvalues[
         Inverse@KineticEnergyMatrix . PotentialEnergyMatrix][[2]])|>]

For the same code in a Mathematica notebook, see Mathematica Cloud Notebook

My goal is the following:

• expand it to any number of springs and masses. For example, if the user wanted to enter 7 springs and 6 masses.

• Make a function that shows steps

Answer 1

There are two parts to your question. The first part is how do you form the mass and stiffness matrices for a dynamic system. The second part is how do you solve the eigenvalue problem to get the natural frequencies and mode shapes. I will answer the the second part first and then the first part.

Solving eigenvalue problems for dynamic system.

I am going to do this in a general way, not just for your chain of masses and springs. All mass and stiffness matrices, however complicated the dynamic system, are symmetric positive definite. The symmetric part follows from linearity and work concepts. The positive definite part follows from the fact that kinetic and potential energies must always be positive.

First we make a module that will generate symmetric positive definite matrices.

ClearAll[pdMatrix];
pdMatrix[n_] := Module[{d, a},
  d = DiagonalMatrix[RandomReal[{0, 1}, n]];
  a = RandomReal[{-1, 1}, {n, n}];
  Transpose[a] . d . a]

Now we generate a mass and stiffness matrix and solve the eigenvalue problem using Eigensystem. This gives us eigenvalues (natural frequencies) and eigenvectors (mode shapes).

n = 5; (* Number of degrees of freedom *)
mm = pdMatrix[n];
kk = pdMatrix[n];
{vals, vecs} = Eigensystem[{kk, mm}];

The vals are the eigenvalues and these are the natural frequencies squared in radians per second. The vecs are the eigenvectors. The important property of the eigenvectors is that they can be used to diagonalize the mass and stiffness matrices simultaneously. Thus pre and post multiplying the mass and stiffness matrices by the eigenvectors produces diagonal matrices; I demonstrate this below. (Note the vectors are in row form so they have to be transposed to get them to the usual configuration).

m1 = vecs . mm . Transpose[vecs] // Chop; MatrixForm[m1]
k1 = vecs . kk . Transpose[vecs] // Chop; MatrixForm[k1]

Generating the mass and stiffness matrix for a chain of springs and masses.

Generating mass and stiffness matrices must be done for all dynamic systems. See here for a finite element example also look up NDEigensystm in help.

For your case of a chain of masses and springs there are three cases. 1) Springs attach the end masses to a fixed foundation, 2) the end masses are free and 3) one end mass is attached by a spring to a foundation and the other end is free. I will do the first case. Note that in this case there are n+1 springs for the n masses.

First I generate springs of random stiffness and masses of random mass. I then form the potential energy and kinetic energy summations. These are then turned into matrices by taking the derivatives with respect to displacement and velocity.

n = 9;(* Number of degrees of freedom *)
ClearAll[x, v];
k = RandomReal[{0, 100}, n + 1];
m = RandomReal[{0, 10}, n];
PE = Plus @@ 
   Join[{1/2 k[[1]] x[1]^2}, 
    Table[1/2 k[[j]] (x[j] - x[j - 1])^2, {j, 2, n}], {1/
      2 k[[n + 1]] x[n]^2}];
KE = Plus @@ Table[1/2 m[[j]] v[j]^2, {j, 1, n}];
{a, kk} = 
 CoefficientArrays[Table[D[PE, x[j]], {j, n}], Table[x[j], {j, n}]] //
   Normal; MatrixForm[kk]
{a, mm} = 
 CoefficientArrays[Table[D[KE, v[j]], {j, n}], Table[v[j], {j, n}]] //
   Normal; MatrixForm[mm]

As expected the stiffness matrix is tridiagonal and the mass matrix diagonal. We can now solve the eigenvalue problem and look at the results

{vals, vecs} = Eigensystem[{kk, mm}];
Column[Reverse@
  Table[ListLinePlot[{Transpose[vecs[[j]]], -Transpose[vecs[[j]]]}, 
    AspectRatio -> 1/4, PlotRange -> All, ImageSize -> 10 72, 
    PlotLabel -> 
     "Mode = " <> ToString[n - j + 1] <> " Frequency (Hz) = " <> 
      ToString[Sqrt[vals[[j]]]/(2 π)]], {j, n}]]

The pictures show the eigenvectors or mode shapes at the two extreme ends of their vibration cycle.

I hope that helps.