# Free vibrations of square plate with bi-harmonic equation

I am still playing with the wave bi-harmonic equation. My question is very similar to another question I posted few weeks ago on this site (Solving rectangular plates vibrations wave equation). I want to visualize the free vibrations of a steel square plate measuring [-a,a;-b,b] and thickness = h (origin = center of the plate).

The issue I met today is the following one : the bi-harmonic equation I am using to NDSolve my problem is working only with some initial conditions parameters sets. With other sets, the computation is never ending (or maybe too long, I don't know) and I get error messages like "DSolve::mxsst: Using maximum number of grid points 100 allowed by the MaxPoints or MinStepSize options for independent variable x." Idem for y and z.

Of course, I searched and found some explanations on this site and some others, but the given solutions never work for me. I think my problem is a method tuning matter in NDSolve. Too hard for me. Help. Thanks in advance.

Below is my code with some helping comments :

(* Free vibrations of square plate *)
Remove["Global*"]; // Quiet
Needs["DifferentialEquationsNDSolveProblems"];
Needs["DifferentialEquationsNDSolveUtilities"];
Needs["NDSolveFEM"];
Ey = SetPrecision[2.078*10^11, Infinity];(*N/m^2*)
\[Nu] = Rationalize[0.317756];(*unitless*)
\[CurlyRho]d = 8166;(*kg/m^3*)
h = 1/1000;
Df = (Ey h^3)/(12 (1 - \[Nu]^2)); a = 1; b = 1;
eqn = {D[u[x, y, t], {x, 4}] + 2 D[u[x, y, t], {x, 2}, {y, 2}] +
D[u[x, y, t], {y, 4}] + (\[CurlyRho]d h)/
Df D[u[x, y, t], {t, 2}] ==
0}; (* for plates, governing equation = bi-harmonic equation *)
m = 1; n = 2;
ic = {u[x, y, 0] ==
Cos[m \[Pi] x/a] Cos[n \[Pi] y/b] +
Cos[n \[Pi] x/a] Cos[m \[Pi]  y/b],
Derivative[0, 0, 1][u][x, y, 0] == 0}; (* initial impulse *)
bc1fe = {Derivative[2, 0, 0][u][-a, y, t] == 0,
Derivative[3, 0, 0][u][-a, y, t] == 0,
Derivative[2, 0, 0][u][a, y, t] == 0,
Derivative[3, 0, 0][u][a, y, t] ==
0}; (* free edges/ simplified *)
bc2fe = {Derivative[0, 2, 0][u][x, -b, t] == 0,
Derivative[0, 3, 0][u][x, -b, t] == 0,
Derivative[0, 2, 0][u][x, b, t] == 0,
Derivative[0, 3, 0][u][x, b, t] ==
0}; (* free edges/ simplified *)
numsol = NDSolve[{eqn, ic, bc1fe, bc2fe},
u, {x, -a, a}, {y, -b, b}, {t, 0, 1},
Method -> {"StiffnessSwitching"},
MaxSteps -> 10^6] (* must show the plate behavior in time *)
(* Results visualization *)
uu[x_, y_, t_] := u[x, y, t] /. numsol;
ListAnimate[
Table[ContourPlot[Evaluate[uu[x, y, t] == 0], {x, -a, a}, {y, -b, b},
PlotLabel -> Style[t, 30, Bold]], {t, 0, 1, .05}]]

*)
(* I also tried the following possibilities without any success : *)
(* numsol=NDSolve[{eqn,ic,bc1fe,bc2fe},u,{x,-a,a},{y,-b,b},{t,0,1},\
Method\[Rule]{"MethodOfLines","SpatialDiscretization"\[Rule]{\
"TensorProductGrid","MaxPoints"\[Rule]101}}]*)

(* numsol=NDSolve[{eqn,ic,bc1fe,bc2fe},u,{x,-a,a},{y,-b,b},{t,0,1},\
Method\[Rule]{"MethodOfLines","TemporalVariable"\[Rule]t,\
"SpatialDiscretization"\[Rule]"FiniteElement"}] *)

(* one important thing : if I replace " Cos[m \[Pi] x/a] Cos[n \[Pi] \
y/b]+Cos[n \[Pi] x/a] Cos[m \[Pi](\y)/b]" by " Cos[1 x] Cos[3 \
y]+Cos[3x] Cos[1y] " the resolution is working well and relatively fast. *)
(* No error messages hapening *)


Hereafter is the computation result for Cos[1 x] Cos[3y]+Cos[3x] Cos[1y]. It is also working well for Cos[1 x] Cos[2y]+Cos[2x] Cos[1y]:

• Remove["Global*"]; // Quiet; Needs["DifferentialEquationsNDSolveProblems"]; Needs["DifferentialEquationsNDSolveUtilities"]; Needs["NDSolveFEM"]; are redundant. Also, the FiniteElement method isn't called in this case, because it cannot handle this type of problem directly at the moment: mathematica.stackexchange.com/a/199369/1871 May 20, 2021 at 6:37
• @ xzczd : thnak you. May 20, 2021 at 13:54

Several issues here:

1. Since ibcinc warning pops up and you b.c.s all involve derivative, manual adjustion of "DifferentiateBoundaryConditions" option discussed in this post is necessary.

2. Certain silent change (backslide?) has happened on the ODE solver of NDSolve in recent versions of Mathematica. When the spatial grid is not too dense, v8.0.4 can handle the problem without adjustion of the ODE solver:

but it's not the case for v12.2. One possible work-around I find is to explicitly set Method -> Adams.

The following is the fixed code for v12.2. jianshi is modified from the tools in this post:

showStatus[status_]:=LinkWrite[\$ParentLink,
SetNotebookStatusLine[FrontEndEvaluationNotebook[],ToString[status]]];

clearStatus[]:=showStatus[""];
clearStatus[]

jianshi[t_]:=EvaluationMonitor:>showStatus["t = "<>ToString[CForm[t]]]

mol[n:_Integer|{_Integer..}, o_:"Pseudospectral"] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,
"MinPoints" -> n, "DifferenceOrder" -> o}}

mol[tf:False|True,sf_:Automatic]:={"MethodOfLines",
"DifferentiateBoundaryConditions"->{tf,"ScaleFactor"->sf}}

numsol = NDSolveValue[{eqn, ic, bc1fe, bc2fe}, u, {x, -a, a}, {y, -b, b}, {t, 0, 1},
Method -> Union[mol[25, 4], mol[True, 100], {Method -> Adams}], jianshi[t],
MaxSteps -> Infinity]; // AbsoluteTiming
(* {61.7096, Null} *)

Table[
Plot3D[numsol[x, y, t], {x, -a, a}, {y, -b, b}, PlotLabel -> Style[t, 30, Bold],
PlotRange -> All], {t, 0, 1, .1}] // Partition[#, 5]& // Grid


Don't be worried about the eerri and eerr warnings, they're merely warnings in this case.

• @ xzczd, thank you for your help. It is very complex for me, but I will try to understand. :) May 20, 2021 at 13:53
• @ xzczd : I am impressed ! I applied your solution and the solving time was very very short : {95.3622, Null}. May 20, 2021 at 15:49
• @ xzczd : I begin to (try to) understand. Your solution is great and multiple. It's amazing. Thank you. I really appreciate. The secret of NDSolve is hiding in "Method" settings. May I ask you a question : how did you get all your knowledge about this? Because, personnally I feel it very very tricky, isn't it? May 21, 2021 at 8:26
• @Pascal77 The options adjusted in this post are mentioned in the tutorial The Numerical Method of Lines, but I admit it's obscure, I myself dare not say I've fully understood it even now. Things help me most are the answers in this site, for example, I had probably given up using NDSolve` if I didn't receive user21's answer here 9 years ago. Knowledge about PDE and core language of Mathematica also helps me to understand the tutorial. The rest is trial & error. May 21, 2021 at 10:42
• @ xzczd : you practiced Mathematica for years. Congratulations for your perseverance. I admit that it is fascinating. Personally, I recently discovered Mathematica (november 2020). I also tried Maple, but I prefer Mathematica. So, my knowledge is very low, and I try to learn through subjects which are near and dear to me, like "Chladni simulation". I began with strange attractors and some electronic circuits modelization. At last, I am trying to understand more deeply Mathematica. In any case, it is a pleasure to discuss with you. May 21, 2021 at 18:09