(************************************************************************* This program makes a 210 x 210 matrix and verifies that it has nonzero determinant. Below is a run of the program, and one may verify that the determinant computed is a rational number whose numerator is -1 and whose denominator is a 302 digit positive integer. In particular, the determinant is nonzero. ----------------------------------------------------------------- adams@obsidian % math Mathematica 5.0 for Linux Copyright 1988-2003 Wolfram Research, Inc. -- Motif graphics initialized -- In[1]:= < 121455478350865788963744826837190849488014250533056686482842058460763\ > 182801234658814997671716912459939840000000000000000000000000000000000\ > 000000000000000000000000000000000000000000000000000000000000000000000\ > 00000000000000000000000000000000) In[2]:= adams@obsidian % ----------------------------------------------------------------- Explanation of the connection between the 210 x 210 matrix and the 162 degrees of freedom specified in the paper: Let \R = the set of real numbers. First, 210 = 6 x 35 = 6 x (7 choose 3) is the dimension of the vector space V of \S-valued polynomials in three variables of degree \le 4, where \S is the 6 dimensional space of 3 x 3 symmetric matrices. Row divergence rdiv carries V into the vector space W of \R^3-valued polynomials in three variables of degree \le 3. The dimension of W is 3 x (6 choose 3) = 3 x 20 = 60. Within W is the vector space X of \R^3-valued polynomials in three variables of degree \le 1. The dimension of X is 3 x (4 choose 1)= 3 x 4 = 12. The preimage Y of X under rdiv : V --> W has dimension 210 - (60 - 12) = 210 - 48 = 162 and this is the space on which the 162 degrees of freedom are specified in the paper. Note that Y is a vector subspace of V. The 162 degrees of freedom all extend from Y to V, and, putting them all together, we get a map F : V --> \R^{162}. The goal is then to show that F | Y : Y --> \R^{162} is a bijective linear map, i.e., that the degrees of freedom on Y are "unisolvent". Let Z = \R^{48} and let f : V --> X be the map which, given a polynomial v in V computes the 48 coefficients of the quadratic and cubic terms in rdiv(v). This map is surjective because rdiv : V --> W is surjective. Moreover, the kernel of f is the 162-dimensional space Y. Then F | Y : Y --> \R^{162} is bijective iff F x f : V --> \R^{162} x \R^{48} = \R^{210} is. Fixing the basis of monomials of V, we identify V with \R^{210} and identify F x f with a linear map \R^{210} --> \R^{210}. The matrix of this map is 210 x 210, and the program verifies that this matrix is invertible. ----------------------------------------------------------------- Some words about the program itself: The program begins by making a 3 x 3 symmetric matrix M whose i,j entry is a generic polynomial of degree \le 4 in three variables x,y,z. The coefficient on x^k * y^l * z^m in the (i,j) entry is denoted a[i,j,k,l,m]. Note that, as defined, a[i,j,k,l,m] is symmetric in i and j, so the matrix M is symmetric. It then sets vars to be a flattened list of all the a[i,j,k,l,m]. This list has 210 entries, if one takes into account symmetry in i and j. These 210 variables a[i,j,k,l,m] are to be thought of as unknowns and we wish to place 210 conditions on them. We eventually calculate the determinant of the resulting 210 x 210 matrix. divM is the row divergence of M. It is a 3-vector of polynomials of degree \le 3 in x,y,z. If i is 1,2 or 3, if k,l,m are all 0,1,2 or 3, then the coefficient on the x^k * y^l * z^m term in the i-th entry of divM is a linear combination of the 210 variables a[i,j,k,l,m]. quotientdivlinear is the collection of coefficients of entries of divM that are on terms of the form x^k * y^l * z^m with k + l + m > 1. This is then a collection of 48 linear combinations of the 210 variables a[i,j,k,l,m] and the simultaneous vanishing of all these linear combinations is equivalent to the requirement that all three entries of divM are of degree \le 1. alphatetrahedron is a collection of 6 linear combinations of the 210 variables a[i,j,k,l,m] whose vanishing is equivalent to the vanishing of the chamber degrees of freedom of M. alphafaces is a collection of 36 linear combinations of the 210 variables a[i,j,k,l,m] whose vanishing is equivalent to the vanishing of the facet degrees of freedom of M. alphaedges is a collection of 96 linear combinations of the 210 variables a[i,j,k,l,m] whose vanishing is equivalent to the vanishing of the edge degrees of freedom of M. alphavertices is a collection of 24 linear combinations of the 210 variables a[i,j,k,l,m] whose vanishing is equivalent to the vanishing of the vertex degrees of freedom of M. alpha is the 48+6+36+96+24=210 linear combinations of the 210 variables a[i,j,k,l,m] obtained by concatenating all the linear combinations described above. So it has 210 linear combinations, and each one has 210 coefficients. solutionmatrix is the 210 x 210 matrix of the coefficients in alpha. The only output from the program is the determinant of this matrix. *********************************************************************) M=Table[ Sum[ a[Min[i,j],Max[i,j],k,l,m] * x^k * y^l * z^m ,{k,0,4},{l,0,4-k},{m,0,4-k-l}] ,{i,1,3},{j,1,3}]; vars=Flatten[Table[ a[i,j,k,l,m], {i,1,3},{j,i,3},{k,0,4},{l,0,4-k},{m,0,4-k-l}]]; flatM={M[[1,1]],M[[1,2]],M[[1,3]],M[[2,2]],M[[2,3]],M[[3,3]]}; (* general functions *) Coef[function_,var_,power_]:= If[power==0,function/.var->0, Coefficient[function,var^power] ]; (* for cubic quadratic divergence *) divM={D[M[[1,1]],x]+D[M[[1,2]],y]+D[M[[1,3]],z], D[M[[2,1]],x]+D[M[[2,2]],y]+D[M[[2,3]],z], D[M[[3,1]],x]+D[M[[3,2]],y]+D[M[[3,3]],z]}; CubicsQuadratics= Module[{vars}, vars=Flatten[Table[{i,j,k},{i,0,3},{j,0,3-i},{k,0,3-i-j}] ,2]; Complement[vars, Flatten[Table[{i,j,k},{i,0,1},{j,0,1-i},{k,0,1-i-j}] ,2] ] ] NonLinearCoeffs[f_]:= Map[Coef[ Coef[ Coef[f,x,#[[1]] ],y,#[[2]] ], z,#[[3]]]&,CubicsQuadratics]; quotientdivlinear=Flatten[Map[NonLinearCoeffs,divM]]; (* for tetrahedron *) OldSimplexIntegrate[monomial_]:=MonSimplexIntegrate[monomial]= Integrate[monomial,{x,0,1},{y,0,1-x},{z,0,1-x-y}]; SimplexIntegrate[f_]:= Sum[ monomial=x^i * y^j * z^k; Coef[Coef[Coef[f,x,i],y,j],z,k]* OldSimplexIntegrate[monomial], {i,0,4},{j,0,4-i},{k,0,4-i-j}]; alphatetrahedron=SimplexIntegrate[flatM]; (* for faces *) OldstIntegrate[monomial_]:=OldstIntegrate[monomial]= Integrate[monomial,{s,0,1},{t,0,1-s}]; stIntegrate[f_]:= Sum[ monomial=s^i*t^j; Coef[Coef[f,s,i],t,j]* OldstIntegrate[monomial] ,{i,0,5},{j,0,5-i}] n[i_]:=If[i==0, {1,1,1}, IdentityMatrix[3][[i]] ]; alphaface1= Join[ stIntegrate[((M/.{x->0,y->s,z->t}).n[1])], stIntegrate[(s*((M/.{x->0,y->s,z->t})).n[1])], stIntegrate[(t*((M/.{x->0,y->s,z->t})).n[1])] ]; alphaface2= Join[ stIntegrate[((M/.{x->s,y->0,z->t}).n[2])], stIntegrate[(s*((M/.{x->s,y->0,z->t})).n[2])], stIntegrate[(t*((M/.{x->s,y->0,z->t})).n[2])] ]; alphaface3= Join[ stIntegrate[((M/.{x->s,y->t,z->0}).n[3])], stIntegrate[(s*((M/.{x->s,y->t,z->0})).n[3])], stIntegrate[(t*((M/.{x->s,y->t,z->0})).n[3])] ]; alphaface0= Join[ stIntegrate[((M/.{x->s,y->t,z->1-s-t}).n[0])], stIntegrate[(s*((M/.{x->s,y->t,z->1-s-t})).n[0])], stIntegrate[(t*((M/.{x->s,y->t,z->1-s-t})).n[0])] ]; alphafaces=Join[alphaface0,alphaface1,alphaface2,alphaface3]; (* for edges *) OldsIntegrate[monomial_]:=OldsIntegrate[monomial]= Integrate[monomial,{s,0,1}] sIntegrate[f_]:= Sum[ monomial=s^i; Coef[f,s,i]* OldsIntegrate[monomial] ,{i,0,6}] vertex[1]={1,0,0}; vertex[2]={0,1,0}; vertex[3]={0,0,1}; vertex[0]={0,0,0}; alphaedge[i_,j_]:= Module[{p,q,v,w,subst,k,l}, p=n[i]; q=n[j]; complement=Complement[{0,1,2,3},{i,j}]; k=complement[[1]]; l=complement[[2]]; v=vertex[k]; w=vertex[l]; subst=s*v+(1-s)*w; { sIntegrate[ p.(M/.{x->subst[[1]],y->subst[[2]],z->subst[[3]]}).p], sIntegrate[ q.(M/.{x->subst[[1]],y->subst[[2]],z->subst[[3]]}).p], sIntegrate[ (v-w).(M/.{x->subst[[1]],y->subst[[2]],z->subst[[3]]}).p], sIntegrate[ q.(M/.{x->subst[[1]],y->subst[[2]],z->subst[[3]]}).q], sIntegrate[ q.(M/.{x->subst[[1]],y->subst[[2]],z->subst[[3]]}).(v-w)], sIntegrate[ (s*p.(M/.{x->subst[[1]],y->subst[[2]],z->subst[[3]]}).p)], sIntegrate[ (s*q.(M/.{x->subst[[1]],y->subst[[2]],z->subst[[3]]}).p)], sIntegrate[ (s*(v-w).(M/.{x->subst[[1]],y->subst[[2]],z->subst[[3]]}).p)], sIntegrate[ (s*q.(M/.{x->subst[[1]],y->subst[[2]],z->subst[[3]]}).q)], sIntegrate[ (s*q.(M/.{x->subst[[1]],y->subst[[2]],z->subst[[3]]}).(v-w))], sIntegrate[ (s^2*p.(M/.{x->subst[[1]],y->subst[[2]],z->subst[[3]]}).p)], sIntegrate[ (s^2*q.(M/.{x->subst[[1]],y->subst[[2]],z->subst[[3]]}).p)], sIntegrate[ (s^2*(v-w).(M/.{x->subst[[1]],y->subst[[2]],z->subst[[3]]}).p)], sIntegrate[ (s^2*q.(M/.{x->subst[[1]],y->subst[[2]],z->subst[[3]]}).q)], sIntegrate[ (s^2*q.(M/.{x->subst[[1]],y->subst[[2]],z->subst[[3]]}).(v-w))], sIntegrate[ ((v-w).(M/.{x->subst[[1]],y->subst[[2]],z->subst[[3]]}).(v-w))] } ] alphaedges=Join[alphaedge[0,1],alphaedge[0,2],alphaedge[0,3], alphaedge[1,2],alphaedge[1,3],alphaedge[2,3]]; (* for vertices *) alphavertices=Join[flatM/.{x->0,y->0,z->0}, flatM/.{x->1,y->0,z->0}, flatM/.{x->0,y->1,z->0}, flatM/.{x->0,y->0,z->1}]; (* Put them together, form a matrix and calculate its determinant *) alpha=Join[alphavertices,alphaedges,alphafaces,alphatetrahedron, quotientdivlinear]; solutionmatrix= Table[ Coefficient[alpha[[i]],vars[[j]]] ,{i,1,Length[alpha]},{j,1,210}]; Det[solutionmatrix]