algorithms for 'invert', was: program works with maxima 5.29.1 but freezes with maxima 5.30.0
Subject: algorithms for 'invert', was: program works with maxima 5.29.1 but freezes with maxima 5.30.0
From: Barton Willis
Date: Sat, 20 Apr 2013 11:46:55 +0000
> For example, if the matrix is 'almost' singular (i.e., det(M) almost zero), then things can go haywire.
A matrix can be ill conditioned with a determinant of one; famous example:
(%i1) m : matrix([63245986, 102334155], [102334155, 165580141])$
The determinant of m is 1, but its condition number is about 10^16
(%i2) [determinant(m), mat_norm(m,'inf) * mat_norm(m^^-1,'inf)];
(%o2) [1,71778070001175616]
Let m be a coefficient matrix
(%i4) eq : m . matrix([x],[y])- matrix([1],[1])$
(%i5) eq : xreduce('append,args(eq))$
Solve using rational numbers
(%i6) linsolve(eq,[x,y]);
(%o6) [x=63245986,y=-39088169]
Solve using binary64
(%i7) linsolve(float(eq),[x,y]),keepfloat : true;
(%o7) [x=3.1622993000000007*10^7,y=-1.95440845*10^7]
Partial pivoting doesn't save the calculation--at least the estimated condition number will warn the alert user:
(%i8) linsolve_by_lu(m,[1,1],'floatfield);
(%o8) [matrix([4.1475558898394227*10^7],[-2.5633305101605773*10^7]),4.7070743271309312*10^16]
--Barton