Tag: matrix equation

*Upadate:  I’ve implemented a Maxima matrix-vector equation solver with simpler Maxima-specific algorithm in a later post.  That method is based on the built-in function linsolve().  In that post I show an example that breaks linsolve() but that is handled correctly by the backsolve() method.

Is there really not a solver in Maxima that takes matrix A and vector b and returns the solution of ?  Of course we could do invert(A).b, but that ignores consistent systems where isn’t invertible…or even isn’t square.

Here’s a little function matsolve(A,b)  that solves for general  using the built-in Gaussian Elimination routine echelon(), with the addition of a homemade backsolve() function.  The function in turn relies on a little pivot column detector pivot() and my matrix dimension utility matsize(). This should include the possibilities of non-square , non-invertible , and treats the case of non-unique solutions in a more or less systematic way.

matsolve(A,b):=block(
   [AugU],
   AugU:echelon(addcol(A,b)),
   backsolve(AugU)
 );

backsolve(augU):=block(
 [i,j,m,n,b,x,klist,k,np,nosoln:false],
 [m,n]:matsize(augU),
 b:col(augU,n),
 klist:makelist(concat('%k,i),i,1,n-1),
 k:0,
 x:transpose(matrix(klist)),
 for i:m thru 1 step -1 do (
   np:pivot(row(augU,i)), 
   if is(equal(np,n)) then
     (nosoln:true,return())
   else if not(is(equal(np,0))) then
     (x[np]:b[i],
     for j:np+1 thru n-1 do
       x[np]:x[np]-augU[i,j]*x[j])
    ),
 if nosoln then 
    return([])
 else 
   return(expand(x)) 
)$

matsize(A):=[length(A),length(transpose(A))]$

pivot(rr):=block([i,rlen],
 p:0,
 rlen:length(transpose(rr)),
 for i:1 thru rlen do( 
 if is(equal(part(rr,1,i),1)) then (p:i,return())), 
 return(p)
)$