Newton's Method

[Snorre]

Hello Han,

I was thinking about creating my own Newton approximator for the CAS (based on a program I once wrote for my HP42s).
But since you did that already, I'd like to add my thoughts here.

Step 1: Take a vector of functions (expressions) f and a vector of variables x \[ \mathbf{f}(\mathbf{x}) = [f_1(x_1,...,x_n), ..., f_m(x_1,...,x_m)] \] and search for a x* so that f(x*)=0. Thus the iteration becomes \[ \mathbf{d}:=\mathbf{J}^{-1}\mathbf{f}(\mathbf{x}) \qquad \mathbf{x}_\mathrm{new} := \mathbf{x} - \mathbf{d} \] where J is the Jacobian matrix of f(x), i.e. \[ J_{ij} = \frac{\partial f_i(\mathbf{x})}{\partial x_j} \]

Step 2: Since the number of functions may be different from the numbers of variables (m≠n), don't search for an exact solution but for the nearest
\[ \mathbf{d}:=(\mathbf{J}' \mathbf{J})^{-1} \mathbf{J}' \mathbf{f}(\mathbf{x}) \qquad \mathbf{x}_\mathrm{new} := \mathbf{x} - \mathbf{d} \] where J' means transposition of J.

Step 3: Add another threshold (thd) for breaking up iteration, i.e.
1. Stop if iteration count k exceeds maximum: k>kmax
2. Stop if we're near the solution x*: |f(x)|<thf
3. Stop if our current guess doesn't change: |d|<thd

My first naive attempt is

Code:

#cas
newton(exprs,vars,guess,thf,thd,kmax):=
BEGIN
  LOCAL jacobi,X,F,J,D,k;
  jacobi:=transpose(diff(exprs,vars));

  PRINT();
  X:=approx(guess);
  k:=0;
  WHILE (k:=k+1)≤kmax DO
    PRINT("k="+k);

    F:=approx(subst(exprs,vars=X));
    PRINT("|F|="+ABS(F));
    IF ABS(F)<thf THEN BREAK END;

    J:=approx(subst(jacobi,vars=X));
    X:=X-(D:=(trn(J)*J)^-1*trn(J)*F);
    PRINT("|D|="+ABS(D));
    IF ABS(D)<thd THEN BREAK END; 
  END;

  X; 
END;
#end

Usage: newton( [f1(x1,...,xn), ..., fm(x1,...,xn)], [x1,...,xn], [xinit1,...,xinitn], thf, thd, kmax )
Prints out the current values of |f(x)| and |d|. Returns x=[x1,...,xn].

Maybe you want to incorporate the basic approach (I don't even know if it's all correct, but some simple tests where promising) into your program since you're far more experienced with Prime programming (all that fancy interface stuff, error checking, etc.).

Greetings