It seems your code is correct, and what a speedup; now the above for-loop takes only a little more than a second to complete, i.e., it is ten times faster.

Thanks for your link to your function Tensorproduct. I have done a primitive comparison of speed among the four versions of the Kronecker product:

with(LinearAlgebra):
A := Matrix(4,4,(i,j) -> a||i||j):
B := Matrix(4,4,(i,j) -> b||i||j):
t := time():
for i from 1 to 1000 do
    kronprod(A,B)		# Robert Israel
#   Tensorproduct(A,B)	# Sandor Szabo
#   KroneckerProd(A,B)	# Douglas Keenan
#   KroneckerProduct(A,B)	# John Fredsted
end do:
time() - t;

which on my old timer computer results in (subsequently choosing the appropiate line) the following times (in order of appearance in the for-loop): 20.8, 36.4, 11.0, and 10.4. Some comments:

• I guess that Robert Israel's version should not really be compared to the other three because his version is more general.
• I guess (please correct me if I am mistaken) that what slows down your version is (even though quite ingenious) the iterative use of a construction like Y := Matrix([Y,X]).

Thanks for that link, from which I found the source file for the Kronecker product at

http://www.math.ubc.ca/~israel/advisor/advisor6/kroneck5.txt

Thanks for your efforts, but a simpler code seems possible:

KroneckerProduct := proc(a::Matrix,b::Matrix)
	local i,j,aRow,aCol,bRow,bCol,p;
	aRow,aCol := LinearAlgebra:-Dimension(a);
	bRow,bCol := LinearAlgebra:-Dimension(b);
	p := Matrix(aRow * bRow,aCol * bCol);
	for i from 1 to aRow do
	for j from 1 to aCol do
		p[
			(i - 1) * bRow + 1..i * bRow,
			(j - 1) * bCol + 1..j * bCol
		] := a[i,j] * b
	end do;
	end do;
	p
end proc:

An example:

with(LinearAlgebra):
KroneckerProduct(Matrix(3,3,(i,j) -> a||i||j),IdentityMatrix(2));

Matrix([[a11, 0, a12, 0, a13, 0], [0, a11, 0, a12, 0, a13], [a21, 0, a22, 0, a23, 0], [0, a21, 0, a22, 0, a23], [a31, 0, a32, 0, a33, 0], [0, a31, 0, a32, 0, a33]])

PS: I do apologize if I have made you believe that the intention of my post was to obtain an explicit code for calculating the Kronecker products of matrices. The sole purpose of my post was to put forward the wish that such a function would be a built-in one in the next release of Maple.

Thanks for that pleasant surprise.

As far as I am aware, in Maple 11 there is no implementation of something completely basic as a commutator or an anticommutator of matrices, compare my previous comment in Only vector fields. Also, there does not seem to be any implementation of the equally basic Kronecker product, i.e., the tensor product of matrices, so that instead of

with(LinearAlgebra):
M := Matrix(4,4):
M[1..2,1..2] := a11*IdentityMatrix(2):
M[1..2,3..4] := a12*IdentityMatrix(2):
M[3..4,1..2] := a21*IdentityMatrix(2):
M[3..4,3..4] := a22*IdentityMatrix(2):

say, one could just write

with(LinearAlgebra):
A := Matrix([[a11,a12],[a21,a22]]):
M := KroneckerProduct(A,IdentityMatrix(2)):

Note added: This issue seems now to have been resolved, see Order of a PDE. Comforting to hear that it surprised you too. Actually, now having read some of the help text [why I did not do that prior to posting is one of the returning conundrums of my life :-)] select does perform correctly:

The select function selects the operands of the expression e which satisfy the Boolean-valued procedure f...

So select seems to be of no help in filtering out the appropiate terms in the differential equations, or what? P.S.: Now I did it again; misplacing a post. This post should of course have been a reply to Doug Meade. I apologize.

Nice to hear! If that entails being able to perform both abstract index manipulations as well as number crunching, both in curved spacetimes, in one and the same package, then surely I will look forward to that. Regarding the distinction between space and time, I do not mean being able to split space and time apart (working covariantly usually is much more powerful). What I mean is, in line with the author of this thread, a clear distinction between covariant and contravariant indices in spacetimes with indefinite metrics of any general signature (p,q) where p+q need not equal 4.

Nice to hear! If that entails being able to perform both abstract index manipulations as well as number crunching, both in curved spacetimes, in one and the same package, then surely I will look forward to that. Regarding the distinction between space and time, I do not mean being able to split space and time apart (working covariantly usually is much more powerful). What I mean is, in line with the author of this thread, a clear distinction between covariant and contravariant indices in spacetimes with indefinite metrics of any general signature (p,q) where p+q need not equal 4.

You are quite right; this thread is becoming quite difficult to follow. Probably that is the reason why I did not properly notice the comment your are referring to. Anyway, thanks for the answer, which reassures me that it is not me having lost the ability to reason; the problem lies in the physics package, not in my brain :-).

You are quite right; this thread is becoming quite difficult to follow. Probably that is the reason why I did not properly notice the comment your are referring to. Anyway, thanks for the answer, which reassures me that it is not me having lost the ability to reason; the problem lies in the physics package, not in my brain :-).

I think it boils down to the following question: Why does the following code return zero?

restart:
with(Physics):
Setup(
   quantumoperators = {a,b},
   algebrarules = {
      %AntiCommutator(a,a) = 0,
      %AntiCommutator(b,b) = 0,
      %AntiCommutator(a,b) = 1}
):
AntiCommutator(a,b);

To me it should certainly return 1, and indeed it does if I uncomment either the %AntiCommutator(a,a) = 0 line or the %AntiCommutator(b,b) = 0 line, or both lines, which just makes it even more confusing, of course.

I think it boils down to the following question: Why does the following code return zero?

restart:
with(Physics):
Setup(
   quantumoperators = {a,b},
   algebrarules = {
      %AntiCommutator(a,a) = 0,
      %AntiCommutator(b,b) = 0,
      %AntiCommutator(a,b) = 1}
):
AntiCommutator(a,b);

To me it should certainly return 1, and indeed it does if I uncomment either the %AntiCommutator(a,a) = 0 line or the %AntiCommutator(b,b) = 0 line, or both lines, which just makes it even more confusing, of course.

The code you ask for was provided by a link in the previous post. However, I am confident that the issues associated with that post will be resolved if those associated with the post The bare essentials are resolved.

The code you ask for was provided by a link in the previous post. However, I am confident that the issues associated with that post will be resolved if those associated with the post The bare essentials are resolved.