Linear Systems



I don't know if it's elegant, but you might look at linsolve_by_lu.

(%i1) m : matrix([1,2],[3,4])$
(%i2) b : [3,7]$

(%i3) linsolve_by_lu(m,b);
(%o3) [matrix([1],[1]),false]

The first list member is the solution; the second list member is
either 'false' or a lower bound to the matrix condition number.
When the calculation uses non-floating point numbers, the second
argument should always be false.

To solve m x = m^2 for x, do

  (%i4) linsolve_by_lu(m, m.m);
  (%o4) [matrix([1,2],[3,4]),false]

To solve using double floats, do

  (%i5) linsolve_by_lu(m,b, 'floatfield);
  (%o5) [matrix([1.0],[1.0]),21.77777777777778]

Check the condition number

  (%i6) mat_cond(m,'inf);
  (%o6) 21

Solve using CRE expressions:

  (%i7) m : genmatrix(lambda([i,j],1/(i+j+a)),4,4)$

  (%i8) b : makelist(1,i,1,4)$

  (%i9) linsolve_by_lu(m,b, 'crering);
  (%o9) [matrix([-(a^4+14*a^3+71*a^2+154*a+120)/6],..]

Check the solution

  (%i10) m . first(%);
  (%o10) matrix([1],[1],[1],[1])

BW

-----maxima-bounces at math.utexas.edu wrote: -----

>To: maxima at math.utexas.edu
>From: Peter Danenberg <pcd at wikitex.org>
>Sent by: maxima-bounces at math.utexas.edu
>Date: 10/27/2007 03:32AM
>Subject: Linear Systems
>
>Given a coefficient matrix and a right-hand-side vector, this function
>will solve for x_1, ..., x_n:
>
>    msolve(A, b) := block([n, vars, x, eqns],
>      n: length(transpose(A)),
>      vars: makelist(concat(x, i), i, 1, n),
>      x: covect(append(vars, [1])),
>      eqns: flatten(args(echelon(addcol(A, b)) . x)),
>      solve(eqns, vars));
>
>For instance:
>
>    (%i18) A: matrix([1,1,-2], [3,-1,-7], [1,1,1], [2,2,-4])$
>    (%i19) b: covect([-3,2,0,-6])$
>    (%i20) msolve(A, b);
>    (%o20) [[x1 = - 2, x2 = 3, x3 = - 1]]
>
>Problem is, it seems a little ad-hoc; does Maxima have a more elegant,
>native facility?
>_______________________________________________
>Maxima mailing list
>Maxima at math.utexas.edu
>http://www.math.utexas.edu/mailman/listinfo/maxima