Expression equivalent to zero does not simplify to zero

Question

So I have some code which outputs a sparse matrix of all zeroes and this expression:

(Em h (2-\[Gamma]^2 (-1+\[Nu])))/(12 \[Gamma] (-1+\[Nu]) (1+\[Nu]))+(Em h (-2+\[Gamma]^2 (-1+\[Nu])))/(12 \[Gamma] (-1+\[Nu]) (1+\[Nu]))

The matrix is supposed to be all zeroes, though. So I checked the above expression by hand and it is exactly equal to zero.

Furthermore, running Simplify on this expression in a separate notebook shows it evaluates to zero:

Print[Simplify[(Em h (2 - \[Gamma]^2 (-1 + \[Nu])))/(
12 \[Gamma] (-1 + \[Nu]) (1 + \[Nu])) + (
Em h (-2 + \[Gamma]^2 (-1 + \[Nu])))/(
12 \[Gamma] (-1 + \[Nu]) (1 + \[Nu]))]];

This is strange because when I try Simplify on the full matrix in the original notebook, it does not evaluate to zero, it just stays as the expression. This inconsistency holds for Simplify, FullSimplify, and Expand.

Full matrix which is supposed to be all zeroes:

{{0,0,0,0,(Em h (2-\[Gamma]^2 (-1+\[Nu])))/(12 \[Gamma] (-1+\[Nu]) (1+\[Nu]))+(Em h (-2+\[Gamma]^2 (-1+\[Nu])))/(12 \[Gamma] (-1+\[Nu]) (1+\[Nu])),0,0,0},{0,0,0,0,0,0,0,0},{0,0,0,0,0,0,(Em h (2-\[Gamma]^2 (-1+\[Nu])))/(12 \[Gamma] (-1+\[Nu]) (1+\[Nu]))+(Em h (-2+\[Gamma]^2 (-1+\[Nu])))/(12 \[Gamma] (-1+\[Nu]) (1+\[Nu])),0},{0,0,0,0,0,0,0,0},{(Em h (2-\[Gamma]^2 (-1+\[Nu])))/(12 \[Gamma] (-1+\[Nu]) (1+\[Nu]))+(Em h (-2+\[Gamma]^2 (-1+\[Nu])))/(12 \[Gamma] (-1+\[Nu]) (1+\[Nu])),0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0},{0,0,(Em h (2-\[Gamma]^2 (-1+\[Nu])))/(12 \[Gamma] (-1+\[Nu]) (1+\[Nu]))+(Em h (-2+\[Gamma]^2 (-1+\[Nu])))/(12 \[Gamma] (-1+\[Nu]) (1+\[Nu])),0,0,0,0,0},{0,0,0,0,0,0,0,0}}

Any ideas?

Here's the full notebook I'm working:

ClearAll[Evaluate[Context[] <> "*"]];

(* From Figure 17.1 *)

Quad4IsoPShapeFunDer[encoor_, qcoor_] := 
  Module[{Nf, dNx, dNy, dN\[Xi], dN\[Eta], i, J11, J12, J21, J22, 
    Jdet, \[Xi], \[Eta], x, y, x1, x2, x3, x4, y1, y2, y3, 
    y4}, {\[Xi], \[Eta]} = 
    qcoor; {{x1, y1}, {x2, y2}, {x3, y3}, {x4, y4}} = encoor;
   Nf = {(1 - \[Xi])*(1 - \[Eta]), (1 + \[Xi])*(1 - \[Eta]), (1 + \
\[Xi])*(1 + \[Eta]), (1 - \[Xi])*(1 + \[Eta])}/4;
   dN\[Xi] = {-(1 - \[Eta]), (1 - \[Eta]), (1 + \[Eta]), -(1 + \
\[Eta])}/4;
   dN\[Eta] = {-(1 - \[Xi]), -(1 + \[Xi]), (1 + \[Xi]), (1 - \[Xi])}/4;
   x = {x1, x2, x3, x4}; y = {y1, y2, y3, y4};
   J11 = dN\[Xi].x; J12 = dN\[Xi].y; J21 = dN\[Eta].x; 
   J22 = dN\[Eta].y;
   Jdet = Simplify[J11*J22 - J12*J21];
   dNx = (J22*dN\[Xi] - J12*dN\[Eta])/Jdet; dNx = Simplify[dNx];
   dNy = (-J21*dN\[Xi] + J11*dN\[Eta])/Jdet; dNy = Simplify[dNy];
   Return[{Nf, dNx, dNy, Jdet}]
   ];

(* From Figure 17.3 *)

LineGaussRuleInfo[{rule_, numer_}, point_] := 
  Module[{g2 = {-1, 1}/Sqrt[3], g3 = {-Sqrt[3/5], 0, Sqrt[3/5]}, 
    w3 = {5/9, 8/9, 5/9}, 
    w4 = {(1/2) - Sqrt[5/6]/6, (1/2) + Sqrt[5/6]/6, (1/2) + 
       Sqrt[5/6]/6, (1/2) - Sqrt[5/6]/6}, 
    g4 = {-Sqrt[(3 + 2*Sqrt[6/5])/7], -Sqrt[(3 - 2*Sqrt[6/5])/7], 
      Sqrt[(3 - 2*Sqrt[6/5])/7], Sqrt[(3 + 2*Sqrt[6/5])/7]}, 
    g5 = {-Sqrt[5 + 2*Sqrt[10/7]], -Sqrt[5 - 2*Sqrt[10/7]], 0, 
       Sqrt[5 - 2*Sqrt[10/7]], Sqrt[5 + 2*Sqrt[10/7]]}/3, 
    w5 = {322 - 13*Sqrt[70], 322 + 13*Sqrt[70], 512, 
       322 + 13*Sqrt[70], 322 - 13*Sqrt[70]}/900, i = point, p = rule,
     info = {Null, 0}}, If[p == 1, info = {0, 2}];
   If[p == 2, info = {g2[[i]], 1}];
   If[p == 3, info = {g3[[i]], w3[[i]]}];
   If[p == 4, info = {g4[[i]], w4[[i]]}];
   If[p == 5, info = {g5[[i]], w5[[i]]}];
   If[numer, Return[N[info]], Return[Simplify[info]]];
   ];

(* From Figure 17.5 *)

QuadGaussRuleInfo[{rule_, numer_}, point_] := 
  Module[{\[Xi], \[Eta], p1, p2, i, j, w1, w2, m, 
    info = {{Null, Null}, 0}},
   If[Length[rule] == 2, {p1, p2} = rule, p1 = p2 = rule];
   If[p1 < 0, 
    Return[QuadNonProductGaussRuleInfo[{-p1, numer}, point]]];
   If[Length[point] == 2, {i, j} = point, m = point;
    j = Floor[(m - 1)/p1] + 1; i = m - p1*(j - 1)];
   {\[Xi], w1} = LineGaussRuleInfo[{p1, numer}, i];
   {\[Eta], w2} = LineGaussRuleInfo[{p2, numer}, j];
   info = {{\[Xi], \[Eta]}, w1*w2};
   If[numer, Return[N[info]], Return[Simplify[info]]];
   ];

(* From Figure 17.8 *)

Quad4IsoPMembStiffness[encoor_, Emat_, th_, options_] := 
  Module[{i, k, pG = 2, numer = False, h = th, qcoor, c, w, Nf, dNx, 
    dNy, Jdet, Be, Ke = Table[0, {8}, {8}]},
   If[Length[options] == 2, {numer, pG} = options, {numer} = 
     options];
   If[pG < 1 || pG > 4, Print["pG out of range"]; Return[Null]];
   For[k = 1, k <= pG*pG, 
    k++, {qcoor, w} = QuadGaussRuleInfo[{pG, numer}, k];
    {Nf, dNx, dNy, Jdet} = Quad4IsoPShapeFunDer[encoor, qcoor];
    If[Length[th] == 4, h = th.Nf];
    c = w*Jdet*h;
    Be = {Flatten[Table[{dNx[[i]], 0}, {i, 4}]], 
      Flatten[Table[{0, dNy[[i]]}, {i, 4}]], 
      Flatten[Table[{dNy[[i]], dNx[[i]]}, {i, 4}]]};
    Ke += Simplify[c*Transpose[Be].(Emat.Be)];
    ];
   Return[Simplify[Ke]]
   ];

b = a/\[Gamma];
ncoor = {{0, 0}, {a, 0}, {a, b}, {0, b}};
Emat = Em/(1 - \[Nu]^2)*{{1, \[Nu], 0}, {\[Nu], 1, 0}, {0, 
     0, (1 - \[Nu])/2}};
Ke2 = Quad4IsoPMembStiffness[ncoor, Emat, h, {False, 2}];
scaledKe2 = Simplify[Ke2*(24*(1 - \[Nu]^2)*\[Gamma]/(Em*h))];
Print["For p = 2, Ke = Ke2 = ", Em*h/(24*\[Gamma]*(1 - \[Nu]^2)), 
  "*\n", scaledKe2 // MatrixForm];
Print["For p = 2, the eigenvalues of Ke2 are: ", Eigenvalues[Ke2]];

ClearAll[Em, \[Nu], a, b, h, \[Gamma]];

b = a/\[Gamma];
ncoor = {{0, 0}, {a, 0}, {a, b}, {0, b}};
Emat = Em/(1 - \[Nu]^2)*{{1, \[Nu], 0}, {\[Nu], 1, 0}, {0, 
     0, (1 - \[Nu])/2}};
Ke3 = Quad4IsoPMembStiffness[ncoor, Emat, h, {False, 3}];
scaledKe3 = Simplify[Ke3*(24*(1 - \[Nu]^2)*\[Gamma]/(Em*h))];
Print["For p = 3, Ke = Ke3 = ", Em*h/(24*\[Gamma]*(1 - \[Nu]^2)), 
  "*\n", scaledKe3 // MatrixForm];
Print["For p = 3, the eigenvalues of Ke3 are: ", Eigenvalues[Ke3]];

ClearAll[Em, \[Nu], a, b, h, \[Gamma]];


b = a/\[Gamma];
ncoor = {{0, 0}, {a, 0}, {a, b}, {0, b}};
Emat = Em/(1 - \[Nu]^2)*{{1, \[Nu], 0}, {\[Nu], 1, 0}, {0, 
     0, (1 - \[Nu])/2}};
Ke4 = Quad4IsoPMembStiffness[ncoor, Emat, h, {False, 4}];
Ke4=Simplify[Ke4];
scaledKe4 = Simplify[Ke4*(24*(1 - \[Nu]^2)*\[Gamma]/(Em*h))];
Print["For p = 4, Ke = Ke4 = ", Em*h/(24*\[Gamma]*(1 - \[Nu]^2)), 
  "*\n", scaledKe4 // MatrixForm];
Print["For p = 4, the eigenvalues of Ke4 are: ", Eigenvalues[Ke4]];

Print["As shown above, where p=2,3,4, the stiffness matrix does not \
change and three Eigenvalues are zero" ];

(* Ke2, Ke3, Ke4 are identical if inspected item by item. Yet the following evaluates to False *)
If[Simplify[Ke2 == Ke3] && Simplify[Ke2 == Ke4], Print["Ke2 = Ke3 = Ke4"] ];

(* The following should evaluate to zero, but Ke3 - Ke4 spits out the problem expression noted above (which itself is exactly equal to zero, but Mathematica doesn't want to evaluate it to zero for some reason) *)
Print[Simplify[Ke2 - Ke3]];
Print[Simplify[Ke3 - Ke4]];

Answer 1

I believe your question to be: "If Ke3 and Ke4 are fully simplified and equal to each other, why isn't Ke3-Ke4 the zero matrix?". In which case the answer: "because there is no universal most simplified form of any expression". Where you end up depends on where you started. So just because Simplify[a==b] is True does not mean that, in general, Simplify[a]==Simplify[b]. It often is, but there are no guarantees. If you look at the entries of Ke3 and Ke4 you'll notice that they differ in placements of minus signs in certain numerators and denominators. This is enough to prevent automatic evaluation from producing zero. In this case, you only need a very cheap form of simplification: putting things over a common denominator. So Ke3 - Ke4 // Together evaluates to the zero matrix. But just plain subtraction or equality will not.