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

Sorry I'm being spammy.

I made a github project for this:
https://github.com/Mickle-Mouse/Maxima-Matrix-Power




On Fri, Dec 20, 2013 at 11:35 PM, Mike Valenzuela <mickle.mouse at gmail.com>wrote:

> Here is a very crummy basic thing I through together. I know there are
> problems with this code, but I've included it anyways.
>
> /* Michael Valenzuela's original symbolic matrix power */
> /* This very non-robust code is effectively released   */
> /* Under the "New BSD License"                         */
> /* This is compatible with GPL, but is not copy-left   */
> /* I just don't want to be held accountable for this   */
> /* code                                                */
>
> /* Version 1.0 --- not robust, not very efficient version */
> matrix_power(A, t):=block([A_, A_list, B_, B_list, eqs, eigvals, i, j, k,
> numeq, tt, solvedeq, offset],
> /* Initialize Variables */
> /* This is the value we set t to this value */
> offset: 1,
> /* The number of equations */
> numeq: length(A),
> /* The list of relevant variables */
> eqs: makelist(0,i,1,numeq),
> A_list: makelist( A_^^(i+offset)=tt,i,0,numeq-1),
> B_list: makelist( B_[i],i,1,numeq),
>
>
> /* Exact Eigenvalues and calculate powers of A */
> eigvals:(radcan(fullratsimp(eivals(A)))),
> for i:1 thru numeq do block([],
>   A_list[i]: at( A_list[i], tt=A^^(i+offset-1))
> ),
> A_list : reverse(A_list),
>
>
> /* Setup the first equation */
> k:1,
> eqs[1]: 0,
> for i:1 thru length(eigvals[1]) do block([],
>   for j:1 thru eigvals[2][i] do block ([],
>     eqs[1]: eqs[1] + B_list[k] * t^(j-1) * eigvals[1][i]^t,
>     k:k+1
>   )
> ),
> eqs[1]: A_^^t = eqs[1],
>
>
> /* Generate all the necessary equations */
> for i:2 thru numeq do block([],
>   eqs[i]: at(eqs[i-1], [t=1+tt]),
>   eqs[i]: at(eqs[i], [tt=t])
> ),
>
>
> /* Set t=offset and solve the equations */
> solvedeq: linsolve( at(eqs, t=offset), B_list ),
>
>
> /* Substitute the symbolic A_^k with the calculation */
> solvedeq: factorsum(radcan(at( solvedeq, A_list))),
>
>
> /* Substitute the resulting equations into eq0 */
> return( rhs(at(eqs[1], solvedeq ) ) )
> );
>
>
>
>
> On Fri, Dec 20, 2013 at 7:08 PM, Mike Valenzuela <mickle.mouse at gmail.com>wrote:
>
>> 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
>>>
>>
>>
>