6
$\begingroup$

Having installed version 11 I thought I would check an old Stack Exchange vibration problem using NDEigensytem. The old problem was Test a wooden board's vibration mode and I think this was before NDEigensystem and thus very complicated. However I get different answers to the old question and a warning I don't understand.

User21 had done all the hard work in the previous question part of which I copy now to make the mesh. The code is as follows. First define the stress operator for 3D structures

Needs["NDSolve`FEM`"]

ClearAll[stressOperator, u, v, w, x, y, z];
stressOperator[
  Y_, \[Nu]_] := {Inactive[
     Div][{{0, 0, -((Y*\[Nu])/((1 - 2*\[Nu])*(1 + \[Nu])))}, {0, 0, 
       0}, {-Y/(2*(1 + \[Nu])), 0, 0}}.Inactive[Grad][
      w[x, y, z], {x, y, z}], {x, y, z}] + 
   Inactive[
     Div][{{0, -((Y*\[Nu])/((1 - 2*\[Nu])*(1 + \[Nu]))), 
       0}, {-Y/(2*(1 + \[Nu])), 0, 0}, {0, 0, 0}}.Inactive[Grad][
      v[x, y, z], {x, y, z}], {x, y, z}] + 
   Inactive[
     Div][{{-((Y*(1 - \[Nu]))/((1 - 2*\[Nu])*(1 + \[Nu]))), 0, 
       0}, {0, -Y/(2*(1 + \[Nu])), 0}, {0, 
       0, -Y/(2*(1 + \[Nu]))}}.Inactive[Grad][
      u[x, y, z], {x, y, z}], {x, y, z}], 
  Inactive[Div][{{0, 0, 0}, {0, 
       0, -((Y*\[Nu])/((1 - 
              2*\[Nu])*(1 + \[Nu])))}, {0, -Y/(2*(1 + \[Nu])), 
       0}}.Inactive[Grad][w[x, y, z], {x, y, z}], {x, y, z}] + 
   Inactive[
     Div][{{0, -Y/(2*(1 + \[Nu])), 
       0}, {-((Y*\[Nu])/((1 - 2*\[Nu])*(1 + \[Nu]))), 0, 0}, {0, 0, 
       0}}.Inactive[Grad][u[x, y, z], {x, y, z}], {x, y, z}] + 
   Inactive[
     Div][{{-Y/(2*(1 + \[Nu])), 0, 
       0}, {0, -((Y*(1 - \[Nu]))/((1 - 2*\[Nu])*(1 + \[Nu]))), 0}, {0,
        0, -Y/(2*(1 + \[Nu]))}}.Inactive[Grad][
      v[x, y, z], {x, y, z}], {x, y, z}], 
  Inactive[Div][{{0, 0, 0}, {0, 
       0, -Y/(2*(1 + \[Nu]))}, {0, -((Y*\[Nu])/((1 - 
              2*\[Nu])*(1 + \[Nu]))), 0}}.Inactive[Grad][
      v[x, y, z], {x, y, z}], {x, y, z}] + 
   Inactive[
     Div][{{0, 0, -Y/(2*(1 + \[Nu]))}, {0, 0, 
       0}, {-((Y*\[Nu])/((1 - 2*\[Nu])*(1 + \[Nu]))), 0, 0}}.Inactive[
       Grad][u[x, y, z], {x, y, z}], {x, y, z}] + 
   Inactive[
     Div][{{-Y/(2*(1 + \[Nu])), 0, 0}, {0, -Y/(2*(1 + \[Nu])), 0}, {0,
        0, -((Y*(1 - \[Nu]))/((1 - 2*\[Nu])*(1 + \[Nu])))}}.Inactive[
       Grad][w[x, y, z], {x, y, z}], {x, y, z}]}

Next set up the geometry.

base = {0, 0, 0};
h1 = 5;
h2 = 5;
w1 = 40;
l1 = 76;
cw1 = 5;
cl1 = 68;
cw2 = 36;
cl2 = 5;
offset1 = base + {(w1 - cw1)/2, (l1 - cl1)/2, 0};
offset2 = base + {(w1 - cw2)/2, (l1 - cl2)/2, 0};
offset3 = base + {(w1 - cw1)/2, (l1 - cl2)/2, 0};
ClearAll[rect]
rect[base_, w_, l_, h_] := {base + {0, 0, h}, base + {w, 0, h}, 
  base + {w, l, h}, base + {0, l, h}}
coords = ConstantArray[{0., 0., 0.}, 4 + 4 + 12 + 12];
coords[[{1, 2, 3, 4}]] = rect[base, w1, l1, 0];
coords[[{5, 6, 7, 8}]] = rect[base, w1, l1, h1];
coords[[{9, 10, 15, 16}]] = rect[offset1, cw1, cl1, h1];
coords[[{19, 12, 13, 18}]] = rect[offset2, cw2, cl2, h1];
coords[[{20, 11, 14, 17}]] = rect[offset3, cw1, cl2, h1];
coords[[20 + Range[12]]] = ({0, 0, h2} + #) & /@ 
   coords[[8 + Range[12]]];
bmesh = ToBoundaryMesh["Coordinates" -> coords, 
   "BoundaryElements" -> {QuadElement[{{1, 2, 3, 4}, {1, 2, 6, 5}, {2,
         3, 7, 6}, {3, 4, 8, 7}, {4, 1, 5, 8}, {5, 6, 10, 9}, {6, 12, 
        11, 10}, {6, 7, 13, 12}, {7, 15, 14, 13}, {7, 8, 16, 15}, {8, 
        18, 17, 16}, {8, 5, 19, 18}, {5, 9, 20, 19}, 
       Sequence @@ ({{9, 10, 11, 20}, {11, 12, 13, 14}, {14, 15, 16, 
            17}, {17, 18, 19, 20}, {20, 11, 14, 17}} + 12), 
       Sequence @@ (Partition[Join[Range[9, 20]], 2, 1, 
           1] /. {i1_, i2_} :> {i1, i2, i2 + 12, i1 + 12})}]}];

Look at the geometry and then the mesh

Show[bmesh["Wireframe"], 
 bmesh["Wireframe"["MeshElement" -> "PointElements", 
   "MeshElementIDStyle" -> Red]]]

Mathematica graphics

mesh = ToElementMesh[bmesh, "MeshOrder" -> 1, "MaxCellMeasure" -> 10];
mesh["Wireframe"]

Mathematica graphics

Now we use NDEigensystem to find the eigenvalues and vectors

{vals, vecs} = 
  NDEigensystem[
   stressOperator[100, 1/3], {u, v, w}, {x, y, z} \[Element] mesh, 
   14];

This gives the warning

(* Eigensystem::chnpdef: Warning: there is a possibility that the second matrix SparseArray[<<1>>] in the first argument is not positive definite, which is necessary for the Arnoldi method to give accurate results. *)

So something is suspect.

The eigenvalues are as follows

TableForm[vals, TableHeadings -> {Automatic, None}]

Mathematica graphics

The first six are zero which is correct since these are the rigid body modes. The seventh is negative which is very suspect. The next values are possible but are different to the values found in the original question which were

{0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.011403583383327644`, \
0.01526089137692353`, 0.05661022352859022`, 0.07266104128273859`}

Although they may just about agree to a couple of figures. The eigenvectors are

Column@(MeshRegion[
     ElementMeshDeformation[mesh, vecs[[#]], 
      "ScalingFactor" -> -300]] & /@ Range[7, 14])

Mathematica graphics

I also tried to run the code from the previous question directly and I could not get it to work.

In summary what does the warning mean from NDEigensystem and how can this be avoided? Also, it looks like the output from either the previous answer or this one is wrong. Finally, should the previous code work in version 11? Thanks

$\endgroup$
  • $\begingroup$ Are you by any chance coming to the Wolfram Tech conference this year? Might be fun to meet up. $\endgroup$ – user21 Sep 6 '16 at 19:58
  • $\begingroup$ I would like to come and a meet up would be fun but other duties and finance mean I can't. Pity. $\endgroup$ – Hugh Sep 6 '16 at 21:38
  • $\begingroup$ Sure I understand. If you have nice FEM applications you can share I am always curious of what customers do - to better understand in which directions the product needs to be developed. $\endgroup$ – user21 Sep 6 '16 at 21:47
4
$\begingroup$

Looking at the post you mention, you'll see that Eigensystem is called with no method option. That means that Eigensystem will use a "Direct" solver.

For NDEigensystem the default eigensolver is "Arnoldi" - which is much more efficient than the "Direct" solver but has disadvantages. To use the "Direct" solver from NDEigenststem you can call:

{vals, vecs} = 
 NDEigensystem[
  stressOperator[100, 1/3], {u, v, w}, {x, y, z} \[Element] mesh, 14, 
  Method -> "Direct"]

This will be slow but is reliable and always a good option to cross-check results (perhaps on a smaller mesh). One other alternative is to use:

{vals, vecs} = 
 NDEigensystem[
  stressOperator[100, 1/3], {u, v, w}, {x, y, z} \[Element] mesh, 14, 
  Method -> {"FEAST", "Tolerance" -> 10^-6}]

Which is better and probably the best option for this problem. Note that the original post has the same issue. So nothing really changed.

If one looks at the ref page of Eigenvalue you'll find a small note that states:

"If they are numeric, eigenvalues are sorted in order of decreasing absolute value."

I think that's an unfortunate design decision, but it is what it is. In this example we can assume real eigenvalues so we can go ahead and do:

{vals, vecs} = 
 NDEigensystem[
  stressOperator[100, 1/3], {u, v, w}, {x, y, z} \[Element] mesh, 14, 
  Method -> {"Arnoldi", "Criteria" -> "RealPart"}]

This will still give a message but the spurious negative eigenvalue is not relevant any longer but, on the downside, the duplicate 0 eigenvalues get lost too and the result in not terribly accurate.

To understand the message it help to understand how NDEigensytem works. It transforms the equations given into time dependent equations, uses dummy initial conditions and a dummy time integration range. It never actually performs a time integration. But the transient system of PDEs is assembled like this:

Needs["NDSolve`FEM`"]
{dpde, dbcs, vd, sd, md} = 
  ProcessPDEEquations[{Thread[
     D[#[t, x, y, z] & /@ {u, v, w}, t] == 
       stressOperator[100, 1/3] /. 
      Map[#[x, y, z] -> #[t, x, y, z] &, {u, v, w}]], #[0, x, y, z] ==
        0 & /@ {u, v, w}}, {u, v, w}, {t, 0, 1}, {x, y, z} \[Element] 
    mesh];

It then extracts the stiffness and damping matrix and does the eigen values analysis on those matrices like so:

{load, stiff, damp, mass} = dpde["SystemMatrices"];
Eigensystem[{stiff, damp}, -10, Method -> "Arnoldi"];

This give the same message.

Note that

PositiveDefiniteMatrixQ[damp]
True

So Eigensystem's Arnoldi is a bit of a feisty fellow.

Concerning the code from the other post not working, I can not confirm that; perhaps you forgot to adjust the numEigenToCompute to be larger then 7? In any case I'll add a note for that. If there is anything else that does not work could you clarify that a bit?

$\endgroup$
  • $\begingroup$ Many thanks for your fast response. The two methods you suggest do give consistent results with regard to the original post. With method direct it takes 775 seconds and with method FEAST it takes 23 seconds. $\endgroup$ – Hugh Sep 6 '16 at 21:37
  • $\begingroup$ Ultimately, it were good if Eigensystem had a bit more reliable "Arnoldi" method or something similar, but I would not hold my breath for that.... $\endgroup$ – user21 Sep 6 '16 at 21:45
  • $\begingroup$ I have consolidated my issues with the original question and put them as an answer there If you could explain what goes wrong that would be kind and helful. Thanks $\endgroup$ – Hugh Sep 7 '16 at 10:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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