-----maxima-bounces at math.utexas.edu wrote: -----
>Is there a special function to calculate the
>exponential of a matrix?
Yes, there are at least two such Maxima packages. I know
the most about the one in linearalgebra; the other package
is /share/contrib/diag.mac. The linearalgebra package
is ill-suited for numerical calcuation (even with small
matrices). A short demo of the linearalgebra package:
(%i1) load(linearalgebra)$
(%i2) m : matrix([1,2],[2,1]);
(%o2) matrix([1,2],[2,1])
Find the spectral represenation ([list of eigenvalues, list of projections,
nilpotent])
(%i3) spectral_rep(m);
(%o3) [[3,-1],[matrix([1/2,1/2],[1/2,1/2]),matrix([1/2,-1/2],[-1
/2,1/2])],matrix([0,0],[0,0])]
check
(%i4) %[1] . %[2] + %[3];
(%o4) matrix([1,2],[2,1])
compute exp(m * x)
(%i5) matrixexp(m,x);
(%o5)
matrix([(%e^(-x)*(%e^(4*x)+1))/2,(%e^(-x)*(%e^(4*x)-1))/2],[(%e^(-x)*(%e^(4*x)-1))/2,(%e^(-x)*(%e^(4*x)+1))/2])
try a matrix with a nonzero nilpotent part
(%i6) m : matrix([1,2],[3,1-%i*sqrt(24)]);
(%o6) matrix([1,2],[3,1-2*sqrt(6)*%i])
(%i7) spectral_rep(m);
(%o7)
[[1-sqrt(6)*%i],[matrix([1,0],[0,1])],matrix([sqrt(6)*%i,2],[3,-sqrt(6)*%i])]
(%i8) %[1] . %[2] + %[3];
(%o8) matrix([1,2],[3,1-2*sqrt(6)*%i])
(%i9) last(%o7)^^2;
(%o9) matrix([0,0],[0,0])
(%i10) matrixexp(m,x);
(%o10)
matrix([(sqrt(6)*%i*x+1)*%e^(x-sqrt(6)*%i*x),2*x*%e^(x-sqrt(6)*%i*x)],[3*x*%e^(x-sqrt(6)*%i*x),-(sqrt(6)*%i*x-1)*%e^(x-sqrt(6)*%i*x)])
matrixexp fails even for easy problems with floating point entries:
(%i24) m : matrix([1.0,2.0],[2.0,1.0]);
(%o24) matrix([1.0,2.0],[2.0,1.0])
(%i25) matrixexp(%,x), keepfloat;
`rat' replaced 1.0 by 1//1 = 1.0
...
Unable to find the spectral representation
Other matrix functions:
(%i31) m : matrix([1,2],[2,1]);
(%o31) matrix([1,2],[2,1])
(%i32) matrixfun(lambda([x], sqrt(x)),m);
(%o32)
matrix([(%i+sqrt(3))/2,-(%i-sqrt(3))/2],[-(%i-sqrt(3))/2,(%i+sqrt(3))/2])
(%i33) expand(%^^2);
(%o33) matrix([1,2],[2,1])
(%i34) matrixfun(lambda([x], x^(1/5)),m);
(%o34) matrix([(3^(1/5)-1)/2,(3^(1/5)+1)/2],[(3^(1/5)+1)/2,(3^(1/5)-1)/2])
(%i35) expand(%^^5);
(%o35) matrix([1,2],[2,1])
Barton