How do you find the eigenvalues of a PDE (Dynamic Euler-Bernoulli beam)?

Question

I am continuing to work on the vibration of a beam modeled by the Euler-Bernoulli equation. I have had some good answers to simulating the motion which may be found here. Now I wish to calculate the eigenvalues and vectors.

The equation is

$\frac{\partial ^2}{\partial x^2}\left(\text{EI}(x) \frac{\partial ^2v(t,x)}{\partial x^2}\right)+m(x) \frac{\partial ^2v(t,x)}{\partial t^2}=0$

where $v(t,x)$ is the displacement $\text{EI}(x)$ is the bending stiffness and $m(x)$ the mass per unit length.

The standard approach for calculating the eigenvalues is to assume a harmonic solution of the form

$v(t, x) = u(x) sin(t ω)$

which gives us the ordinary differential equation

$\frac{\partial ^2}{\partial x^2}\left(\text{EI}(x) \frac{\partial ^2u(x)}{\partial x^2}\right)-\lambda m(x) u(x)=0$

Here I have written $\lambda$ for $\omega ^2$ for simplicity.

When the bending stiffness and mass per unit length are constants there is an exact solution. However even for that case we cannot use NDEigensystem because that only supports second order equations.

Here is a minimal working example

ClearAll[m, EI, L];
L = 5;
m[x_] := 12.25 (-0.024 + (0.4 - 0.03 x)^2)
EI[x_] := 94.522 (-0.0005772 + (0.4 - 0.03 x)^4)
eqn = -λ m[x] u[x] + D[EI[x] D[u[x], {x, 2}], {x, 2}] == 0;
bc = {u[0] == 0, u'[0] == 0, u''[L] == 0, u'''[L] == 0};

Mathematica is not able to take the Laplace transform so I could not make progress in that direction.

What I have been able to do is to add an oscillating force to the boundary conditions and use ParametricNDSolve

bc1 = {u[0] == 0, u'[0] == 0, u''[L] == 0, u'''[L] == 1};
ps = ParametricNDSolve[Join[{eqn}, bc1], u, {x, 0, L}, λ];
Plot3D[Evaluate[u[λ][x] /. ps], {x, 0, L}, {λ, 0, 30},
  PlotRange -> All]

When $\lambda$ equals an eigenvalue the response should be infinite and this is reflected in the ridges. In particular plotting the displacement at the end of the beam

Plot[Evaluate[u[λ][L] /. ps], {λ, 0, 30}]

shows a nice resonance curve with infinite values at the eigenvalues. I have tried using FindRoot on the reciprocal of this curve, for example,

FindRoot[Evaluate[(1/u[λ][L] /. ps) == 0], {λ, 6.6} ]

and this does give a result together with many warnings from ParametricNDSolve.

I need a better method for finding eigenvalues and vectors. Please can you help? Thanks

Answer 1

The discrete formulation of the vibration movement of an elastic structure is

$$[M]\cdot \ddot x + [K]\cdot x=0$$

where $[M]$ is the mass and [K] is the stiffness Matrix. With the harmonic solution approach the equation can be written as

$$[K]^{-1}[M]-\lambda \cdot [I]=0$$

Length, stiffness and mass density:

ClearAll[m, EI, L];
L = 5;
m[x_] := 12.25 (-0.024 + (0.4 - 0.03 x)^2)
EI[x_] := 94.522 (-0.0005772 + (0.4 - 0.03 x)^4)

Spatial discretization (set number of elements and get middle coordinate of each):

n = 80;(*Elements number*)
dL = ConstantArray[L/n, n];
xElementMiddle = Range[L/(2*n), L, L/n];

Mass matrix $[M]$:

massDensityElement = m /@ xElementMiddle;
massElement = MapThread[Times, {massDensityElement, dL}];
massMatrix = DiagonalMatrix@massElement;

Flexibility matrix $[K]^{-1}$ (calculated with Maxwell-Mohr):

flexibilityMatrix[dl_, eI_] := Module[{n, maxI, minI, eiM},
  n = Length@dl;
  eiM = Map[
    (maxI = Max[#];
      minI = Min[#];
      Sum[dl[[k]]/(6*eI[[k]])*(
         (Plus @@ dl[[1 ;; minI]] - 
             If[k > 1, Plus @@ dl[[1 ;; (k - 1)]], 0])*
           (2*(Plus @@ dl[[1 ;; maxI]] - 
                If[k > 1, Plus @@ dl[[1 ;; (k - 1)]], 0]) + 
             Plus @@ dl[[1 ;; maxI]] - Plus @@ dl[[1 ;; k]])
          + (Plus @@ dl[[1 ;; minI]] - Plus @@ dl[[1 ;; k]])*
           (Plus @@ dl[[1 ;; maxI]] - 
             If[k > 1, Plus @@ dl[[1 ;; (k - 1)]], 0] + 
             2*(Plus @@ dl[[1 ;; maxI]] - Plus @@ dl[[1 ;; k]]))), {k,
         minI}]) &,
    Table[{i, j}, {i, n}, {j, n}], {2}];
  eiM
  ]
eiElement = EI /@ xElementMiddle;
flexMatrix = flexibilityMatrix[dL, eiElement];`

Calculating the eigenvalues results in almost the same solutions as in your plot:

dynamicMatrix = flexMatrix.massMatrix*1.;
eiValues = Eigenvalues[dynamicMatrix];
eiValues^-1

{0.05,1.01,6.5,23.7,63.3,...}

Eigenforms:

eigenForms = Eigenvectors@dynamicMatrix;
ListLinePlot[eigenForms[[1 ;; 4]]]

This approach only takes into account the displacement of the finite masses. You can, however, add the rotation terms to the mass and stiffness matrix to get better results.

Answer 2

1.Package-based solution

This problem can be solved by EigenNDSolve in this package:

{val, fun} = EigenNDSolve[{eqn, bc} // Flatten, {u}, {x, 0, 5}, λ];

plotEValues[val, {{0, 1000}, {-2, 2}}, estateList -> {fun}]

However, as pointed out by Hugh in the comment, the eigenfunctions returned by EigenNDSolve seems to contain significant error, anyway, (at least first several) eigenvalues returned by EigenNDSolve are in good agreements with those you found with ParametricNDSolve.

2. NDSolve`FiniteDifferenceDerivative-based solution

NDSolve`FiniteDifferenceDerivative can return "DifferentiationMatrix", which can be used for calculating eigenvalue and eigenvector with Eigensystem, let's make use of it:

fdd = NDSolve`FiniteDifferenceDerivative;

domain = {0, 5};
points = 100;
xorder = 4;
grid = Array[# &, points, domain];

MapThread[(toindex[#] = #2) &, {{0, 5, x}, {1, -1, 
   All}}];(*N[CGLGrid[domain,points],16];*)
tomatrix[o_] := {((u n_: 0) | Derivative[n_][u])[x_] :> 
   fdd[n, grid, "DifferenceOrder" -> o]["DifferentiationMatrix"][[toindex@x]], x :> grid}

eqmat = Subtract @@ eqn /. tomatrix@xorder;
bcmat = Subtract @@@ bc /. tomatrix@xorder;
(eqmat[[#]] = bcmat[[#]]) & /@ {1, 2, -1, -2};
Sort[Eigenvalues[{-eqmat /. λ -> 0, Coefficient[eqmat, λ]}]][[;; 20]]
(* {0.0500032, 1.04256, 6.76361, 24.5573, 65.491, 144.303, 279.418, 492.946,… *)

{value, vec} = 
  Function[a, With[{order = Ordering@a[[1]]}, #[[order]] & /@ a]]@
   Eigensystem[{eqmat /. λ -> 0, -Coefficient[eqmat, λ]}];

However, further check shows this method seems to suffer the same issue as the first method. Just observe the following:

Table[Plot[{u[x], Subtract @@ eqn} /. {λ -> value[[i]], 
     u -> Interpolation[{grid, vec[[i]]}\[Transpose]]} // Evaluate, {x, 0, L}], {i, 5}]

Modifying points or xorder doesn't help. I'm not sure about the reason. Can this be a sign that the equation itself is something special that can't solve by FDM or spectral expansion? Or something is wrong with the coding of mine and the author of the EigenNDSolve? I may look into this later, but now I'd like to go to bed.