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

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

• Have you tried discretizing the beam? The solution will be probably as good as telling MMA to solve partial differential equations. – Gypaets Oct 8 '16 at 21:40
• @mns This could be a way forward but as the equation is fourth order it will require good interpolation. What form of discretization do you think would work? – Hugh Oct 8 '16 at 21:45
• It depends on how exact are the values you need. I've posted an answer with some code I wrote doing an elasticity course but perhaps it's to simple for your needs. – Gypaets Oct 9 '16 at 14:05

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;
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.

• Thank you for this. I will look into it. It is not clear to me how you are achieving the boundary conditions. Also, mathematica has a generalised form for eigen systems so it is not necessary to take the inverse of the stiffness matrix. Many thanks for proving a solution. – Hugh Oct 9 '16 at 14:18
• @Hugh The BCs are in the flexibility matrix. With other BCs the flexibility matrix would be different. – Gypaets Oct 9 '16 at 14:27
• Your solution has the same problem as mine, try e.g. With[{i = 2}, Plot[{u@x, Subtract @@ eqn} /. {\[Lambda] -> 1/eiValues[[i]], u -> Interpolation[Transpose[{xElementMiddle, eigenForms[[i]]}]]} // Evaluate, {x, L/2/n, L - L/2/n}]] – xzczd Oct 10 '16 at 3:53
• @xzczd I don't know if the error could be shrunk further, but I've often heard that only the values of the first (or last- it depends on the algorithm) two or three eigenfrequencies can be trusted. In my case I don't expect the dynamic matrix to be very well conditioned, as every matrix I've encountered outside college ;). The Subtract @@ eqn term in our answers seems to be almost equal even using different solving methods, so probably the beam is just ill conditioned. – Gypaets Oct 10 '16 at 11:17

# 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. NDSolveFiniteDifferenceDerivative-based solution

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

fdd = NDSolveFiniteDifferenceDerivative;

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.

• Thanks for code. I have tried it and note that the error seems large. If I do Table[Plot[Evaluate[{ u[x], -\[Lambda] m[x] u[x] + D[EI[x] D[u[x], {x, 2}], {x, 2}] } /. {\[Lambda] -> val[[n]], u -> fun[[n]]}], {x, 0, L}], {n, 1, 10}] then I see that the eigenvector is typically much smaller than the residual. Am I doing something wrong? – Hugh Oct 9 '16 at 13:23
• @Hugh I don't think you've done anything wrong, seems the eigen functions given by EigenNDSolve isn't that accurate (or can this be a bug of the package? ). Adjustment for Npoly and precision seem not to help. But I should point out (at least first several) eigenvalues returned by EigenNDSolve are in good agreements with those you found with ParametricNDSolve. – xzczd Oct 9 '16 at 14:39
• I agree thy look correct. Would it be possible to pull out a mass and stiffness matrix from NDSolve using the discretization there? We could then use the standard eigensolver to get frequencies and mode shapes. – Hugh Oct 9 '16 at 14:46
• @Hugh To be honest, due to my limited knowledge for eigenvalue related problems, I don't know how to set boundary conditions in the differential matrix... – xzczd Oct 9 '16 at 14:50
• I agree that would be a challenge – Hugh Oct 9 '16 at 14:51