One missing maxima feature, that I have the math behind, but not a program

Hi Robert Dodier,

I'm not sure matrixfun(lambda([x], x^t), A) can do what I need it to do. So
I'm going to explain, in way too much detail, what I mean.

---A^^t objective
The matrix A^^t returns a symbolic matrix power. I think it requires A to
be a square matrix, otherwise it doesn't make much sense.

---A^^t method
I'm guess there are multiple ways to do this, but a Hungarian mathematician
taught me this method (since its nearly identical for calculating the
matrix exponential).

(i) First, reexpress the general solution in terms of unknown matrices B1,
B2, ... and eigenvalues:
A^^t = B1 * lambda1^t + B2 * lambda2^t +...

When there are repeated eigenvalues you use the same trick as in solving
differential equations. For example, one eigenvalue with an algebraic
multiplicity of four would have the general form:
A^^t = B1_1 * lambda1^t + B1_2 * t * lambda1^t + B1_3 * t^2 * lambda1^t
+ B1_4 * t^3 * lambda1^t + B2 * lambda2^t + ...

(ii) Then we can generate additional equations by incrementing t.
A^^t = B1_1 * lambda1^t + B1_2 * t * lambda1^t + B1_3 * t^2 * lambda1^t
+ B1_4 * t^3 * lambda1^t + B2 * lambda2^t + ...
A^^(t+1) = B1_1 * lambda1^(t+1) + B1_2 * (t+1) * lambda1^(t+1) + B1_3 *
(t+1)^2 * lambda1^(t+1) + ...
A^^(t+2) = B1_1 * lambda1^(t+2) + B1_2 * (t+2) * lambda1^(t+2) + B1_3 *
(t+2)^2 * lambda1^(t+2) + ...
...
Do this to get enough linearly independent equations (where the variables
are actually the matrices B1_1, B1_2, B1_3, B1_4, B2, B3...)

(iii) If there are no zero eigenvalues:
Since these equations are true regardless the value of t, you can set t=0
and then solve the multilinear equation for each of the B matrices in terms
of Identity, A^^1, A^^2, etc.

If zero is an eigenvalue:
Same as above, but you may have to set t=1 instead of zero to avoid
degenerate solutions, such as A^^t = Identity.

(iv) Then solve using linear algebra. You're augmented matrix will look
something like this (assuming we set t=1)

(Header on the top row)
B1_1,   B1_2,   B1_3 ,  B1_4  | A, A^^2, A^^3, A^^4
----------------------------------------------------
  1 ,     1 ,    1   ,   1    | 1,   0 ,   0 ,   0
 L1 ,  2*L1 ,  4*L1  , 8*L1   | 0,   1 ,   0 ,   0
L1^2, 3*L1^2, 9*L1^2 , 27*L1^2| 0,   0 ,   1 ,   0
L1^3, 4*L1^3, 16*L1^3, 64*L1^3| 0,   0 ,   0 ,   1

Thus using linear algebra techniques, you can solve for each B matrix in
terms of the A matrix and its powers.

---A^^t properties when calculated using the above method
When A is invertible, A^^t, t=-1 is the inverse.
When A is not invertible, A^^t, t=-1 is the Drazin inverse.
When A is not invertible, A^^t, t=0 need not be the identity matrix.
A^^t1 . A^^t2 = A^^(t1 + t2)

---How the above method is related to matrixexp
This same method can calculate matrix exponentials with two changes:
(1) instead of multiplying each B matrix by lambda^t, multiply by
exp(lambda*t)
(2) instead of incrementing t, to get more equations, take the derivative
with respect to t.
The rest of the procedure is the same if I recall correctly at the moment.

---Conclusions
I think I could throw together a hodge-podge of methods to go through this
process, but I'm not familiar enough with Maxima to: make it robust enough
to behave nicely when the matrix is not square; allow it to solve a simpler
set of linear equations when the A matrix does not have a zero eigenvalue;
make an example that shows how to use it; make documentation for Maxima so
people know that A^^-1 is the Drazin inverse and not the Moore-Penrose
inverse.

Also, I think the function 'matrixexp' is currently undocumented.



On Tue, Dec 17, 2013 at 12:54 PM, Robert Dodier <robert.dodier at gmail.com>wrote:

> On 2013-12-07, Mike Valenzuela <mickle.mouse at gmail.com> wrote:
>
> > matrixexp(A * t) is often useful for solving linear differential
> equations.
> > Similarly, A^^t is useful for solving (discrete) difference equations.
>
> Not sure if this is the same as the stuff you are working on, but maybe
> matrixexp and matrixfun in the linearalgebra package are what you are
> looking for. For A^^t, maybe you want matrixfun(lambda([x], x^t), A).
>
> HTH
>
> Robert Dodier
>
> _______________________________________________
> Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima
>