Post Reply 
Newton's Method
01-17-2014, 09:45 PM
Post: #1
Newton's Method
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.

Graph 3D | QPI | SolveSys
Find all posts by this user
Quote this message in a reply
03-29-2014, 06:29 PM
Post: #2
RE: Newton's Method
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.
Find all posts by this user
Quote this message in a reply
03-30-2014, 01:47 AM
Post: #3
RE: Newton's Method
(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?

Graph 3D | QPI | SolveSys
Find all posts by this user
Quote this message in a reply
03-30-2014, 02:45 AM
Post: #4
RE: Newton's Method
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)
Find all posts by this user
Quote this message in a reply
03-30-2014, 02:08 PM
Post: #5
RE: Newton's Method
Did you mean one of the other programs? NEWT() takes its formula from F1.

Graph 3D | QPI | SolveSys
Find all posts by this user
Quote this message in a reply
03-30-2014, 06:25 PM
Post: #6
RE: Newton's Method
Sorry, yes it is the NEW2() program from your "HP Prime - Programming & Newton's Method (Part 2)" video.
Find all posts by this user
Quote this message in a reply
04-07-2014, 07:47 PM
Post: #7
RE: Newton's Method
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().

Graph 3D | QPI | SolveSys
Find all posts by this user
Quote this message in a reply
01-10-2015, 09:52 PM (This post was last modified: 01-10-2015 09:57 PM by Snorre.)
Post: #8
RE: Newton's Method
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 (mn), 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
Find all posts by this user
Quote this message in a reply
01-10-2015, 11:01 PM (This post was last modified: 01-10-2015 11:23 PM by Han.)
Post: #9
RE: Newton's Method
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

Graph 3D | QPI | SolveSys
Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 1 Guest(s)