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]

Mathematica graphics

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}]

Mathematica graphics

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

  • $\begingroup$ Have you tried discretizing the beam? The solution will be probably as good as telling MMA to solve partial differential equations. $\endgroup$
    – Gypaets
    Oct 8, 2016 at 21:40
  • $\begingroup$ @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? $\endgroup$
    – Hugh
    Oct 8, 2016 at 21:45
  • $\begingroup$ 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. $\endgroup$
    – Gypaets
    Oct 9, 2016 at 14:05

2 Answers 2


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[#];
         (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}];
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];



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

enter image description here

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.

  • $\begingroup$ 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. $\endgroup$
    – Hugh
    Oct 9, 2016 at 14:18
  • $\begingroup$ @Hugh The BCs are in the flexibility matrix. With other BCs the flexibility matrix would be different. $\endgroup$
    – Gypaets
    Oct 9, 2016 at 14:27
  • $\begingroup$ 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}]] $\endgroup$
    – xzczd
    Oct 10, 2016 at 3:53
  • 1
    $\begingroup$ @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. $\endgroup$
    – Gypaets
    Oct 10, 2016 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}]

Mathematica graphics

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, 
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}]

Mathematica graphics

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.

  • 1
    $\begingroup$ 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? $\endgroup$
    – Hugh
    Oct 9, 2016 at 13:23
  • $\begingroup$ @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. $\endgroup$
    – xzczd
    Oct 9, 2016 at 14:39
  • $\begingroup$ 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. $\endgroup$
    – Hugh
    Oct 9, 2016 at 14:46
  • $\begingroup$ @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... $\endgroup$
    – xzczd
    Oct 9, 2016 at 14:50
  • $\begingroup$ I agree that would be a challenge $\endgroup$
    – Hugh
    Oct 9, 2016 at 14:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.