Mysterious Difference in Equation Formulation for NDSolveValue and NDEigensystem

Question

I have taken a standard equation from Mathematica help from here for a stress operator. I have also formulated the equivalent equations from engineering texts. When I compare the two equations Using FullSimplify I see they are the same. However, when I use them I get different answers. What is happening? Here are the two equations and the comparison using FullSimplify

ps1 = {Inactive[
      Div][({{0, -((Y ν)/(1 - ν^2))}, {-((Y (1 - ν))/(
          2 (1 - ν^2))), 0}}.Inactive[Grad][v[x, y], {x, y}]), {x,
       y}] + Inactive[
      Div][({{-(Y/(1 - ν^2)), 
         0}, {0, -((Y (1 - ν))/(2 (1 - ν^2)))}}.Inactive[
         Grad][u[x, y], {x, y}]), {x, y}],
   Inactive[
      Div][({{0, -((Y (1 - ν))/(2 (1 - ν^2)))}, {-((Y ν)/(
          1 - ν^2)), 0}}.Inactive[Grad][u[x, y], {x, y}]), {x, 
      y}] + Inactive[
      Div][({{-((Y (1 - ν))/(2 (1 - ν^2))), 
         0}, {0, -(Y/(1 - ν^2))}}.Inactive[Grad][
        v[x, y], {x, y}]), {x, y}]};
ps = {Y/(2 (1 + ν)) (D[u[x, y], {x, 2}] + D[u[x, y], {y, 2}]) + 
    Y/(2 (1 - ν)) D[(D[u[x, y], x] + D[v[x, y], y]), x],
   Y/(2 (1 + ν)) (D[v[x, y], {x, 2}] + D[v[x, y], {y, 2}]) + 
    Y/(2 (1 - ν)) D[(D[u[x, y], x] + D[v[x, y], y]), y]};
FullSimplify[Activate[ps1] == -ps]

This gives

*True*

I generate a mesh and then use NDSolveValue and get solutions for both which I plot.

Needs["NDSolve`FEM`"];
Ω = ImplicitRegion[True, {x, y}];
mesh = ToElementMesh[Ω, {{0, 1}, {0, 0.25}}, 
   "MaxCellMeasure" -> 0.001];
{uif1, vif1} = NDSolveValue[{
     ps1 == {0, NeumannValue[1, x == 1]},
     DirichletCondition[u[x, y] == 0, x == 0],
     DirichletCondition[v[x, y] == 0, x == 0]
     } /. {Y -> 10^3, ν -> 33/100},
   {u, v},
   {x, y} ∈ mesh];
{uif, vif} = NDSolveValue[{
     ps == {0, -NeumannValue[1, x == 1]},
     DirichletCondition[u[x, y] == 0, x == 0],
     DirichletCondition[v[x, y] == 0, x == 0]
     } /. {Y -> 10^3, ν -> 33/100},
   {u, v},
   {x, y} ∈ mesh];
Plot3D[{vif1[x, y], vif[x, y]}, {x, y} ∈ mesh, 
 BoxRatios -> {2, 1, 1}, PlotRange -> All]

The two solutions are different why? If I look at the ratio, and ignore diving by zero, we can see that there is a factor of about 1.5 but not a constant value.

   Plot3D[Evaluate[vif1[x, y]/vif[x, y]], {x, y} ∈ mesh, 
     BoxRatios -> {2, 1, 1}, PlotRange -> {All, All, {1, 2}}]

Comparison of the other direction gives a similar ratio. What is happening? Is there some subtle issue over NeumannValues that I am missing? Thanks

EDIT: A bit more

I have also looked at the eigenvalues and vectors. These do not involve the Neumann boundary conditions.

{vals1, vecs1} = 
  NDEigensystem[{ps1, DirichletCondition[u[x, y] == 0, x == 0], 
     DirichletCondition[v[x, y] == 0, x == 0]} /. {Y -> 10^3, \[Nu] ->
       33/100}, {u, v}, {x, y} \[Element] mesh, 10];
{vals, vecs} = 
  NDEigensystem[{-ps, DirichletCondition[u[x, y] == 0, x == 0], 
     DirichletCondition[v[x, y] == 0, x == 0]} /. {Y -> 10^3, \[Nu] ->
       33/100}, {u, v}, {x, y} \[Element] mesh, 10];
TableForm[Transpose[{vals1, vals}], 
 TableHeadings -> {None, {"Help Eqn.", "Textbook Eqn."}}]

The table comparing the eigenvalues gives

where I have called equation ps1 the Help equation and ps the text book equation. Clearly they are very different values. Looking at the first eigenvector I again compare the ratio

{vif1, vif} = {vecs1[[1, 2]], vecs[[1, 2]]};
Plot3D[vif1[x, y]/vif[x, y], {x, y} \[Element] mesh, 
 BoxRatios -> {2, 1, 1}, PlotRange -> {All, All, {0.8, 1.2}}]

These are more similar than the deflection calculation but still significantly different.

I am not sure what to conclude but the stiffness matrices must be different. Is there a good reason for this? Version 10.3 on Windows 7.

Can anyone confirm?

Answer 1

Let's look at a simple example:

pde = Inactive[
   Div][{{0, 1}, {2, 0}}.Inactive[Grad][u[x, y], {x, y}], {x, y}]

the coefficient is

{{0, 1}, {2, 0}}

Now when you activate that you get

Activate[pde]
3*Derivative[1, 1][u][x, y]

This is then what NDSolve parses in the two cases

{state} =
  NDSolve`ProcessEquations[{pde == 0,
    DirichletCondition[u[x, y] == 0, True]},
   u, {x, y} \[Element] Rectangle[]];
state["FiniteElementData"][
  "PDECoefficientData"]["DiffusionCoefficients"]
{{{{0, 1}, {2, 0}}}}

{state} =
  NDSolve`ProcessEquations[{Activate[pde] == 0,
    DirichletCondition[u[x, y] == 0, True]},
   u, {x, y} \[Element] Rectangle[]];
state["FiniteElementData"][
  "PDECoefficientData"]["DiffusionCoefficients"]
{{{{0, 3/2}, {3/2, 0}}}}

So you see this a different PDE model. Inactive is a way to prevent Mathematica to evaluate the PDE too early until there is a chance to parse the PDE before it is evaluated. There is no way to go back from the evaluated coefficients like

{{{{0, 3/2}, {3/2, 0}}}} to {{{{0, 1}, {2, 0}}}}

So, yes, Inactive is needed if either

1) you have unsymmetrical coefficients (which we have for plane stress/strain) That's what we have here or

2) if you need NeumannValue that work with the divergence part of the PDE which is explained here:

FEMDocumentation/tutorial/FiniteElementBestPractice

in the section:

NeumannValue and Formal Partial Differential Equations

More info here and there is a question about this as well. That the inactive stress operator is correct has been discussed here.

I hope this helps.