Matrices.wxm

/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 22.04.0 ] */
/* [wxMaxima: input   start ] */
/* Date: Mon Jun 4 07:39:27 WEST 2001 */
/* Contributor: William Schelter  */
/* Description: List of matrix entries */
ListMatrixEntries(m):=block(
    [ans:[]],
    for v in m do ans:append(ans,v),
    ans
)$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
/* Date: Mon Jun 4 15:56:51 WEST 2001 */
/* Contributor: Richard Fateman  */
/* Description: List of matrix entries */
/* Note: If used this function should be renamed (it is named the same as the previous one) */
ListMatrixEntries (M):=  apply (append, substinpart("[",M,0))$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
/* Date: Thu Feb 27 22:26:40 WET 2003 */
/* Contributor: Barton Willis  */
/* Description: returns the size of a matrix */
size(m) := if matrixp(m) then [length(m),length(first(m))] else false$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
/* Date: Fri Mar 23 14:58:27 WET 2007 */
/* Contributor: Jaime E. Villate */
/* Description: Implementation of cross product */
cross (v1, v2) := block(
    [M],
    if length(transpose(v1)[1])=1 then v1: transpose(v1),
    if length(transpose(v2)[1])=1 then v2: transpose(v2),
    M: addcol(v1,v2),
    makelist((-1)^(i-1)*determinant(submatrix(i,M)),i,1,3)
)$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
/* Date: Sun Aug 10 17:34:13 WEST 2008 */
/* Contributor: Volker van Nek  */
/* Description: Code for gaussian elimination */
gaussian_elimination (A) := block(
    [i, j, m, n, maxi, listarith:true],
    m:length(args(A)),
    n:length(A[1])-1,
    i:1,
    j:1,
    while (i<=m and j<=n) do (
        maxi:i,
        for k:i+1 thru m do
         if abs(A[k,j]) > abs(A[maxi,j]) then maxi:k,
        if A[maxi,j]#0 then (
            [A[i],A[maxi]]: [A[maxi],A[i]],
            A[i]: A[i]/A[i,j],
            for u:i+1 thru m do A[u]: A[u] - A[u,j] * A[i],
            i: i+1 
        ),
        j: j+1 ),
    A )$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
/* Date: Mon Aug 11 09:25:04 WEST 2008 */
/* Contributor: Angelique Sta maria  */
/* Description: Code for gaussian elimination */
gaussian_elimination (myMatrix) := block(
    [A, i, j, n, m, maxi, listarith:true],
    A: copymatrix(myMatrix),
    n:length(args(A)),
    m:length(A[1])-1,
    i:1,
    j:1,
    while (i<=m and j<=n) do (
        maxi:i,
        for k:i+1 thru m do if abs(A[k,j]) > abs(A[maxi,j]) then maxi:k,
        if A[maxi,j]#0 then (
            [A[i],A[maxi]]: [A[maxi],A[i]],
            A[i]: float(A[i]/A[i,j]),
            for u:i+1 thru m do A[u]: A[u] - A[u,j] * A[i],
            i: i+1 ),
        j: j+1 ),
    A)$    
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
/* Date: Wed Nov 18 11:20:39 WET 2009 */
/* Contributor: Stefano Ferri */
/* Description: Tries to handle submatrices */
subm ( _m, _rows, _cols ) := block(
    [i,j,_x],
    if not matrixp(_m) then error("First input must be a matrix"),
    if not listp(_rows) and not integerp(_rows) then error("Second input must be a list of rows or a single integer"),
    if not listp(_cols) and not integerp(_cols) then error("Third input must be a list of columns or a single integer"),
    if not listp(_rows) and integerp(_rows) then _rows : [_rows],
    if not listp(_cols) and integerp(_cols) then _cols : [_cols],
    for i in _rows do if not integerp(i) or i<=0 then error("Elements in rows list must be integers > 0"),
    for i in _cols do if not integerp(i) or i<=0 then error("Elements in columns list must be integers > 0"),
    for i in _rows do if i>matrix_size(_m)[1] then error("Row indices must be less or equal than the total number of rows"),
    for i in _cols do if i>matrix_size(_m)[2] then error("Column indices must be less or equal than the total number of columns"),
    _rows : sort(unique(_rows)),
    _cols : sort(unique(_cols)),
    return(genmatrix(lambda([i,j],_m[_rows[i],_cols[j]]),length(_rows),length(_cols)))
)$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
/* Date: Sat Oct 27 11:21:07 WEST 2007 */
/* Contributor: Peter Danenberg */
/* Description: Insert the matrix B after the i-th row of A */
insertrows(A, B, i) := block(
    [n],
    head(A, n) :=if n > 0 then head(submatrix(length(A), A), n - 1) else A,
    tail(A, n) :=if n > 0 then tail(submatrix(1, A), n - 1) else A,
    n: length(A),
    addrow(head(A, n - i), B, tail(A, i))
);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
/* Date: Thu Jan 29 16:09:08 WET 2004 */
/* Contributor: Andrei Zorine */
/* Description: nondestructively deletes the i-th row from the matrix m, returns a new matrix */
delrow(m,i):=block(
    [m1,j],
    if i>1 then (
        m1:row(m,1),
        for j:2 thru i-1 do m1:addrow(m1,row(m,j))
    ),
    if i<length(m) then (
        if symbolp(m1) then (
            m1:row(m,i+1),
            i:i+1),
        for j:i+1 thru length(m) do m1:addrow(m1,row(m,j))
    ),
    m1
)$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
/* Date: Thu Jan 29 16:09:08 WET 2004 */
/* Contributor: Andrei Zorine */
/* Description:  nondestructively deletes the i-th column from the matrix m, returns a new matrix */
delcol(m,i):=block(
    [m1,j,sz],
    sz:length(m[1]),
    if i>1 then (
        m1:col(m,1),
        for j:2 thru i-1 do m1:addcol(m1,col(m,j))
    ),
    if i<sz then (
        if symbolp(m1) then (
            m1:col(m,i+1),
            i:i+1
        ),
        for j:i+1 thru sz do m1:addcol(m1,col(m,j))
    ),
    m1
)$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
/* Date: Thu Jan 29 16:59:47 WET 2004 */
/* Contributor: Stavros Macrakis */
/* Description:  Deletes the row r */
removerow(m,r):=genmatrix(lambda([i,j],m[ if i>=r then i+1 else i, j ]),length(m)-1,length(m[1]))$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
/* Date: Thu Jan 29 16:59:47 WET 2004 */
/* Contributor: Stavros Macrakis */
/* Description:  Deletes the column c */
removecol(m,c):=genmatrix(lambda([i,j],m[ i, if j>=c then j+1 else j ]),length(m),length(m[1])-1)$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
/* Date: Thu Jan 29 18:47:17 WET 2004 */
/* Contributor: Stavros Macrakis */
/* Description:  Deletes the row r */
/* Note:If the previous functions are used this must be renamed */
removerow(m,r):=part(m,allbut(r))$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
/* Date: Thu Jan 29 18:47:17 WET 2004 */
/* Contributor: Stavros Macrakis */
/* Description:  Deletes the column c */
/* Note:If the previous functions are used this must be renamed */
removecol(m,c):=transpose(removerow(transpose(m),c))$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
/* Date: Thu Apr 7 12:31:22 WEST 2005 */
/* Contributor: Stavros Macrakis */
/* Description: Returns a row vector with the same number of entries */
flatten_matrix(m):=matrix(apply(append, args(transpose(m))))$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
/* Date: Tue May 13 15:28:50 WEST 2008 */
/* Contributor: Stavros Macrakis */
/* Description: Returns a matrix with the same number of entries but different size (numbers of rows and columns) */
resize_matrix(m,rows,cols) :=block(
    [list: xreduce(append,args(m))],
    genmatrix( lambda([i,j],list[j+(i-1)*cols]), rows, cols)
)$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
/* Date: Sat Dec 21 06:35:10 WET 2013 */
/* Contributor: Mike Valenzuela */
/* Description: Tries to return the power of a matrix */
/* Known issues: Yes */
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 ) ) )
)$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
/* Date: Sat Aug 28 20:42:56 WEST 2004 */
/* Contributor: Larry Prevett */
/* Description: an implementation of Kronecker product */
/* Note: it is better to use block (that is suggested by the parser) */
directprod(a,b):=(
    [p],
    if not (matrixp(a) and matrixp(b)) then error("Invalid arguments to directprod()") 
    else (
        p:zeromatrix(length(a)*length(b),length(a[1])*length(b[1])),
        for i thru length(a) do for j thru length(a[1]) do for k thru length(b) do for l thru length(b[1]) do p[k+length(b)*(i-1),l+length(b[1])*(j-1) ]:a[i,j]*b[k,l],
        p
    )
)$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
/* Date: Tue Jul 13 20:30:58 WEST 2004 */
/* Contributor: Barton Willis */
/* Description: an implementation of random matrix */
randommatrix(n) := apply(matrix,makelist(makelist(?random(100),k,1,n),j,1,n))$
/* [wxMaxima: input   end   ] */



/* Old versions of Maxima abort on loading files that end in a comment. */
"Created with wxMaxima 22.04.0"$