# Mysterious Difference in Equation Formulation for NDSolveValue and NDEigensystem

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))),
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["NDSolveFEM"];
Ω = 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?

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} =
NDSolveProcessEquations[{pde == 0,
DirichletCondition[u[x, y] == 0, True]},
u, {x, y} \[Element] Rectangle[]];
state["FiniteElementData"][
"PDECoefficientData"]["DiffusionCoefficients"]
{{{{0, 1}, {2, 0}}}}

{state} =
NDSolveProcessEquations[{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