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

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
>>
>
>