I've been looking at animating the motion of some (rather nonstandard) robots: this involves interpreting the kinematic constraints as a bunch of simultaneous polynomial equations in several (7) variables with nullity 1, then tracing the motion path via continuous solution of the equations as a function of the time parameter. If possible, I want to solve these in real time --- i.e. fast enough to keep pace with the graphical display (say 30 frames per second). The classical Adams-Bashforth-Moulton (predictor-corrector) method for numerical solution of a differential equation seems a suitable vehicle for the task: the third order version costs only 2 iterations for acceptable accuracy. [For some reason I had thought first of Runge-Kutta-Merson, which however requires the differentials to be available as an explicit function of time.] To cast the motion equations as an ODE system, they must be differentiated and the resulting linear equations for the derivatives solved --- in practice numerically, at each step. Now, we already have to hand an approximate solution to these linear equations --- the derivatives at the previous step --- so an iterative algorithm looks worth trying. The textbook method is Gauss-Seidel, which (with luck, and a diagonally dominant matrix) converges geometrically. I have to admit I was already prejudiced; and after spending an evening establishing how cumbersome and ineffective this idea proves to be [at any rate, in this application], in desperation I cooked up a home-grown alternative. It's a matrix variation of an old trick assembly-coders used in the days before hardware division: starting from an approximation b to 1/a, iterate b := 2 b - a b^2. In the matrix context this becomes simply B := 2 B - B A B. It's easy to show that if I - A B = D, then for the next iterate B' we have I - A B' = D^2; so provided the spectral radius (maximum eigenvalue modulus) of D is d < 1, the algorithm converges and is second-order. If A is rectangular or singular, it converges instead to the Moore-Penrose pseudo-inverse. This is dead simple to program, and fast --- 2 iterations suffice for the application. Admittedly it is still slightly slower (maybe 25%) than solving the equations directly; though if matrix multiplication were hardwired --- as already in some graphics accelerators --- it might well be faster. Yet I don't recall ever seeing it mentioned in print --- can anybody out there enlighten me? To use it for matrix inversion in general, we'd need to construct an initial approximation B for which d < 1. This unfortunately looks nontrivial. [In the original real number context, it was worked around using the fact that binary representation gives an easy approx. base 2 log; and for reduced a ~ 1 we can just choose b = 1.] Fred Lunnon