0
$\begingroup$

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]];
$\endgroup$
  • $\begingroup$ "on the full matrix" - where is this matrix? $\endgroup$ – J. M. is away Nov 1 '17 at 5:03
  • $\begingroup$ @J.M. I added it to the post $\endgroup$ – Casey Baker Nov 1 '17 at 5:06
  • $\begingroup$ What version are you on? Applying Simplify[] to that matrix you posted does give the zero matrix in 11.2. $\endgroup$ – J. M. is away Nov 1 '17 at 5:08
  • $\begingroup$ @J.M. I'm on 11.2 as well. And yes, that's my point. Applying Simplify[] to the matrix by itself outside of my code gives zero. But applying simplify[] within my code doesn't change it at all and outputs the expression instead. $\endgroup$ – Casey Baker Nov 1 '17 at 5:19
  • $\begingroup$ So your troublesome matrix is Ke3 - Ke4? Simplify[Ke3 - Ke4] comes out all zeros for me. I assume you've tried restarting your kernel. $\endgroup$ – aardvark2012 Nov 1 '17 at 5:28
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ good answer, thanks $\endgroup$ – Casey Baker Nov 4 '17 at 5:49

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.