We assume that the homogeneous Euler-Bernoulli cantilever with constant cross-section length is 1, the material constant EI is 1 and ρS is 1.

enter image description here

According to the constraints, we can establish the following differential equations:

DSolve[{D[D[ω[x, t], {x, 2}], {x, 2}] + 
    D[ω[x, t], {t, 2}] == 0, ω[0, 0] == 
   0, ω[x, 0] == 0, ω[0, t] == 0, 
  D[ω[0, t], {t, 1}] == 
   0(*Constraints on simply supported ends*), 
  D[ω[1, t], {t, 2}] == 0, 
  D[ω[1, t], {t, 3}] == 
   0(*Constraints on fixed end*)}, ω[x, t], {x, t}]

But I still can't find the only analytical solution, while the textbook can determine the natural frequencies of each stage of the cantilever according to these conditions (In addition, I find that $\rho$ and $S$ cannot be replaced effectively in the following code).

Table[(((2. i - 1) π)/(
    2 l))^2 Sqrt[(\[DoubleStruckCapitalE]*\[DoubleStruckCapitalI])/(\
ρ*S)], {i, 1, 10}] /. {l -> 
   1., \[DoubleStruckCapitalE]*\[DoubleStruckCapitalI] -> 
   1., ρ -> 1., S -> 1}

How can I solve the previous differential equation to get similar frequency results?

When I use the relevant code of this post, I can't even get the numerical solution of the cantilever vibration:

tau = 10;
L = 1;
Elastic = 1;
Imoment = 1;
ρ = 1;
S = 1;
sol = NDSolveValue[{D[Elastic*Imoment*D[w[x, t], {x, 2}], {x, 2}] + 
     S*ρ*D[w[x, t], {t, 2}] == 0, w[0, t] == 0, 
   w[x, 0] == x^2/6 (3 - x), Derivative[0, 1][w][0, t] == 0, 
   Derivative[0, 2][w][L, t] == Derivative[0, 3][w][L, t] == 0}, 
  w, {x, 0, L}, {t, 0, tau}, 
  Method -> {"MethodOfLines", 
    "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 100},
     "SpatialDiscretization" -> {"TensorProductGrid", 
      "MaxPoints" -> 100, "MinPoints" -> 100, 
      "DifferenceOrder" -> 2}}, MaxSteps -> 10^6]
  • 4
    $\begingroup$ As per last part you need //. or ReplaceRepeated. Like this: Table[(((2. i - 1) \[Pi])/(2 l))^2 Sqrt[(\[DoubleStruckCapitalE]*\ \[DoubleStruckCapitalI])/(\[Rho]*S)], {i, 1, 10}] //. {l -> 1., \[Rho] -> 1., S -> 1., \[DoubleStruckCapitalE]*\[DoubleStruckCapitalI] -> 1.}. $\endgroup$
    – m0nhawk
    Commented Apr 5, 2020 at 3:12
  • 4
    $\begingroup$ As per first part, you have some strange boundary conditions, in equation you have 2nd derivative for time, but 3rd for boundary condition. I'd recommend to add equation and boundary conditions as-is. $\endgroup$
    – m0nhawk
    Commented Apr 5, 2020 at 3:17
  • 4
    $\begingroup$ Indeed, as m0nhawk pointed out, the constraints at x=1 are weird. Overall, your post is highly inconsistent: From the drawing, x=1 is a free end, but you call it fixed edn in your code. If x=1 is supposed to be the free end, then conditions like D[\[Omega][x, t], {x, 2}] == 0 /. x -> 1 or D[\[Omega][x, t], {x, 3}] == 0 /. x -> 1 should play a role. Also from the picture, the left end appears to be clamped which would result in D[\[Omega][x, t], {x, 1}] == 0 /. x -> 0 instead of D[ω[0, t], {t, 1}] == 0. $\endgroup$ Commented Apr 5, 2020 at 6:26
  • 6
    $\begingroup$ And since this a hyperbolic equation of order two in time, you would also need initial conditions for D[\[Omega][x, t], {t, 1}] /. t -> 0. But this all would be for simluating the dynamics of the beam. I have no idea why you would expect that you could read oof the eigenfrequencies for that. Maybe this is a specific trick from your textbook, but we cannot tell since you did not even cite it. (Frankly, I do not want to read it.) Please, for your own sake: Meet us half way and at least state the precise PDE that you want to solve in $\LaTeX$. $\endgroup$ Commented Apr 5, 2020 at 6:27
  • 3
    $\begingroup$ As is, this question does not seem to be Mathematica related: It is not about entering the PDE correctly (we cannot check because you did not provide the PDE in the first place), but to understand your textbook. This is not the purpose of this site. $\endgroup$ Commented Apr 5, 2020 at 6:33

3 Answers 3


To find eigenmodes, time is not involved. Following Wikipeida You just need to solve an ODE and not a PDE

enter image description here

This below is the code to reproduce the eigenmodes shown at the above page.

I tried to use Mathematica's NDEigesystem on this, but I could not make it do it, so did it by "hand".

The trick to finding eigenvalues is not to put all 4 boundaries conditions in at once, else Mathematica will return trivial solution of course, since all BC are zero and there is no load (it is free vibration). So we put 3 BC. You pick which 3 to put in and which one to keep out. Below I kept out the last one. (the third derivative on the free end).

After you get solution, it will have one constant of integration in it. Then by inspection, we can see now the condition to make the solution non-trivial. Now we set up an equation (it will be nonlinear) to solve for eigenvalues (Letting length of beam = 1).

Now we find the eigenvalues, going back to the ODE and plugging each eigenvalue at a time and putting back now all the 4 BC in. This results in corresponding eigenmode solution. To plot them all on same plot, I normalized them all to 1.

This is the final plot

enter image description here

This is what Wikipedia gives

enter image description here

(ps. For some reason, I had to flip the sign on some of the eigenmodes below to get same plot as Wikipedia. otherwise it was upside down. I have not found out why).

Here is a diagram of the BC (diagram thanks to this page)

For a cantilevered beam, the boundary conditions are as follows: (Thanks to this page)

Mathematica graphics

w(0)=0 . This boundary condition says that the base of the beam (at the wall) does not experience any deflection. w'(0)=0 . We also assume that the beam at the wall is horizontal, so that the derivative of the deflection function is zero at that point. w''(L)=0 . This boundary condition models the assumption that there is no bending moment at the free end of the cantilever. w'''(L)=0 . This boundary condition models the assumption that there is no shearing force acting at the free end of the beam.

Here is the code. This solve the ODE shown above. Notice it has 3 BC and not 4.

ClearAll[w, wHat, x, L];
pde = D[wHat[x], {x, 4}] == w^2 wHat[x];
leftEndBc = {wHat[0] == 0, wHat'[0] == 0};
rightEndBc = {wHat''[L] == 0};
sol = wHat[x] /. First@DSolve[{pde, leftEndBc, rightEndBc}, wHat[x], x]

This is what DSolve gives

$$ -\frac{c_1 e^{-\sqrt{w} x} \left(e^{L \sqrt{w}+2 \sqrt{w} x} \sin \left(L \sqrt{w}\right)+e^{2 L \sqrt{w}+\sqrt{w} x} \sin \left(\sqrt{w} x\right)-e^{L \sqrt{w}+2 \sqrt{w} x} \cos \left(L \sqrt{w}\right)-e^{2 L \sqrt{w}+\sqrt{w} x} \cos \left(\sqrt{w} x\right)+2 e^{L \sqrt{w}+\sqrt{w} x} \cos \left(L \sqrt{w}\right) \sin \left(\sqrt{w} x\right)-2 e^{L \sqrt{w}+\sqrt{w} x} \sin \left(L \sqrt{w}\right) \cos \left(\sqrt{w} x\right)+e^{2 L \sqrt{w}}+e^{L \sqrt{w}} \sin \left(L \sqrt{w}\right)+e^{L \sqrt{w}} \cos \left(L \sqrt{w}\right)-e^{2 \sqrt{w} x}+e^{\sqrt{w} x} \sin \left(\sqrt{w} x\right)+e^{\sqrt{w} x} \cos \left(\sqrt{w} x\right)\right)}{e^{2 L \sqrt{w}}+2 e^{L \sqrt{w}} \sin \left(L \sqrt{w}\right)-1} $$

This makes up equation to satisfy the missing BC at right end (the one left out)

  eq = FullSimplify[(D[sol, {x, 3}] /. x -> L)]

$$ \frac{2 c_1 w^{3/2} \left(\cos \left(L \sqrt{w}\right) \cosh \left(L \sqrt{w}\right)+1\right)}{\sin \left(L \sqrt{w}\right)+\sinh \left(L \sqrt{w}\right)} $$

Here is the "manual" step. Looking at the above, since the above is zero (it is the last B.C.), then for non-trivial solution we do not want $c_1=0$, then only other choice to make the above zero is that

$$ \cos \left(L \sqrt{w}\right) \cosh \left(L \sqrt{w}\right)+1=0 $$

Now NSolve is used to find eigenvalues $w$ (solutions of the above equation). We can use $L=1$ now. Then the equation to solve for frequencies is the above.

frequencies = w /. NSolve[1 + Cos[Sqrt[w]] Cosh[Sqrt[w]] == 0 && 0 < w < 200, w, Reals]


  {3.51602, 22.0345, 61.6972, 120.902, 199.86}

These are the first 5 frequencies. Now for each, we solve the ODE again.

pde = D[wHat[x], {x, 4}] == frequencies[[1]]^2 *wHat[x];
L = 1;
leftEndBc = {wHat[0] == 0, wHat'[0] == 0};
rightEndBc = {wHat''[L] == 0, wHat'''[L] == 0};
sol = NDSolveValue[{pde, leftEndBc, rightEndBc}, wHat, {x, 0, 1}];
max = Maximize[{Abs[sol[x]], 0 < x < 1}, x][[1]];
p1 = Plot[Callout[-sol[x]/max, "mode 1", 0.7], {x, 0, 1}, PlotStyle -> Blue]

pde = D[wHat[x], {x, 4}] == frequencies[[2]]^2 *wHat[x];
L = 1;
leftEndBc = {wHat[0] == 0, wHat'[0] == 0};
rightEndBc = {wHat''[L] == 0, wHat'''[L] == 0};
sol = NDSolveValue[{pde, leftEndBc, rightEndBc}, wHat, {x, 0, 1}];
max = Maximize[{Abs[sol[x]], 0 < x < 1}, x][[1]];
p2 = Plot[Callout[-sol[x]/max, "mode 2", .4, .4], {x, 0, 1}, PlotStyle -> Magenta]

pde = D[wHat[x], {x, 4}] == frequencies[[3]]^2 *wHat[x];
L = 1;
leftEndBc = {wHat[0] == 0, wHat'[0] == 0};
rightEndBc = {wHat''[L] == 0, wHat'''[L] == 0};
sol = NDSolveValue[{pde, leftEndBc, rightEndBc}, wHat, {x, 0, 1}];
max = Maximize[{Abs[sol[x]], 0 < x < 1}, x][[1]];
p3 = Plot[Callout[-sol[x]/max, "mode 3", .2], {x, 0, 1}, PlotStyle -> Brown]

pde = D[wHat[x], {x, 4}] == frequencies[[4]]^2 *wHat[x];
L = 1;
leftEndBc = {wHat[0] == 0, wHat'[0] == 0};
rightEndBc = {wHat''[L] == 0, wHat'''[L] == 0};
sol = NDSolveValue[{pde, leftEndBc, rightEndBc}, wHat, {x, 0, 1}];
max = Maximize[{Abs[sol[x]], 0 < x < 1}, x][[1]];
p4 = Plot[Callout[sol[x]/max, "mode 4", Below], {x, 0, 1}, PlotStyle -> {Thick, Green}]

 Show[{p1, p2, p3, p4}, PlotRange -> {{0, 1}, {-1.2, 1.2}}]

enter image description here

  • $\begingroup$ Thank you very much for your answer. When I give the cantilever a certain deformation at the initial moment, it still cannot be solved. $\endgroup$ Commented Apr 6, 2020 at 0:59
  • 3
    $\begingroup$ @PleaseCorrectGrammarMistakes The purpose of this answer is to give the natural frequencies for the beam, which it does and that was your question. It agrees with Wikipeia's β1 L/π = 0.59686..., β2 L/π = 1.49418..., β3 L/π = 2.50025..., β4 L/π = 3.49999..., ... on the above page if you compare them. (L=1 and E,I are all unity). I do not think we should start mixing this with other solutions. If you have specific question other than finding the frequencies, it will be better to post a separate question I think if you have something else and keep this just on finding the frequencies and modes. $\endgroup$
    – Nasser
    Commented Apr 6, 2020 at 1:12
  • 1
    $\begingroup$ Also notice that Wikipeida gives an expression for the mode shape as function of frequencies. I did not use that as I did not know how they found it. I used NDSolve to solve the ODE again using each specific frequency found earlier. The modes shapes agree. $\endgroup$
    – Nasser
    Commented Apr 6, 2020 at 1:16
  • 2
    $\begingroup$ @PleaseCorrectGrammarMistakes If you post the code you used when you say "When I give the cantilever a certain deformation at the initial moment, it still cannot be solved", I'l try to see what I can do. But please note, there is no time involved in the ODE to find mode shapes. You could have different boundary conditions OK, but time plays no role here when finding modal shapes. The boundary conditions I used are the ones to use for the cantilever beam. Other type of beam ofcourse is different story and need completely different solution. $\endgroup$
    – Nasser
    Commented Apr 6, 2020 at 3:19
  • $\begingroup$ I've posted questions here. But this time I just want to solve the vibration process, so it's still related to time. $\endgroup$ Commented Apr 6, 2020 at 7:55

I think Bill Watts answer is the correct approach. Just to add a bit more that can be generalised to other boundary conditions. Note that this is an eigenvalue problem. We are looking for frequencies that make the differential equation work. Thus we ought to form a determinant and look for roots.

Starting with Bill's equation but without the greek

X[x_] := c1 Cos[b x] + c2 Sinh[b x] + c3 Sin[b x] + c4 Cosh[b x]

We now introduce the boundary conditions for each end of the beam

  eqns = {
  X[0] == 0,  (* No displacment at x = 0 *)
  (D[X[x], x] /. x -> 0) == 0, (* No slope at x = 0 *)
  (D[X[x], {x, 2}] /. x -> L) == 0, (* No bending moment at x = L *)
  (D[X[x], {x, 3}] /. x -> L) == 0 (* No shear force at x = L *)

We can make this a matrix equation as follows. First we make a vector out of the unknown coefficients c1, c2, c3, c4 and then we use CoefficientArrays to get a matrix

vec = {c1, c2, c3, c4};
{rhs, mat} = Normal[CoefficientArrays[eqns, vec]];

enter image description here

We have a matrix equation in the form mat.vec == 0. So either we have the trivial solution vec == 0 or the determinant of the matrix must be zero. This is the standard eigenvalue problem. Taking the determinant we get

det = Det[mat] // Simplify

(*2 b^6 (1 + Cos[b L] Cosh[b L]) *)

Which is the same as Bills solution. We can see that because of the cosh term the roots are going to be approximately (2 n - 1) pi/2.

First we tidy up the determinant by combining L and b.

det1 = det/b^6 /. b -> Lb/L

Now we can get the roots of the determinate using the approximate solution as a starting point.

 roots = {b -> Lb/L} /. 
  Table[FindRoot[det1 == 0, {Lb, (2 n - 1) \[Pi]/2}], {n, 5}]

    (* {{b -> 1.8751/L}, {b -> 4.69409/L}, {b -> 7.85476/L}, {b -> 10.9955/
   L}, {b -> 14.1372/L}} *)

The eignevectors are found by determining the NullSpace of the matrix when each root is substituted in. We can do this as follows

evecs = Thread[vec -> #] & /@ (NullSpace[#][[1]] & /@ (mat /. roots));

Now I join all the root and corresponding eigenvector in one set of replacement rules.

valsVecs = Join[#[[1]], #[[2]]] & /@ Transpose[{roots, evecs}];

Here are the mode shapes

Plot[Evaluate[X[L x] /. valsVecs], {x, 0, 1}]

enter image description here

By applying other boundary conditions one can deal with all beam configurations.


Slightly different way for constant EI and ρS.

Start with the pde.

pde = D[EI*D[w[x, t], x, x], x, x] + \[Rho]S*D[w[x, t], t, t] == 0

Assume sinusoidal in time.

w[x_, t_] = X[x] Sin[ω t]

which simplifies the pde

(*EI*Derivative[4][X][x]*Sin[t*ω] - ρS*ω^2*X[x]*Sin[t*ω] == 0*)

dsol = DSolve[pde, X[x], x] // Flatten

Assign the constants to a new constant, β and solve for ω

ωRule = Solve[(ρS^(1/4) Sqrt[ω])/EI^(1/4) == β, ω] // Flatten

Clean up a bit

X[x_] = X[x] /. dsol /. {C[1] -> c1, C[2] -> c2, C[3] -> c3, 
     C[4] -> c4} /. ωRule // PowerExpand

(*c1 Cos[β x] + c2 E^(β (-x)) + c3 Sin[β x] + c4 E^(β x)*)

Get in Trig form.

X[x] // ExpToTrig // Collect[#, {Sinh[β x], Cosh[β x]}] &
(*c1 Cos[β x] + (c4 - c2) Sinh[β x] + (c2 + c4) Cosh[β x] + c3 Sin[β x]*)

X[x_] = % /. {c4 - c2 -> c2, c2 + c4 -> c4}
(*c1 Cos[β x] + c2 Sinh[β x] + c3 Sin[β x] + c4 Cosh[β x]*)

Apply some boundary conditions for a cantilever. The displacement is zero at the fixed end x = 0.

X[0] == 0

c4 = c4 /. Solve[%, c4][[1]]

Derivative of displacement is zero at the fixed end x = 0.

(D[X[x], x] /. x -> 0) == 0

c3 = c3 /. Solve[%, c3][[1]]

The Shear is zero at the free end.

(D[X[x], x, x, x] /. x -> L) == 0 //Simplify

c2 = c2 /. Solve[%, c2][[1]]

Bending Moment is zero at the free end.

(D[w[x, t], x, x] /. x -> L) == 0//Simplify

(β c1 Sin[t ω] (Cos[β L] Cosh[β L] + 1))/(
 Cos[β L] + Cosh[β L]) == 0

This gives us the condition for β.

Use βL so we can divide by L later.

F[βL_] = Cos[βL] Cosh[βL] + 1

Find the first five roots of βL with FindRoot. I didn't have much luck with NSolve. Plot to get starting values.

Plot[Evaluate[F[βL]], {βL, 0, 10}, PlotRange -> {-10, 10}]

enter image description here

n = 5;
a = FindRoot[F[α] == 0, {α, 2}]
βL[1] = α /. a // N;
a = FindRoot[F[α] == 0, {α, 4.8}]
βL[2] = α /. a // N;
For[i = 3, i <= n,
    a = FindRoot[ F[α] == 0, {α, 2 βL[i - 1] - βL[i - 2]}];
   βL[i] = α /. a // N; 

The first five frequencies.

Table[ω /. ωRule /. β -> βL[i]/L, {i, 5}]
{(3.51602 Sqrt[EI])/(L^2 Sqrt[ρS]), (22.0345 Sqrt[EI])/(
 L^2 Sqrt[ρS]), (61.6972 Sqrt[EI])/(L^2 Sqrt[ρS]), (
 120.902 Sqrt[EI])/(L^2 Sqrt[ρS]), (199.86 Sqrt[EI])/(
 L^2 Sqrt[ρS])}

You don't have to plug in numbers for the frequencies until the end this way. We have found all the constants except for c1 which is found from the initial conditions of the beam and orthogonality. The displacement of beam w = X[x] Sin[ω t] is found by an infinite sum over all the ω`s.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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