# Dynamic Euler–Bernoulli beam equation

I am trying to solve for the vibration of a Euler–Bernoulli beam. The equation is

$$\frac{\partial ^2u(t,x)}{\partial t^2}+\frac{\partial ^4u(t,x)}{\partial x^4}=0$$

For the boundary conditions I would like the displacement to be zero at the ends and with zero second derivative. This corresponds to pinned-pinned conditions. For time I will start with a displacement and no velocity.

In the future I would like to solve for a beam that is not uniform in thickness along the x-axis and for general initial conditions.

There is a similar problem in the NDEigensystem documentation here but this is for the standard wave equation which is only second order in space. However, I follow that example. First I define an initial displacement and try to solve the pde.

ClearAll[f];
f[x_] := x (1 - x)

tu = NDSolveValue[{
D[u[t, x], {t, 2}] + D[u[t, x], {x, 4}] == 0,
u[0, x] == f[x],
Derivative[1, 0][u][0, x] == 0,
DirichletCondition[u[t, x] == 0, True],
DirichletCondition[D[u[t, x], {x, 2}] == 0, True]
}, u, {t, 0, 1}, {x, 0, 1},
Method -> {"PDEDiscretization" -> "MethodOfLines"}];


This gives me the error

NDSolveValue::femcmsd: The spatial derivative order of the PDE may not exceed two.


Thus I proceed to supply two coupled differential equations one for displacement one for the second derivative (which is the bending moment). Thus I try to solve

tu = NDSolveValue[{
D[u[t, x], {t, 2}] + D[m[t, x], {x, 2}] == 0,
D[u[t, x], {x, 2}] == m[t, x],
u[0, x] == f[x],
Derivative[1, 0][u][0, x] == 0,
DirichletCondition[u[t, x] == 0, True],
DirichletCondition[m[t, x] == 0, True]
}, {u, m}, {t, 0, 1}, {x, 0, 1},
Method -> {"PDEDiscretization" -> "MethodOfLines"}];


However this also gives an error

NDSolveValue::ivone: Boundary values may only be specified for one independent variable. Initial values may only be specified at one value of the other independent variable.


I don't understand this error because I think I have done as asked... Can you help? Thanks

This post contains several code blocks, you can copy them easily with the help of importCode.

# Analytic Solution

The analytic solution can be obtained with LaplaceTransform and FourierSinCoefficient. First, make a Laplace transform on the equation and b.c.s and plug in the i.c.s:

Clear[f];
f[x_] = x (1 - x);

eqn = D[u[t, x], {t, 2}] + D[u[t, x], {x, 4}] == 0;
ic = {u[0, x] == f@x, Derivative[1, 0][u][0, x] == 0};
bc = {u[t, 0] == 0, u[t, 1] == 0, Derivative[0, 2][u][t, 0] == 0,
Derivative[0, 2][u][t, 1] == 0};
teqn = LaplaceTransform[{eqn, bc}, t, s] /. Rule @@@ ic


Now we have an ODE, solve it with DSolve:

tsol = u[t, x] /. First@DSolve[teqn/.
HoldPattern@LaplaceTransform[a_, __] :> a, u[t, x], x] // Simplify


Notice the replacement HoldPattern@LaplaceTransform[a_, __] :> a is necessary because DSolve has trouble in handling expression containing LaplaceTransform. The last step is to transform the solution back, but sadly InverseLaplaceTransform can't handle tsol. At this point, one work-around is to turn to numeric inverse Laplace transform, you can use this or this package for this task. But for your specific problem, we can circumvent the issue by expanding tsol with Fourier sine series:

easyFourierSinCoefficient[expr_, t_, {a_, b_}, n_] :=
FourierSinCoefficient[expr /. t -> t + a, t, n,
FourierParameters -> {1, Pi/(b - a)}] /. t -> t - a

easyTerm[t_, {a_, b_}, n_] := Sin[Pi/(b - a) n (t - a)]

term = easyTerm[x, {0, 1}, n];

coe = easyFourierSinCoefficient[tsol, x, {0, 1}, n]


$$-\left(i\left(\frac{(1+i) (-1)^n e^{i \sqrt{2} \sqrt{s}}}{(1+i) \pi n+i \sqrt{2} \sqrt{s}}\right.\right....$$

coe still looks complex, but inspired by those (-1)^ns in it, we split it to odd and even part and simplify:

oddcoe =
Simplify[coe /. n -> 2 m - 1, m > 0 && m ∈ Integers] /. m -> (1 + n)/2
(* (8 s)/(n^3 π^3 (n^4 π^4 + s^2)) *)

evencoe = Simplify[coe /. n -> 2 m, m ∈ Integers] /. m -> n/2
(* 0 *)


InverseLaplaceTransform can handle the series form of the transformed solution without difficulty:

soloddterm = Function[n, #] &@InverseLaplaceTransform[oddcoe term, s, t]
(* Function[n, (8 Cos[n^2 π^2 t] Sin[n π x])/(n^3 π^3)] *)


To find the final solution, just summate:

solgenerator[n_] := Compile[{t, x}, #] &@Total@soloddterm@Range[1, n, 2];

sol = solgenerator[200];

Animate[Plot[sol[t, x], {x, 0, 1}, PlotRange -> .3], {t, 0, 1}]


The animation is similar to the one in the subsequent solution so I'd like to omit it.

# Fully NDSolve-based Numeric Solution

Go back to the old-fashioned "TensorProductGrid", set "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 100} (or NDSolve will set "ScaleFactor" to 0 so the inconsistent b.c.s will be severely ignored, for more information, check the obscure tutorial) and DifferenceOrder -> 2:

mol[n_Integer, o_:"Pseudospectral"] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,
"MinPoints" -> n, "DifferenceOrder" -> o}}
mol[tf:False|True,sf_:Automatic]:={"MethodOfLines",
"DifferentiateBoundaryConditions"->{tf,"ScaleFactor"->sf}}

tu = NDSolveValue[{eqn, ic, bc}, u, {t, 0, 10}, {x, 0, 1},
Method -> Union[mol[27, 2], mol[True, 100]], MaxSteps -> Infinity];

Animate[Plot[tu[t, x], {x, 0, 1}, PlotRange -> .3], {t, 0, 10}]


NDSolve spits out the NDSolveValue::eerr warning, but in many cases NDSolveValue::eerr isn't a big deal, and the result indeed looks OK:

# Partly NDSolve-based Numeric Solution

Theoretically we can also set "DifferentiateBoundaryConditions" -> False to avoid the inconsistent b.c.s being ignored, but strangely NDSolve spits out the icfail warning and fails. I'm not sure about reason, but found that we can manually discretize the spatial derivative and solve the obtained ODE set with NDSolve to avoid the issue.

First, let's define a function pdetoode that discretizes PDEs to ODEs (Additionally, though not related to OP's problem, I've also define a function pdetoae that discretizes differential equations to algebraic equations based on pdetoode. A rebuild function is also created to combine the list of InterpolatingFunctions to a single InterpolatingFunction):

Clear[fdd, pdetoode, tooderule, pdetoae, rebuild]
fdd[{}, grid_, value_, order_, periodic_] := value;
fdd[a__] := NDSolveFiniteDifferenceDerivative@a;

pdetoode[funcvalue_List, rest__] :=
pdetoode[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[[1]], rest];
pdetoode[{func__}[var__], rest__] := pdetoode[Alternatives[func][var], rest];
pdetoode[front__, grid_?VectorQ, o_Integer, periodic_: False] :=
pdetoode[front, {grid}, o, periodic];

pdetoode[func_[var__], time_, {grid : {__} ..}, o_Integer,
periodic : True | False | {(True | False) ..} : False] :=
With[{pos = Position[{var}, time][[1, 1]]},
With[{bound = #[[{1, -1}]] & /@ {grid}, pat = Repeated[_, {pos - 1}],
spacevar = Alternatives @@ Delete[{var}, pos]},
With[{coordtoindex =
Function[coord,
{-1, PossibleZeroQ[# - #2[[-1]]]}}, All] &, {coord, bound}]]},
tooderule@Flatten@{
((u : func) | Derivative[dx1 : pat, dt_, dx2___][(u : func)])[x1 : pat, t_,
x2___] :> (Sow@coordtoindex@{x1, x2};
fdd[{dx1, dx2}, {grid}, Outer[Derivative[dt][u@##]@t &, grid],
"DifferenceOrder" -> o, PeriodicInterpolation -> periodic]),
inde : spacevar :>
With[{i = Position[spacevar, inde][[1, 1]]}, Outer[Slot@i &, grid]]}]]];

tooderule[rule_][pde_List] := tooderule[rule] /@ pde;
tooderule[rule_]@Equal[a_, b_] :=
Equal[tooderule[rule][a - b], 0] //. eqn : HoldPattern@Equal[_, _] :> Thread@eqn;
tooderule[rule_][expr_] := #[[Sequence @@ #2[[1, 1]]]] & @@ Reap[expr /. rule]

pdetoae[funcvalue_List, rest__] :=
pdetoae[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[[1]], rest];
pdetoae[{func__}[var__], rest__] := pdetoae[Alternatives[func][var], rest];

pdetoae[func_[var__], rest__] :=
Module[{t},
Function[pde, #[
pde /. {Derivative[d__][u : func][inde__] :>
Derivative[d, 0][u][inde, t], (u : func)[inde__] :> u[inde, t]}] /. (u : func)[
i__][t] :> u[i]] &@pdetoode[func[var, t], t, rest]]

rebuild[funcarray_, grid_?VectorQ, timeposition_: 1] :=
rebuild[funcarray, {grid}, timeposition]

rebuild[funcarray_, grid_, timeposition_?Negative] :=
rebuild[funcarray, grid, Range[Length@grid + 1][[timeposition]]]

rebuild[funcarray_, grid_, timeposition_: 1] /; Dimensions@funcarray === Length /@ grid :=
With[{depth = Length@grid},
ListInterpolation[
Transpose[Map[DeveloperToPackedArray@#["ValuesOnGrid"] &, #, {depth}],
Append[Delete[Range[depth + 1], timeposition], timeposition]],
Insert[grid, Flatten[#][[1]]["Coordinates"][[1]], timeposition]] &@funcarray]


The syntax of pdetoode is as follows: 1st argument is the function to be discretized (which can be a list i.e. pdetoode can handle PDE system), 2nd argument is the independent variable in the resulting ODE system (usually it's the variable playing the role of "time" in the underlying model), 3rd argument is the list of spatial grid, 4th argument is difference order, 5th argument is to determine whether periodic b.c. should be set or not. (5th argument is optional, the default setting is False. )

Notice pdetoode is a general purpose function. You may feel some part of the source code confusing. To understand it, just notice the following truth:

1. a /. a | b[m_] :> {m} outputs {}.
2. Derivative[][u] outputs u.

Then discretize eqn, ic and bc and remove redundant equations:

lb = 0; rb = 1;

torder = 2;
(* Difference order of x: *)
xdifforder = 2;

points = 25;
grid = Array[# &, points, {lb, rb}];

(* There're 4 b.c.s, so we need to remove 4 equations from every PDE/i.c.,
usually the difference equations that are the "closest" ones to the b.c.s
are to be removed: *)
removeredundant = #[[3 ;; -3]] &;
(* Use pdetoode to generate a "function" that discretizes the spatial derivatives of
PDE(s) and corresponding i.c.(s) and b.c.(s): *)
ptoofunc = pdetoode[u[t, x], t, grid, xdifforder];

odeqn = eqn // ptoofunc // removeredundant;
odeic = removeredundant/@ptoofunc@ic;
odebc = bc // ptoofunc;
(* Another possible treatment for i.c.s and b.c.s: *)
(* odeic = ptoofunc@ic; *)
(* odebc = With[{sf=100},
Map[D[#, {t,torder}]+sf#&,bc//ptoofunc,{2}]];*)

sollst = NDSolveValue[{odebc, odeic, odeqn}, u /@ grid, {t, 0, 10}, MaxSteps -> Infinity];
(* Rebuild the solution for the PDE from the solution for the ODE set: *)
sol = rebuild[sollst, grid];

Animate[Plot[sol[t, x], {x, 0, 1}, PlotRange -> .3], {t, 0, 10}]


The animation is similar to the one in the aforementioned solution so I'd like to omit it. This approach seems to be more robust than the fully NDSolve-based one, because even if the xordereqn i.e. the difference order for spatial derivative is set to 4, it's still stable, while the fully NDSolve-based one becomes wild when t is large.

• Thank you for this. It works but I get different values to you with non physical displacements of about 20 maximum which is larger than the maximum starting displacement ( 0.25). So there is something suspect. I am using version 11. The tutorial is something to read in detail thanks.
– Hugh
Oct 5, 2016 at 17:19
• @hugh I tested the code on Wolfram Cloud and got the same result as in v9: i.stack.imgur.com/ppP43.png Oct 6, 2016 at 3:28
• @Hugh I just noticed that "ScaleFactor"->1 is still too small, see my update, now the solution should be reasonable. I've also added an analytic solution to my answer. Oct 6, 2016 at 6:43
• This is excellent work -thank you. The analytic solution is always tricky since there are two branch cuts to be avoided. This is a minimum working example. The problem is how to extend this to more complex problems. I hope you don't mind if I add a more complicated version...
– Hugh
Oct 6, 2016 at 8:30
• @xzczd Thanks for keeping this uptodate. It seems to be a problem of continued interest.
– Hugh
Mar 24, 2019 at 18:51

Direct Analytic Solution by Separation of Variables

Off[General::wrsym]
Clear["Global*"]

f[x_] = x (1 - x);

pde = D[u[t, x], {t, 2}] + D[u[t, x], {x, 4}] == 0;

ic = {u[0, x] == f[x], Derivative[1, 0][u][0, x] == 0};
bc = {u[t, 0] == 0, u[t, 1] == 0, Derivative[0, 2][u][t, 0] == 0,
Derivative[0, 2][u][t, 1] == 0};


Separate variables in the form

u[t_, x_] = T[t] X[x]

pde = pde/u[t, x] // Apart
(* D[T[t],t,t]/T[t]+D[X[x],{x,4}]/X[x]==0 *)


The first term depends on t, the other x. That can only happen if each term is equal to a constant. We want sinusoidal in t so that we set

teq = D[T[t], t, t]/T[t] == -w^2;

T[t_] = T[t] /. (DSolve[teq, T[t], t][[1]]) /. {C[1] -> c1, C[2] -> c2}
(* c1 Cos[t w]+c2 Sin[t w] *)

From ic[[2]], we can eliminate c2 right away

c2 = c2 /. Solve[ic[[2]], c2][[1]]
(* 0 *)


The x equation

 xsol = ((DSolve[pde, X[x], x] // Flatten) /. {C[1] -> c3, C[2] -> c4,
C[3] -> c5, C[4] -> c6}) // ExpToTrig // Simplify

(* {X[x]->c3 Cos[Sqrt[w] x]+(c6-c4) Sinh[Sqrt[w] x]+(c4+c6)
Cosh[Sqrt[w] x]+c5 Sin[Sqrt[w] x]} *)

X[x_] = X[x] /. % /. {c6 - c4 -> c4, c4 + c6 -> c6};

u[t,x]
(* c1 Cos[t w] (c3 Cos[Sqrt[w] x]+c4 Sinh[Sqrt[w] x]+c5 Sin[Sqrt[w]
x]+c6 Cosh[Sqrt[w] x]) *)

bc[[1]]
(* c1 (c3+c6) Cos[t w]==0 *)


From which

c6 = -c3;


And to consolidate constants

c1 = 1;


Now

bc[[3]]
(* -2 c3 w Cos[t w]==0 *)


From which

c3 = 0;

bc[[2]]
(* Cos[t w] (c4 Sinh[Sqrt[w]]+c5 Sin[Sqrt[w]])==0 *)

c4 = c4 /. Solve[bc[[2]], c4][[1]]
(* -c5 Sin[Sqrt[w]] Csch[Sqrt[w]] *)

bc[[4]]
(* -2 c5 w Sin[Sqrt[w]] Cos[t w]==0 *)


Instead of solving for the trivial c5 = 0 solution, we will solve for w.

Reduce[ {Sin[Sqrt[w]] == 0, w > 0}, w]

(*C[1]\[Element]Integers&&((C[1]>=1&&w\[Equal]4 \[Pi]^2 \
C[1]^2)||(C[1]>=0&&w==4 \[Pi]^2 C[1]^2+4 \[Pi]^2 C[1]+\[Pi]^2))*)


All boils down to

w = n^2 Pi^2;

\$Assumptions = n \[Element] Integers && n > 0;

u[t_, x_] = u[t, x] // Simplify
(* c5 Cos[Pi^2 n^2 t] Sin[PI n x] *)

(ic[[1, 1]] // Simplify) == ic[[1, 2]]
(* c5 Sin[Pi n x]==(1-x) x *)


We use orthogonality to solve for c5. Multiply each side of the above by Sin[n Pi x] and Integrate over the length of the beam

    Integrate[%[[1]] Sin[n Pi x], {x, 0, 1}] ==
Integrate[%[[2]] Sin[n Pi x], {x, 0, 1}] // Simplify
(* Pi^3 c5 n^3+4 (-1)^n==4 *)

c5 = c5 /. Solve[%, c5][[1]] // Simplify
(* -((4 ((-1)^n-1))/(Pi^3 n^3)) *)

u[t, x]
(* -((4 ((-1)^n-1) Cos[Pi^2 n^2 t] Sin[Pi n x])/(Pi^3 n^3))*)


We can see that (-1)^n-1 will make all even n terms equal to 0. Rather than just choosing odd terms in the series for u, we can change n to 2m-1, which will give us only the odd n terms.

    um[t_, x_] = ((u[t, x] /. n -> 2 m - 1) //
Simplify[#, m > 0 && m \[Element] Integers] &)
(* (8 Cos[Pi^2 (1-2 m)^2 t] Sin[Pi (2 m-1) x])/(Pi^3 (2 m-1)^3) *)


The analytic solution is

 u[t_, x_] :=
8/Pi^3 Sum[(
Cos[Pi^2 (1 - 2 m)^2 t] Sin[Pi (2 m - 1) x])/(2 m - 1)^3, {m,  1, \[Infinity]}]


which is not practical for computing. Similar to the method of xzczd above:

term = Function[m, #] &@um[t, x]

(* Function[m,(8 Cos[Pi^2 (1-2 m)^2 t] Sin[Pi (2 m-1) x])/(Pi^3 (2 m-1)^3)]*)

mterms[m_] := Compile[{t, x}, #] &@Total@term@Range[1, m]

U = mterms[100];

Animate[Plot[U[t, x], {x, 0, 1}, PlotRange -> .3], {t, 0, 1}]


We get the same plot as the Laplace Transform solution.

• Thank you. This needs to be here as the classic text book solution of the equation. I agree that this is simpler than going into Laplace transforms. Unfortunately it only works for uniform beams. I don't think we can use this for beams where the mass or stiffness varies along the length. This is why I was hoping to use the finite element method. Very useful.
– Hugh
Nov 13, 2017 at 9:35

Update:

I have answered a similar question here.

Here are two (partial) ideas:

One could try to make use of the TensorProductGrid as a discretization method.

ClearAll[f];
f[x_] := x (1 - x)

tu = NDSolveValue[{D[u[t, x], {t, 2}] + D[u[t, x], {x, 4}] == 0,
u[0, x] == f[x], Derivative[1, 0][u][0, x] == 0,
u[t, 0] == 0, u[t, 1] == 0
(*,
Derivative[0,2][u][t,0]\[Equal]0,
Derivative[0,2][u][t,1]\[Equal]0
*)
}, u, {t, 0, 1}, {x, 0, 1},
Method -> {"PDEDiscretization" -> "MethodOfLines"}];


DirichletCondition will trigger a FEM attempt, which does not work as the FEM can not handle 4th order spatial derivatives (V11). Note that I disabled the derivatives as that gave inconsistent initial and boundary conditions. Perhaps you know what needs to be done.

The second idea is to treat this as a pure spatial problem.

ClearAll[f];
f[x_] := x (1 - x)

tu = NDSolveValue[{
D[u[t, x], {t, 2}] + D[m[t, x], {x, 2}] == 0,
D[u[t, x], {x, 2}] == m[t, x],
DirichletCondition[u[t, x] == f[x], t == 0],
DirichletCondition[u[t, x] == 0, x == 1 || x == 0],
DirichletCondition[m[t, x] == 0, True]}, {u, m}, {t, 0, 1}, {x, 0,
1}, Method -> {"PDEDiscretization" -> {"FiniteElement"}}];


The issue with your decoupling of the equations is that the second equation is no longer time dependent. So I was thinking of making this a purely spatial problem. Have a look and see if the solutions are any good. Maybe the DirichletCondition on m does not need to be True but something more specific. I did not check. Hope this gives you a starting point.

• Sadly both attempts don't work... NDSolve probably ignores the inconsistent b.c. in the first case, and ignores the inconsistent i.c. in the second case… Oct 5, 2016 at 14:51
• @xzczd See the update for the second idea. Oct 5, 2016 at 15:35
• The second idea works thank you. However, I don't understand the initial conditions and boundary conditions. The first DirichletCondition is the initial displacement, the second is a boundary condition for each end, the third is unclear as you suggest. However, there is no initial condition for velocity. Is this assumed to be zero?
– Hugh
Oct 5, 2016 at 17:05
• I have continued to look at your second solution. Using DirichletCondition[m[t, x] == 0, x == 1 || x == 0] instead of the last boundary condition gives a different solution. Also, how would I implement an initial condition with zero displacement and an initial velocity? Thanks
– Hugh
Oct 5, 2016 at 17:59
• I think it's OK to keep it undeleted, you've claimed it's "a starting point" in the body of answer anyway. Look forward to your complete answer :) Oct 6, 2016 at 12:42

The analytical solution of E-B beam after assuming harmonic dependency,

    L = 1;
sol = Flatten[DSolve[(D[y[x], {x, 4}] - b^4 y[x]) == 0, y[x], x]];
a = y[x] /. sol;
beamsol = Simplify[ExpToTrig[a]];
(*below is the text book form general sol*)
generalsol =
Flatten[beamsol /. {C[1] -> C1, (C[2] + C[4]) -> C2,
C[3] -> C3, (C[4] - C[2]) -> C4}];
(*BC*)
e[1] = beamsol /. x -> 0;
e[2] = D[beamsol, {x, 2}] /. x -> 0;
e[3] = beamsol /. x -> L;
e[4] = D[beamsol, {x, 2}] /. x -> L;
eq = Table[e[i], {i, 1, 4}];
var = Table[C[i], {i, 1, 4}];
R = Normal@CoefficientArrays[eq, var][[2]];
MatrixForm[R];
P = Det[R];
s1 = NSolve[P == 0 && 0 < b < 10];
s2 = b /. s1;
NN = Flatten[NullSpace[R /. b -> s2[[1]]]];
beamsol = (beamsol /.
Table[var[[i]] -> NN[[i]], {i, 1, Length[NN]}]) /. b -> s2[[1]];
Plot[beamsol, {x, 0, L}]
`
• This isn't what OP wants, is it? Feb 3, 2019 at 8:10
• Thank you this is the standard solution for harmonic motion. However, I want a solution from initial conditions. This could be built-up from a Fourier series of your solutions. However, I really wish to be general and be able to solve non-uniform beams for general initial conditions.
– Hugh
Feb 4, 2019 at 14:40