# Two-directional vibration of Euler–Bernoulli beam with Lagrange multiplier

## Background

I'm going to investigate a beam-pendulum coupling system (in this question I won't consider the pendulum though), that is, a spherical pendulum is suspended on the tip of a cantilever beam. therefore it may be essential to write the vibrational displacement in two functions, which are constrained by inextensibility of the beam.

I'm following this article, which starts from a Lagrangian with a Lagrange multiplier.

https://doi.org/10.1119/10.0001389

I upload 2 pages of this paper here.

I've set $$T_1,T_2=0$$ to ignore effects from pendulum etc., i.e. there is the beam only.

## PDE

(* PDE *)
eqn = {
\[Mu] Derivative[2,0][u][t,s] + c Derivative[1,0][u][t,s] - D[\[Lambda][t, s](1+Derivative[0,1][u][t,s]), s] == 0,
\[Mu] Derivative[2,0][v][t,s] + c Derivative[1,0][v][t,s] - D[\[Lambda][t, s]Derivative[0,1][v][t,s], s] +
EI(D[(1+Derivative[0,1][v][t,s]^2+Derivative[0,1][v][t,s]^4)Derivative[0,2][v][t,s], {s, 2}] +
D[(1+Derivative[0,1][v][t,s]^2/2)Derivative[0,2][v][t,s]^2, s]) == 0
}

(* Boundary Condition *)
bcrhs = D[(1+D[v[t,s], s]^2+D[v[t,s], s]^4)D[v[t,s], s, s], s] + (1+D[v[t,s], s]^2/2)D[v[t,s], s, s]^2 /. s->l
bc = {
u[t, 0] == 0,
Derivative[0, 1][u][t, l] == -1,
v[t, 0] == 0,
Derivative[0, 1][v][t, 0] == 0,
Derivative[0, 2][v][t, l] == 0,
\[Lambda][t, l]Derivative[0, 1][v][t, l] == EI bcrhs
};

(* Initial Condition, a trivial one though *)
ic = {
u[0, s] == 0,
Derivative[1, 0][u][0, s] == 0,
v[0, s] == 0,
Derivative[1, 0][v][0, s] == 0
}

(* Constraint *)
constraint = (1+Derivative[0,1][u][t,s])^2+Derivative[0,1][v][t,s]^2==1

(* Parameters *)
l = 1;
EI = 1;
\[Mu] = 1;
c = 1;


## My attempts

I've tried manual discretization of PDEs and BCs for spatial variable $$s$$(with the help of pdetoode from 127997). Redundant ODEs are removed (2 of u, 4 of v). Finally I got a correct amount of DEs, which is $$2m$$ where $$m$$ is the grid size.

For the constraint, as its sign is not fixed, it cannot be directly solved. Thus I followed the instruction of the wiki, differentiated the constraint for 2 times and substituted high-order derivatives by terms including $$\lambda$$ to obtain a constraint with $$\lambda$$. Then differentiate it once more to obtain a DE version constraint.

For the initial condition for the DE version constraint, I don't know what to do. After substitute ICs into the constraint equation with $$\lambda$$, I got a equation which evaluates to True. Thus, being inspired by some observation of arbitrariness, I tried to set constraint IC $$\lambda(0,s)=0$$. Then it's discretized to $$m$$ ICs.

Finally, I got $$2m$$ second order ODEs, $$m$$ first order ODEs, and $$5m$$ initial conditions. When I throw them into NDSolve(without set any option), it says:

1. NDSolve::parpiv: Zero pivot was detected during the numerical factorization or there was a problem in the iterative refinement process. It is possible that the matrix is ill-conditioned or singular.
2. NDSolve::ivres: NDSolve has computed initial values that give a zero residual for the differential-algebraic system, but some components are different from those specified. If you need them to be satisfied, giving initial conditions for all dependent variables and their derivatives is recommended.
3. NDSolve::index: The DAE solver failed at t = 0. The solver is intended for index 1 DAE systems and structural analysis indicates that the DAE index is 3. The option Method->{"IndexReduction"->Automatic} may be used to reduce the index of the system.

Set Method -> {"IndexReduction" -> True} gives NDSolve::ivres and NDSolve::nderr.

## My question

I want to know what's the correct way to solve this PDE system. (Use built-in discretizer is also OK, but I'm afraid that it may be not so handy for introduce the pendulum coupling)

If you think that my code for attempts is a needed info, please comment.

• The paper is behind a pay wall so I can't read it. The Bernoulli Euler beam has forth derivatives of the displacement for the PDE. I can't immediately see those. Is this a beam that has its axis horizontal? Have you also got damping in there? Is your equation for the Lagrangian or the equation of motion?
– Hugh
Commented Jan 17 at 20:18
• @Hugh 4th derivative can be found in the equation for v (where there're two nested $\partial^2$). The beam axis is horizontal, and gravity is ignored and we only care for the horizontal dimensions. There's a linear damping. The equation is a Euler-Lagrange equation of a Lagrangian. Commented Jan 18 at 0:08
• @Hugh I've uploaded two main pages of that paper. Commented Jan 18 at 0:23
• There is a typo in constraint, it should be (1+Derivative[0,1][u][t,s])^2+Derivative[0,1][v][t,s]^2==1 - see equation (3) in the paper. Commented Jan 18 at 1:14
• @AlexTrounev Thanks. I'll edit, though fixing it does not change any complaint of NDSolve. Commented Jan 18 at 1:20

To solve this problem we can use FDM code based on explicit Verlet step for the pendulum bob motion and implicit algorithm for the cantilever beam deformation. We use standard simplification to derive parameters $$u,\lambda$$ from equations (3) and (8) in a form

$$u_s=\sqrt{1-v_s^2}-1, u=\int_0^su_sds$$

$$\lambda (1+u_s)=\lambda_0+\int_0^s(\mu u_{tt}+c u_t)ds$$

Note that code with $$\lambda(t,s)$$ in a form above is terribly slow, while numerical results are not differ from that computed with $$\lambda(t)$$ only. Therefore we use simplified code as follows

(*Parameters *)

l = .815; n = 15; h = l/(n - 1);
EI = .096; g = 9.8;
\[Mu] = 8.41 10^-3; mp = 0.033; R = 0.225; gamma = 2.02;
c = 1.51 10^-2; cp = 5.31 10^-4; alpha = 42/180 Pi;
dt = 5 10^-2; tmax = 20;
sgrid = Range[0, l, h];
tgrid = Range[0, tmax, dt]; jmax = nt = Length[tgrid]; nn =
Length[sgrid]; ds =
Table[NDSolveFiniteDifferenceDerivative[i, sgrid,
"DifferenceOrder" -> 4]["DifferentiationMatrix"], {i, 4}];

U = Array[uu, {nn, nt}]; V = Array[vv, {nn, nt}]; lambda =
Array[ll, {nn, nt}]; lambda0 = Array[ll0, {nt}]; xp =
Array[xx, {nt}]; yp = Array[yy, {nt}]; zp = Array[zz, {nt}]; T1 =
Array[tt1, {nt}]; T2 = Array[tt2, {nt}]; Vpx =
Array[vpx, {nt}]; Vpy =
Array[vpy, {nt}]; fu = -T2 Cos[alpha] + T1 xp/R;(*External force fu*)
fv = -T2 Cos[alpha] V[[-1, All]]/l + T1 yp/R;(*External force fv*)
us = Sqrt[1 - (ds[[1]] . V)^2] - 1;(*inextensibility constraint*)
Do[uu[i, j] =
If[i == 1, 0,
h/2 Sum[(us[[k - 1, j]] + us[[k, j]]), {k, 2, i}]];, {i, nn}, {j,
nt}];(*u[t,s]= Integrate[us,{x,0,s}] *)
lambda0 = fu/(1 + us[[-1, All]]);
(*the Lagrange multiplier \
\[Lambda](t,s)=\[Lambda]0(t)+Integrate[\[Mu] D[u[t,x],t,t]+c \
D[u[t,x],t],{x,s,l}]/(1+D[u[t,s],s])*)
bcrhs[j_] := (1 + 1/2 (ds[[1]] . V[[All, j]])^2) (ds[[2]] .
V[[All, j]])^2 +
ds[[2]] .
V[[All, j]] (2 ds[[1]] . V[[All, j]] ds[[2]] . V[[All, j]] +
4 (ds[[1]] . V[[All, j]])^3 ds[[2]] .
V[[All, j]]) + (1 + (ds[[1]] . V[[All, j]])^2 + (ds[[1]] .
V[[All, j]])^4) ds[[3]] . V[[All, j]];
(*Initial Condition*){x0, vx0, y0, vy0} = {.0, .1, 0,
0.01}; z0 = -Sqrt[R^2 - x0^2 - y0^2];
Do[uu[i, 1] = 0; uu[i, 2] = 0; vv[i, 1] = 0;
vv[i, 2] = 0;, {i, nn}]; Do[xx[i] = x0; vpx[i] = vx0; yy[i] = y0;
vpy[i] = vy0; zz[i] = z0; tt1[i] = -mp g R/z0;
tt2[i] = mp g/Sin[alpha], {i, 2}]; ll0[1] = -tt2[1] Cos[alpha];
ll0[2] = -tt2[2] Cos[alpha];
test = Table[.1, {nn}];

Do[(*Verlet explicit step*)
fx = If[j <= 2, 0,
mp (U[[-1, j]] + U[[-1, j - 2]] - 2 U[[-1, j - 1]])/dt^2 +
cp (U[[-1, j]] - U[[-1, j - 1]])/dt ];
fy = If[j <= 2, 0,
mp (V[[-1, j]] + V[[-1, j - 2]] - 2 V[[-1, j - 1]])/dt^2 +
cp (V[[-1, j]] - V[[-1, j - 1]])/dt ];
xx[j + 1] =
xx[j] + vpx[j] dt - (fx + cp vpx[j] + tt1[j] xx[j]/R) dt^2/2 ;
vpx[j + 1] = vpx[j] - (fx + cp vpx[j] + tt1[j] xx[j]/R) dt;
yy[j + 1] =
yy[j] + vpy[j] dt - (fy + cp vpy[j] + tt1[j] yy[j]/R) dt^2/2 ;
vpy[j + 1] = vpy[j] - (fy + cp vpy[j] + tt1[j] yy[j]/R) dt;
zz[j + 1] = -Sqrt[R^2 - xx[j + 1]^2 - yy[j + 1]^2];
tt1[j + 1] = -R/
zz[j + 1] (mp g + (mp (zz[j + 1] + zz[j - 1] - 2 zz[j])/dt^2 +
cp /dt/2 (zz[j + 1] - zz[j - 1])));
tt2[j + 1] = -tt1[j + 1] zz[j + 1]/R/ Sin[alpha];
(*beam implicit step*)
r2 = EI ds[[1]] . bcrhs[j + 1] -
lambda0[[j + 1]] (ds[[2]] . V[[All, j + 1]]);
eq2 = \[Mu] (V[[All, j + 1]] + V[[All, j - 1]] - 2 V[[All, j]]) +
c dt/2 (V[[All, j + 1]] - V[[All, j - 1]]) + dt^2 r2;
eqs = Table[Drop[Drop[eq2, 2], -2][[i]] == 0, {i, nn - 4}];
varAll = V[[All, j + 1]];
bC0 = {(ds[[1]] . V[[All, j + 1]])[[1]] == 0, V[[1, j + 1]] == 0};
bC1 = {lambda0[[j + 1]] (ds[[1]] . V[[All, j + 1]])[[-1]] -
EI bcrhs[j + 1][[-1]] ==
fv[[j + 1]], (ds[[2]] . V[[All, j + 1]])[[-1]] == 0};
eqAll = Join[eqs, bC0, bC1];
ini = Table[{varAll[[i]], test[[i]]}, {i, Length[varAll]}];
sol = Quiet@
FindRoot[eqAll, ini,
Method -> {"Newton", "StepControl" -> "TrustRegion"},
MaxIterations -> 20]; test = varAll /. sol;
Do[vv[i, j + 1] = test[[i]], {i, nn}];, {j, 2,
jmax - 1}] // AbsoluteTiming


It takes about 205s on my laptop. Visualization of the cantilever beam tip and the pendulum bob movement

ListLinePlot[Transpose[{tgrid, V[[-1, All]]}],
PlotStyle -> {Dashed, Red}]

ListLinePlot[{Transpose[{tgrid, xp}], Transpose[{tgrid, yp}],
Transpose[{tgrid, zp}]}, PlotLegends -> {"x", "y", "z"}]


Forces acting on the cantilever beam and pendulum

ListLinePlot[{Transpose[{tgrid, T1}], Transpose[{tgrid, T2}],
Transpose[{tgrid, fu}], Transpose[{tgrid, fv}]},
PlotLegends -> {"T1", "T2", "fu", "fv"}]


Lagrange multiplier $$\lambda$$

ListLinePlot[Evaluate[Transpose[{tgrid, lambda0}]],
PlotStyle -> {Dashed, Red}]


Animation

frame3D2 =
Table[Show[
Graphics3D[{Red, PointSize[Large],
Point[{{xp[[j]] + l, yp[[j]] + V[[-1, j]], zp[[j]]}, {0, 0,
0}}], Blue,
Line[{{l, V[[-1, j]], 0}, {xp[[j]] + l, yp[[j]] + V[[-1, j]],
zp[[j]]}}], Tube[{{0, 0, -l}, {0, 0, l Tan[alpha]}}, 0.05],
Black, Line[{{0, 0, l Tan[alpha]}, {l, V[[-1, j]], 0}}]},
PlotRange -> {{0, l + R}, {-l, l}, {-l, l}}, Boxed -> False],
ListLinePlot3D[Table[{sgrid[[i]], V[[i, j]], 0}, {i, nn}],
ColorFunction -> Hue, Axes -> False]], {j, 1, jmax, 2}];

ListAnimate[frame3D2]


Update 1. With initial data {x0, vx0, y0, vy0} = {.0, .01, 0, 0.05} we have Animation

Update 2. In general case we can compute $$\lambda(t,s)$$ as well using integral term in the Lagrange multiplayer definition. The corresponding code is terribly slow due to FindRoot usage. Also this code is unstably at large t. In a case of transvers excitation (see Update 1) we have

(*Parameters *)

l = .815; n = 15; h = l/(n - 1);
EI = .096; g = 9.8;
\[Mu] = 8.41 10^-3; mp = 0.033; R = 0.225; gamma = 2.02;
c = 1.51 10^-2; cp = 5.31 10^-4; alpha = 42/180 Pi;
dt = 5 10^-2; tmax = 15;
sgrid = Range[0, l, h];
tgrid = Range[0, tmax, dt]; jmax = nt = Length[tgrid]; nn =
Length[sgrid]; ds =
Table[NDSolveFiniteDifferenceDerivative[i, sgrid,
"DifferenceOrder" -> 4]["DifferentiationMatrix"], {i, 4}];

U = Array[uu, {nn, nt}]; V = Array[vv, {nn, nt}]; lambda =
Array[ll, {nn, nt}]; lambda0 = Array[ll0, {nt}]; xp =
Array[xx, {nt}]; yp = Array[yy, {nt}]; zp = Array[zz, {nt}]; T1 =
Array[tt1, {nt}]; T2 = Array[tt2, {nt}]; Vpx =
Array[vpx, {nt}]; Vpy =
Array[vpy, {nt}]; fu = -T2 Cos[alpha] + T1 xp/R;(*External force fu*)
fv = -T2 Cos[alpha] V[[-1, All]]/l + T1 yp/R;(*External force fv*)
us = Sqrt[1 - (ds[[1]] . V)^2] - 1;(*inextensibility constraint*)
Do[uu[i, j] =
If[i == 1, 0,
h/2 Sum[(us[[k - 1, j]] + us[[k, j]]), {k, 2, i}]];, {i, nn}, {j,
nt}];(*u[t,s]= Integrate[us,{x,0,s}] *)
Do[ll0[j] = fu[[j]]/(1 + us[[-1, j]]);, {j, nt}];
Do[ll[i, j] =
If[i == 1,
ll0[j], (ll0[j] +
h/2 ( Sum[\[Mu] (uu[k, j] + uu[k, j - 2] - 2 uu[k, j - 1])/
dt^2 + c/dt /2 (uu[k, j] - uu[k, j - 2]), {k, 1,
i - 1}] +
Sum[\[Mu] (uu[k, j] + uu[k, j - 2] - 2 uu[k, j - 1])/dt^2 +
c/dt /2 (uu[k, j] - uu[k, j - 2]), {k, 2, i}]))/(1 +
us[[i, j]])];, {i, nn}, {j, 3, nt}];
(*the Lagrange multiplier \
\[Lambda](t,s)=\[Lambda]0(t)+Integrate[\[Mu] D[u[t,x],t,t]+c \
D[u[t,x],t],{x,s,l}]/(1+D[u[t,s],s])*)
bcrhs[j_] := (1 + 1/2 (ds[[1]] . V[[All, j]])^2) (ds[[2]] .
V[[All, j]])^2 +
ds[[2]] .
V[[All, j]] (2 ds[[1]] . V[[All, j]] ds[[2]] . V[[All, j]] +
4 (ds[[1]] . V[[All, j]])^3 ds[[2]] .
V[[All, j]]) + (1 + (ds[[1]] . V[[All, j]])^2 + (ds[[1]] .
V[[All, j]])^4) ds[[3]] . V[[All, j]];
(*Initial Condition*){x0, vx0, y0, vy0} = {.0, .01, 0,
0.05}; z0 = -Sqrt[R^2 - x0^2 - y0^2];
Do[uu[i, 1] = 0; uu[i, 2] = 0; vv[i, 1] = 0;
vv[i, 2] = 0;, {i, nn}]; Do[xx[i] = x0; vpx[i] = vx0; yy[i] = y0;
vpy[i] = vy0; zz[i] = z0; tt1[i] = -mp g R/z0;
tt2[i] = mp g/Sin[alpha], {i, 2}]; ll0[1] = -tt2[1] Cos[alpha];
ll0[2] = -tt2[2] Cos[alpha]; Do[ll[i, j] = ll0[1], {i, nn}, {j, 2}];
test = Table[.1, {nn}];

Do[(*Verlet explicit step*)
fx = If[j <= 2, 0,
mp (U[[-1, j]] + U[[-1, j - 2]] - 2 U[[-1, j - 1]])/dt^2 +
cp (U[[-1, j]] - U[[-1, j - 1]])/dt ];
fy = If[j <= 2, 0,
mp (V[[-1, j]] + V[[-1, j - 2]] - 2 V[[-1, j - 1]])/dt^2 +
cp (V[[-1, j]] - V[[-1, j - 1]])/dt ];
xx[j + 1] =
xx[j] + vpx[j] dt - (fx + cp vpx[j] + tt1[j] xx[j]/R) dt^2/2 ;
vpx[j + 1] = vpx[j] - (fx + cp vpx[j] + tt1[j] xx[j]/R) dt;
yy[j + 1] =
yy[j] + vpy[j] dt - (fy + cp vpy[j] + tt1[j] yy[j]/R) dt^2/2 ;
vpy[j + 1] = vpy[j] - (fy + cp vpy[j] + tt1[j] yy[j]/R) dt;
zz[j + 1] = -Sqrt[R^2 - xx[j + 1]^2 - yy[j + 1]^2];
tt1[j + 1] = -R/
zz[j + 1] (mp g + (mp (zz[j + 1] + zz[j - 1] - 2 zz[j])/dt^2 +
cp /dt/2 (zz[j + 1] - zz[j - 1])));
tt2[j + 1] = -tt1[j + 1] zz[j + 1]/R/ Sin[alpha];
(*beam implicit step*)
r2 = EI ds[[1]] . bcrhs[j + 1] -
lambda[[All, j + 1]] (ds[[2]] . V[[All, j + 1]]);
eq2 = \[Mu] (V[[All, j + 1]] + V[[All, j - 1]] - 2 V[[All, j]]) +
c dt/2 (V[[All, j + 1]] - V[[All, j - 1]]) + dt^2 r2;
eqs = Table[Drop[Drop[eq2, 2], -2][[i]] == 0, {i, nn - 4}];
varAll = V[[All, j + 1]];
bC0 = {(ds[[1]] . V[[All, j + 1]])[[1]] == 0, V[[1, j + 1]] == 0};
bC1 = {lambda[[-1, j + 1]] (ds[[1]] . V[[All, j + 1]])[[-1]] -
EI bcrhs[j + 1][[-1]] ==
fv[[j + 1]], (ds[[2]] . V[[All, j + 1]])[[-1]] == 0};
eqAll = Join[eqs, bC0, bC1];
ini = Table[{varAll[[i]], test[[i]]}, {i, Length[varAll]}];
sol = Quiet@
FindRoot[eqAll, ini,
Method -> {"Newton", "StepControl" -> "TrustRegion"},
MaxIterations -> 20]; test = varAll /. sol;
Do[vv[i, j + 1] = test[[i]], {i, nn}];, {j, 2,
jmax - 1}] // AbsoluteTiming


It takes about 1210 s on my laptop, finally we have similar pictures as above
Visualization $$\lambda(t,s)$$ (solid lines computed at s=sgrid) together with $$\lambda_0(t)$$ from Update 1 (red dashed line)

• Actually OP has made a mistake when interpreting the system in the paper, $\lambda$ isn't a function of $s$, see my comment above for more info. Commented Jan 25 at 3:15
• @xzczd Thank you for your comments. Actually in general case $\lambda=\lambda(t,s)$ as it defined by equation in my answer. But for this particular case before resonance we take $\lambda=\lambda_0(t)$. To compute resonance we need take into account nonlinear integral terms $\int_0^s(\mu u_{tt}+c u_t)ds$. I have tested this code as well. It terribly slow due to FindRoot usage. Commented Jan 25 at 9:34
• @xzczd Please, see Update 2 to my answer with $\lambda(t,s)$ computation. Commented Jan 26 at 11:16