for my master thesis in physics, I am currently working on the following problem:
Given four symmetric $n\times n$ matrices $A_1,A_2,B_1,B_2$ yielding the commutation relations
\begin{equation}\label{1}
[A_i,B_j]=0
\end{equation}
I get an algebraic subset $X\subset\mathbb{R}^{2n(n+1)}$, i.e., the space spanned by the components of $A_i,B_j$. Note that there are no restrictions on the commutators $[A_1,A_2],[B_1,B_2]$. I want to obtain the irreducible subsets of $X$ and therefore set up the scenario in Macaulay2:
reVars = (x,n) -> deepSplice for i from 1 to n list x_(i,i)..x_(i,n)
generateRing = (Field,x,xx,z,zz,n) -> Field[reVars(x,n),reVars(xx,n),reVars(z,n),reVars(zz,n)]
commutator = (A,B) -> AB - BA
rePart = (R,x,n) -> genericSymmetricMatrix(R,x_(1,1),n)
reCom = (R,x,z,n) -> trim ideal unique flatten entries commutator(rePart(R,x,n),rePart(R,z,n))
completeCom = (R,x,xx,z,zz,n) -> reCom(R,x,z,n) + reCom(R,x,zz,n) + reCom(R,xx,z,n) + reCom(R,xx,zz,n)
kk = ZZ/101
R = generateRing(kk,x,xx,z,zz,2)
I = completeCom(R,x,xx,z,zz,2)
Think of the x_(k,l),xx_(k,l) and z_(k,l),zz(k,l) as components of the $A_i$ and $B_j$, respectively. In the end, by skew symmetry of the commutator, I get an ideal generated by $2n(n-1)$ quadrics.
Now my problem: In case of $n=2$, if I use the command
primaryDecomposition I
everything works out well and I get a result in the blink of an eye. But already when in three dimensions, the above command seems to need an infinite amount of time. I am pretty sure my knowledge of commutative algebra and Macaulay2 is very basic and I hope anyone here has more experience and can help me out on this.
How do I get the primary decomposition in an acceptable amount of time? Or is it daily business that you run this code for hours up to one day?