Newton's Method

[Han]

Newton's Method computes the root of a function \( f(x) \) using linear approximations of \( f(x) \) via tangent lines. The formula is
\[ x_{new} = x_{old} - \frac{f(x_{old})}{f'(x_{old})} \]
We use F1 to store the formula for \( f(X) \) (note the capital letter for \( X \) ), and F0 to set up a formula for the right hand side.

CAS computations currently do not recognize local variables. And global variables are generally evaluated during parsing in a program, since the programs run in the Home view. This causes some issues with passing values and expressions to CAS commands. In the Home view, a symbolic expression is surrounded by single-quotes. So in the first example, we create a symbolic expression and "insert" other symbolic expressions (such as the derivative of a function) by using strings. In the later examples, we pass a single string to the CAS command, which parses the string as if typed from the CAS command line. Since the CAS recognizes functional algebra, we create a formula for F0 by simply operating on the functions themselves (the term id is the identity function \( id(x)=x \)) and completely ignore the dummy variables. Since our desired formula is
\[ F0(x):= x - \frac{F1(x)}{F1'(x)} \]
then the equivalent function using functional algebra is
\[ F0 := id - \frac{F1}{F1'} \]

Style 1: Initial attempt as a solution to a specific problem. The function is pre-stored into F1 in the Function App. The user then runs the command NEWT() from the command line, or from the program catalog, or from the [Toolbox] menu (press the "User" menu option).

Code:


export NEWT()
begin
  local n,xold,xnew,err;

  err:=.000001;
  n:=0;
  xnew:=2;
  xold:=xnew-2*err;
  F0:=expr("'X-F1(X)/(" + diff(F1(X),X) + ")'");

  L1:={};

  while (abs(xnew-xold)>err and n<100) do
    n:=n+1;
    L1(n):=xnew;
    xold:=xnew;
    xnew:=F0(xold);
  end;

  L1(n+1):=xnew;

end;

Syle 2: Creating a user interface The INPUT() command is used to allow the user to enter their formula, rather than having it already pre-stored in F1.

Code:


export NEWT2()
begin
  local n,xold,xnew,err,N,f;

  N:=100; err:=.00001; xnew:=1;

  if input(
    {f,xnew,err,N},
    "Newton's Method",
    {"f(X)=", "Guess=", "Error=", "Max Iter.="},
    {
      "Enter the function surrounded by single quotes",
      "Enter the initial guess",
      "Enter the tolerance",
      "Enter the maximum number of iterations"
    },
    {f,xnew,err,N}
  ) then
    F1:=f;

    CAS("F0:=id-F1/F1'");
    L1:={}; L1(1):=xnew;
    for n from 2 to N+1 do
      xold:=xnew;
      xnew:=F0(xold);
      L1(n):=xnew;
      if abs(xnew-xold)<err then break; end;
    end;
    editlist(L1);

  end;

end;

Style 3: Function-like command This style uses functional notation to send inputs to the command NEWT3(). Due to how variables are initialized, this command must be run from the command line. Otherwise, the built-in input screen will prompt for the arguments, and only accept real-valued inputs, and an algebraic expression cannot be entered for f. Usage: NEWT3("X^2-5", 2.1, .0001, 100);

Code:


export NEWT3(f,guess,tol,maxiter)
begin
  local n,xold,xnew,err,N;

  N:=maxiter;
  err:=tol;
  xnew:=guess;

    F1:=f;

    CAS("F0:=id-F1/F1'");
    L1:={}; L1(1):=xnew;
    for n from 2 to N+1 do
      xold:=xnew;
      xnew:=F0(xold);
      L1(n):=xnew;
      if abs(xnew-xold)<err then break; end;
    end;
    editlist(L1);

end;

Style 4: Mix of Style 1 and 3 Usage: Place formula into F1 in the Function App, and run the program NEWT3 either from the command line, the program catalog, or the [Toolbox] interface.

Code:


export NEWT3(guess,tol,maxiter)
begin
  local n,xold,xnew,err,N;

  N:=maxiter;
  err:=tol;
  xnew:=guess;

    CAS("F0:=id-F1/F1'");
    L1:={}; L1(1):=xnew;
    for n from 2 to N+1 do
      xold:=xnew;
      xnew:=F0(xold);
      L1(n):=xnew;
      if abs(xnew-xold)<err then break; end;
    end;
    editlist(L1);

end;

REMARKS:

A more user-friendly program would actually preserve the existing formula for F0, if it exists, and restore it so that the user does not lose their previously stored formulas after running our program. To do this, one could implement:

Code:


LOCAL oldfunc=""; // ensure string type
IFFERR
  oldfunc:=STRING(F0);
THEN
  oldfunc:="";
END;

The error trap is to account for the possibility that there may be no pre-existing F0 formula. We restore the user's old F0 formula by simply using:

Code:


F0:=oldfunc;

Fortunately, this works even if oldfunc is an empty string.

[FrankiD]

I have left a few messages on the NEWT program on youtube if you want to catch up on some of the details there.

I am noticing that the NEWT program behaves differently depending on the Entry type selected in the Home Settings. The Entry method selected affects whether I have to use 2 sets of single quotes or only 1 set.

[Han]

(03-29-2014 06:29 PM)FrankiD Wrote:  I have left a few messages on the NEWT program on youtube if you want to catch up on some of the details there.

I am noticing that the NEWT program behaves differently depending on the Entry type selected in the Home Settings. The Entry method selected affects whether I have to use 2 sets of single quotes or only 1 set.

I wasn't able to find any comments on Youtube that seem related to this post. Could you provide a specific example of what you have to do in each entry mode?

[FrankiD]

My calculator (HP Prime) was rebooting when I would run the NEWT program from your 2nd video. Once I formatted the calculator (something I found in the hpmuseum forum), I guess this wiped the original firmware and fixed the reboot (for the most part...it's happened one time after the formatting).

Since then, I have been able to figure out that the NEWT program works differently depending on the Entry method.

In settings (Shift-Home), the Entry method can be Textbook, Algebraic or RPN

RPN Entry method:
Must enter the function as ''X^2-31'' (double quotes or 2 sets of single quotes)
Textbook or Algebraic Entry method:
Can enter 'X^2-31' (single quotes)

[Han]

Did you mean one of the other programs? NEWT() takes its formula from F1.

[FrankiD]

Sorry, yes it is the NEW2() program from your "HP Prime - Programming & Newton's Method (Part 2)" video.

[Han]

I see what you mean now. The issue is with how INPUT() behaves. Right now, it is quite particular in how it works due to parser used for the Home view, which of course depends on your settings. I imagine that future updates (esp. if they adjust the parser) will probably add in new quirks to commands like INPUT().

[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

[Han]

Snorre,

I think your program is fine, but it doesn't match exactly with the math above your code. The math presumes f(x) is a column vector whereas it appears you have a row vector. Strangely the CAS seems to not mind this and allows the multiplication to go through as is. Nice find!

I was planning to write something similar (to incorporate into the equation library program), but basically had a few extra lines in mind (that now seem unnecessary). I was thinking it was necessary to transpose the f vector to ensure valid multiplication. The CAS seems to be rather flexible in this respect!

Lastly, defining the jacobi as you did vs leaving it untransposed (how I thought about doing it) comes down to deciding between two transposes vs a single transpose (needed prior to taking inverse). I think the latter choice would give a slight speedup in cases where the number of iterations is significantly large.

As for experience, I doubt I have any more experience than anyone here -- perhaps more chunks of free time to tinker, maybe.

Edit: I think there may be a few other ways to exit the iterations. You've listed a few, but we could also exit based on whether max|f(x_i)| is within a certain tolerance as suggested by: http://ocw.usu.edu/Civil_and_Environment...Matlab.pdf

[Gene222]

Regarding Style 2 INPUT example in the first post, the use of the INPUT command will not allow an expression to be inputted. The solution is discussed by Han on a different forum, which states:

"Han
01-17-2017, 04:14 PM
INPUT() has gone through a number of changes since that post. There are two ways to get around the issue you have:

1) Initialize your f and f1 to objects of the type you wish to use in the INPUT() command. When you do not initialize them, then they are set to the value 0 by default (type: real number). Thus, when a user then enters something that is not a real number, the INPUT() command complains of invalid input. You can, for example, set f:='X^2-5' and f1:='2*X' before the INPUT() command, and then it will expect an algebraic expression.

2) The other method use the more advanced form of INPUT() to enforce the input type for each field (i.e. f and f1 and be required to expect an expression). The help screen for INPUT() explains everything you need to make use of the advanced form."