# Dynamic, nonlinear, damped Euler–Bernoulli beam equation

I would like to solve the 3 coupled PDEs describing a damped, nonlinear (i.e displacements in the $$x$$ direction along the beam need to be considered along with the $$y$$ displacements normally considered), inextensible, euler-bernoulli cantilever beam's motion. The equations are as follows:

1: where $$u(t,s)$$ is the displacement in the $$x$$ direction $$\mu\frac{\partial^2 u}{\partial t^2}+k\frac{\partial u}{\partial t}-\frac{\partial}{\partial s}\left(\lambda\left(1+\frac{\partial u}{\partial s}\right)\right)=0$$ 2: where $$v(t,s)$$ is the displacement in the $$y$$ direction

EDIT: To make things simpler, we expand the elasticity term to first order only. $$\mu\frac{\partial^2 v}{\partial t^2}+k\frac{\partial v}{\partial t}-\frac{\partial}{\partial s}\left(\lambda\frac{\partial v}{\partial s}\right)+EI\frac{\partial^4 v}{\partial s^4}=0$$ 3: inextensibility of beam condition $$\left(1+\frac{\partial u}{\partial s}\right)^2+\left(\frac{\partial v}{\partial s}\right)^2-1=0$$ $$\mu$$, $$EI$$ and $$k$$ are constants representing mass density, stiffness and drag coefficient, while EDIT: $$\lambda(t,s)$$ is the lagrange multiplier needed to maintain the inextensibility condition.

My boundary conditions are as follows (for a beam of unit length):

EDIT: New b.c.s. that include $$\lambda(t,s)$$

2 for $$u$$

$$u(t,0)=0$$

$$\lambda(t,1)\left(1+u^{(0,1)}(t,1)\right)=0$$

4 for $$v$$

$$v(t,0)=0$$

$$v^{(0,1)}(t,0)=0$$

$$v^{(0,2)}(t,1)=0$$

$$EI\left(v^{(0,3)}(t,1)\right)-\lambda(t,1)v^{(0,1)}(t,1)=1$$

EDIT: I use a point load of 1N at the end of the beam to perturb it instead so the beam is completely stationary in the i.c.s

Next, my initial conditions are of a fully stationary, straight beam:

$$u(0,x)=0$$

$$u^{(1,0)}(0,x)=0$$

$$v(0,x)=0$$

$$v^{(1,0)}(0,x)=0$$

My approach is to use the pdetoode function by @xzczd. I had no problems following and implementing the solution for the linearized Euler-Bernoulli equation as in this earlier question (that is, only small deflections in the $$y$$ direction are considered and inextensibility is ignored), but I have issues with generating the correct number of ODEs for NDSolve to solve for the nonlinear system.

Here is my implementation within Mathematica:

EDIT: Here is the code with the updated functions and b.c.s.

μ = 1;
EI = 10;
k = 1;
eqn1 = μ*D[u[t, s], {t, 2}] + k*D[u[t, s], t] -
D[λ[t, s]*(1 + D[u[t, s], s]), s] == 0;
eqn2 = μ*D[v[t, s], {t, 2}] + k*D[v[t, s], t] -
D[λ[t, s]*D[v[t, s], s], s] + EI*D[v[t, s], {s, 4}] == 0;
eqn3 = (1 + D[u[t, s], s])^2 + D[v[t, s], s]^2 - 1 == 0;
bc1 = {u[t, 0] == 0, λ[t, 1]*(1 + Derivative[0, 1][u][t, 1]) ==
0};
bc2 = {v[t, 0] == 0, Derivative[0, 1][v][t, 0] == 0,
Derivative[0, 2][v][t, 1] == 0,
EI*Derivative[0, 3][v][t, 1] - λ[t, 1]*
Derivative[0, 1][v][t, 1] == 1};
ic1 = {u[0, s] == 0, Derivative[1, 0][u][0, s] == 0};
ic2 = {v[0, s] == 0, Derivative[1, 0][v][0, s] == 0};


Generate the finite-difference grid (i use 9 points for now to speed things up but more would be needed for accuracy):

lb = 0;
rb = 1;
torder = 2;
sdifforder = 2;

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


Generate ODEs using pdetoode:

removeredundant1 = #[[2 ;; -2]] &;
removeredundant2 = #[[3 ;; -3]] &;

ptoofunc = pdetoode[{u, v, λ}[t, s], t, grid, sdifforder];

odeqn1 = eqn1 // ptoofunc // removeredundant1;
odeqn2 = eqn2 // ptoofunc // removeredundant2;
odeqn3 = eqn3 // ptoofunc;
odeic1 = Flatten[removeredundant1 /@ ptoofunc@ic1];
odeic2 = Flatten[removeredundant2 /@ ptoofunc@ic2];
odebc1 = bc1 // ptoofunc;
odebc2 = bc2 // ptoofunc;


EDIT: For the ODEs generated from eqn 1 and i.c. 1, I remove 2 equations to make space for 2 b.c.s. For the ODEs generated from eqn 2 and i.c. 2, I remove 4 equations to make space for 4 b.c.s. I remove no equations from the ODEs generated from eqn 3. The system is now correctly determined and I no longer get the error NDSolve::overdet:

tEnd = 1;
sollst = NDSolveValue[
Join[odebc1, odebc2, odeic1, odeic2, odeqn1, odeqn2, odeqn3],
Join[u /@ grid, v /@ grid, λ /@ grid], {t, 0, tEnd},
MaxSteps -> Infinity];


Now, the issue is as follows:

NDSolveValue::icfail: Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions.

• How was this system obtained? – Alex Trounev Dec 6 '18 at 23:31
• As far as I understand, the 3rd unknown function to be calculated is the lagrange multiplier $\lambda(t)$. The first two equations are derived from the Euler-Lagrange equation with the lagrange multiplier being the force required to maintain inextensibility. Drag is accounted for with the extended hamiltonian principle. The last term in eqn 2 is obtained by expanding the elasticity term in the lagrangian to 4th order. – Unevaluated-Sequence Dec 7 '18 at 5:42
• It is only a comment, maybe an irrelevant one. There is a statement of the problem for such a rod, that is more simple and yields equations of the second, rather than of the fourth order. It works like this due to a nonlinear transform done from the very beginning. If you are interested, have a look at L. D. Landau and E. M. Lifshitz, Theory of Elasticity, 3 ed. (Pergamon Press, Oxford, 1986), Problems 1 to 4 for § 19. These problems are for the static case, but the generalization to the dynamic one seems rather straightforward. – Alexei Boulbitch Dec 7 '18 at 9:04
• @xzczd thanks, I'm looking into it. That said, time-domain discretization is rather unideal for my eventual application goal because I will be needing to trigger piecewise conditions and possibly even use WhenEvent (these methods are already working for me for the lineaerized equtions). – Unevaluated-Sequence Dec 7 '18 at 9:28
• At this juncture, I am still stuck with the same error. I am more inclined to think this problem is related to the solving method rather than an incorrect formulation of equations, because I found a reference that formulates similar equations to myself for a similar system: Yabuno, Hiroshi, et al. "Suppression of parametric resonance in cantilever beam with a pendulum (Effect of static friction at the supporting point of the pendulum)." Journal of vibration and acoustics 126.1 (2004): 149-162. – Unevaluated-Sequence Dec 7 '18 at 16:02

For those values u, v that arise in this problem, we can use eqn3 for finding u. The equation for v can be solved by the usual "MethodOfLines", for example

f[x_] := x;
μ = 1;
EI = 10;
k = 1;
λ[t_] := 1
eqn1 = μ*D[u[t, s], {t, 2}] + k*D[u[t, s], t] -
D[λ[t]*(1 + D[u[t, s], s]), s] == 0;
eqn2 = μ*D[v[t, s], {t, 2}] + k*D[v[t, s], t] -
D[λ[t]*D[v[t, s], s], s] +
EI*(D[(1 + D[v[t, s], s]^2 + D[v[t, s], s]^4)*
D[v[t, s], {s, 2}], {s, 2}] +
D[(1 + (1/2)*D[v[t, s], s]^2)*D[v[t, s], {s, 2}]^2, s]) == 0;
eqn3 = (1 + D[u[t, s], s])^2 + D[v[t, s], s]^2 - 1 == 0;
ic1 = {u[0, s] == 0};
ic2 = {v[0, s] == 0, Derivative[1, 0][v][0, s] == f[s]};
bc1 = {u[t, 0] == 0, Derivative[0, 1][u][t, 0] == 0};
bc2 = {v[t, 0] == 0, Derivative[0, 1][v][t, 0] == f[0],
Derivative[0, 3][v][t, 1] == 0, Derivative[0, 2][v][t, 1] == 1};
tEnd = 1;

sol1 = NDSolveValue[{eqn2, ic2, bc2}, v, {t, 0, tEnd}, {s, 0, 1},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 25, "MaxPoints" ->25,
"DifferenceOrder" -> 2}}, MaxSteps -> 10^6];


Use pdetoode to solve this problem and compare two solutions.

lb = 0; rb = 1;

torder = 2;

xdifforder = 2;

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

removeredundant = #[[3 ;; -3]] &;

ptoofunc = pdetoode[v[t, s], t, grid, xdifforder];

odeqn = eqn2 // ptoofunc // removeredundant;
odeic = removeredundant /@ ptoofunc@ic2;
odebc = bc2 // ptoofunc;

sollst = NDSolveValue[{odebc, odeic, odeqn}, v /@ grid, {t, 0, tEnd},
MaxSteps -> Infinity];

sol = ListInterpolation[
DeveloperToPackedArray@#["ValuesOnGrid"] & /@ sollst //
Transpose, {Flatten@sollst[[1]]["Grid"], grid}];

{ContourPlot[sol1[t, s], {t, 0, tEnd}, {s, 0, 1}, Contours -> 20,
ColorFunction -> Hue, FrameLabel -> {"t", "s"},
PlotLegends -> Automatic, PlotLabel -> "Method Of Lines"],
ContourPlot[sol[t, x], {t, 0, tEnd}, {x, 0, 1}, Contours -> 20,
ColorFunction -> Hue, FrameLabel -> {"t", "s"},
PlotLegends -> Automatic, PlotRange -> All,
PlotLabel -> "pdetoode"]}


Solutions do not match. Which method is wrong? For "MethodOfLines", we have two messages.

NDSolveValue::ibcinc: Warning: boundary and initial conditions are inconsistent.

NDSolveValue::eerr: Warning: scaled local spatial error estimate of 1642.781029898344 at t = 1. in the direction of independent variable s is much greater than the prescribed error tolerance. Grid spacing with 25 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or a smaller grid spacing can be specified using the MaxStepSize or MinPoints method options.


For pdetoode there is no messages. As advised by xzczd, we will add options to

sol1 = NDSolveValue[{eqn2, ic2, bc2}, v, {t, 0, tEnd}, {s, 0, 1},
Method -> {"MethodOfLines",
"DifferentiateBoundaryConditions" -> {True,
"ScaleFactor" -> 100},
"SpatialDiscretization" -> {"TensorProductGrid",
"MaxPoints" -> 25, "MinPoints" -> 25, "DifferenceOrder" -> 2}},
MaxSteps -> 10^6];


then the results become similar, although they do not match

Now we have to solve the original problem.The system of equations of the third order in powers of amplitude is given by Pramod Malatkar. Nonlinear Vibrations of Cantilever Beams and Plates - see eqs (2.54), (2.55), (2, 56), (2.57). We put $$w = 0, Q_u = 0, D_\zeta = EI, Q_v = f (t)$$ and use the method of the false transient and method of lines with pdetoode

Clear[f];
f[t_] := Sin[2*Pi*t];
μ = 1;
EI = 10;
k = 1;
c = 0;
dif = 1/100; tm = 8;
eqn1 = D[F[t, s], {s, 2}] - dif*D[F[t, s], t] -
D[D[v[t, s], s]^2, {t, 2}] == 0;
eqn2 = μ*D[v[t, s], {t, 2}] + k*D[v[t, s], t] +
c*D[v[t, s], t]*Abs[D[v[t, x], t]] + EI*D[v[t, s], {s, 4}] +
EI*(D[D[v[t, s], s]*D[D[v[t, s], s]*D[v[t, s], {s, 2}], s],
s]) + μ/2*D[D[v[t, s], s]*F[t, s], s] - f[t] == 0;

ic1 = {F[0, s] == 0};
ic2 = {v[0, s] == 0, Derivative[1, 0][v][0, s] == 0};
bc1 = {F[t, 0] == 0, F[t, 1] == 0};
bc2 = {v[t, 0] == 0, Derivative[0, 1][v][t, 0] == 0,
Derivative[0, 3][v][t, 1] == 0, Derivative[0, 2][v][t, 1] == 0};
lb = 0; rb = 1;

torder = 4;

xdifforder = 2;

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

removeredundant1 = #[[2 ;; -2]] &;
removeredundant2 = #[[3 ;; -3]] &;

ptoofunc = pdetoode[{F[t, s], v[t, s]}, t, grid, xdifforder];
odeqn1 = eqn1 // ptoofunc // removeredundant1;
odeic1 = removeredundant1 /@ ptoofunc@ic1;
odebc1 = bc1 // ptoofunc;
odeqn2 = eqn2 // ptoofunc // removeredundant2;
odeic2 = removeredundant2 /@ ptoofunc@ic2;
odebc2 = bc2 // ptoofunc;

sollst = NDSolveValue[{odebc1, odeic1, odeqn1, odebc2, odeic2,
odeqn2}, v /@ grid, {t, 0, tm}];

sol = ListInterpolation[
DeveloperToPackedArray@#["ValuesOnGrid"] & /@ sollst //
Transpose, {Flatten@sollst[[1]]["Grid"], grid}];

sollst1 =
NDSolveValue[{odebc1, odeic1, odeqn1, odebc2, odeic2, odeqn2},
F /@ grid, {t, 0, tm}]; sol1 =
ListInterpolation[
DeveloperToPackedArray@#["ValuesOnGrid"] & /@ sollst1 //
Transpose, {Flatten@sollst1[[1]]["Grid"], grid}];

{ContourPlot[sol1[t, x], {t, 0, tm}, {x, 0, 1}, Contours -> 20,
ColorFunction -> Hue, FrameLabel -> {"t", "s"},
PlotLegends -> Automatic, PlotRange -> All, PlotLabel -> "F"],

ContourPlot[sol[t, x], {t, 0, tm}, {x, 0, 1}, Contours -> 20,
ColorFunction -> Hue, FrameLabel -> {"t", "s"},
PlotLegends -> Automatic, PlotRange -> All, PlotLabel -> "v"],
Plot[sol1[t, .5], {t, 0, tm}, PlotRange -> All, PlotLabel -> "F"],
Plot[sol[t, 1], {t, 0, tm}, PlotRange -> All, PlotLabel -> "v"]}


here F[] is similar to $$\lambda$$

• How is eqn 1 accounted for? Also as per my edits in the OP, $\lambda$ is a function of both $s$ and $t$ and is the 3rd unknown parameter to be solved for. The system is now correctly determined. – Unevaluated-Sequence Dec 7 '18 at 8:43
• Yes, but this is a completely different problem. Who will be in time for your thought? – Alex Trounev Dec 7 '18 at 12:46
• Well, to be honest, I think you've answered the question too quickly. The question is indeed somewhat unclear and involves mistakes at beginning, but that's not a reason for posting hasty answer. Leaving a comment asking for clarification is the better way to go, I think. – xzczd Dec 7 '18 at 13:17
• @xzczd Note that I answered another question that the author deleted. Should I delete my answer or put a question about deleting a post? – Alex Trounev Dec 7 '18 at 14:50
• Just see your update. The result given by pdetoode is more reliable, if you add Method -> {"MethodOfLines", "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 100}, "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 25, "MinPoints" -> 25, "DifferenceOrder" -> 2}} to sol1`, the result will be similar. The underlying issue is discussed in this post. – xzczd Dec 9 '18 at 4:50