OS/2 Shareware BBS: 22 gnu/22-gnu.zip/fweb140x.zip/demos/Newton.web

Text File | 1993-10-30 | 10KB | 289 lines

% Newton.web, fweb, version 1.13, January 17, 1991
% Bart Childs, Tom McCurdy, Clarissa Wilson  [e-mail]
% Texas A&M University, 409-845-5470
% Department of Comp Science, TAMU, College Station, TX 77843-3112
%
%   REVISION HISTORY
%   1.0 January 17, 1990  a companion to Intro.tex which uses this
%       file's line numbers <---====  NOTICE, don't change line numbers
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   LIMBO MATERIAL
%
@x
\font\ninerm=cmr9
\let\mc=\ninerm % medium caps for names like UNIX
\def\WEB{{\tt WEB}}
\def\FWEB{{\tt FWEB}}
\def\Fortran{F{\sc ORTRAN}}
\def\PASCAL{Pascal}
\def\<{$\langle\,$}
\def\>{$\,\rangle$}
\def\C{{\bf C}}
\def\Unix{{\mc UNIX}}  \let\unix=\Unix
\def\today{\ifcase\month\or
   January\or February\or March\or April\or
   May\or June\or July\or August\or
   September\or October\or November\or December\fi
   \space\number\day, \number\year}
\def\title{{The Newton Method as a Short \WEB{} Example}}
\def\topofcontents{\hsize 6in
  \vglue -30pt plus 1fil minus 1.5in
  \centerline{\title}
  \vskip 15pt
  \centerline{\today}
{\bigskip\parskip6pt plus2pt \parindent20pt
The following output has been changed to achieve some efficiency
in paper.  The major modules that appear in the following table of
contents would normally start a new page.  In the interest of
brevity, this has been suppressed in modules 7, 13, and 18.

This is written in John Krommes' \FWEB{} because
\Fortran{} is such a simple language.  It is to show
features of \WEB{}s and is not purported to be the absolutely
best way of coding a solution of this problem.

This is therefore an expository note about the \WEB{} style
of Literate Programming.
\vskip0.5in}
\vfil
\centerline{\bf Table of Contents}
\medskip
\def\?##1]{\hbox to 1in{\hfil##1.\ }}}
\def\botofcontents{\vskip 0pt plus 1fil minus 1.5in}
%
%   END OF LIMBO MATERIAL
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   BEGINNING OF WEB
%
\def\FWEAVE{{\tt FWEAVE}}
\def\sc{\eightrm}

@n

@* The Newton Method.  The Newton method is based on
the Taylor's Series:
$$ f(x) = \sum_{n=0}^{+\infty}{{f^{(n)}(a)}\over {n!}}(x-a)^n =
    f(a) + f'(a)(x-a) + {f''(a)\over{2!}}(x-a)^2 +\ldots+
    {f^{(n)}(a)\over {n!}}(x-a)^n +\ldots $$

  In order to use the Newton method a number of conditions must be met:
$f(x)$ must be analytic and the second derivative of $f(x)$ must be
monotone in the neighborhood of $x$.  If these conditions are met,
then the Newton method will converge quadratically and monotonically
({\it i.e.} with each iteration the number of significant digits will
double).  The Newton method is an iterative method using the first two
terms of the aforementioned Taylor Series:
$$f(x_{k+1})= f(x_k) + f'(x_k)(x_k-x_{k+1})$$
Since $f(x_{k+1}) = 0 $, the formula is then transformed into the Newton
Method:
@^analytic@> @^Quadratic@> @^Monotone@>
$$ x_{k+1} = x_k - {{f(x_k)}\over {f'(x_k)}} $$

@ This program is written to run under the Cray {\tt cft77} and
the usual compilers for SUNs.

@ This is controlled by the environment variables |CRAY| and |unix|.
It should be straightforward to accomodate other environments.

The {\tt ftangle} command line to produce |CRAY| code is:\hfil\break
\centerline{\tt ftangle Newton -m"CRAY" -d}
The {\tt -d} converts  |do|/|end do| constructions into labeled |do|
loops. Cray, like everybody else, has extended their compiler beyond
the standard, but their |do while| requires a labeled ending
statement.  Our coding
style assumes that the standards extended by MIL-STD 1752 are
present.  This added the |do|/|end do| and |implicit none|
statements to \Fortran.
@^32-bit 64-bit@>

@ Charles Karney of Princeton University furnished macros
for machine dependencies as part of a large body of code released
with John Krommes' \FWEB{} system.  The default is a system
that is typical of |SUN|'s |unix| systems. \FWEB{} actually takes care
of the |do/end do| but we use these macros to handle |implicit none|
and differences between systems whose native modes are 32-bit or 64-bit.
We would define another environment variable, |VMS|, if we were also
considering VAX-VMS targets.
@^32-bit 64-bit@>

@#if defined(CRAY)
@m CRAY 1
@m unix 0
@#else
@m CRAY 0
@m unix 1
@#endif


@ We start each \Fortran{} module with a command that defeats
the archaic default typing of variable names.  We don't suggest
that it should be removed from the compiler because of the
need for upward compatability from millions of lines of
existing code.  We question it remaining the default mode.

Most |unix| systems use |implicit undefined(a-z)|
while the other target environments
have a more reasonable form, |implicit none|.  We define |implicit_none|
in terms of a macro.  Our program uses this in its
\WEB{} form and the macro substitution affects the resulting
\Fortran{} code.

@f implicit_none implicit

@#if unix
@m implicit_none implicit undefined(a-z)
@#else
@m implicit_none implicit none
@#endif

@ We introduce another macro to make it easy to run the tests
with 64~bit precision in all the environments.  The Cray uses
64~bit precision for |real| variables by default.  The other
machines will use this when |real*8| is specified.  However,
the actual storage will vary among most of the machines.

Another macro is used to convert floating point constants
by appending the ``{\tt d0}'' exponent on 32-bit machines.
@^32-bit 64-bit@>

@f floating real

@#if CRAY
  @m floating real
  @m const(x) x
@#else
  @m floating real*8
  @m const(x) x##d0
  @#endif

@* A Newton method program.
The following program will use the Newton method to solve the equation
$\cos(x)=0$.  The {\it semicolons} at the end of each statement were
added automatically by \FWEAVE{} and denote the end of a statement.
These were added to make the code more sentence-like
({\it i.e. Literate Programming}) and only appear in the typeset
output.  We will add the |implicit_none| as previously described to avoid
errors in the |implicit| type casting found in \Fortran.

@a
    program cosine
      implicit_none
      @<Declare Variables@>
      @<Initialization of Variables@>
      @<Iteration Loop@>
      @<Print Results@>
      stop
      end

@ We declare the following variables: |x_0| will represent $x_k$
from the Taylor Series, |x| will correspond to $x_{k+1}$, and
|delta_x| will be calculated as the difference between $x_{k+1}$ and $x_k$.
These will be defined as |floating| which will be translated to the
|real| type in the \Fortran{} source code for the machine specified
in the {\tt FTANGLE} command line.

@<Declare Variables@>=
    floating x_0, x, delta_x

@ We will also declare two |integer| variables.  These are
an iteration index, |k|, and a |limit| for the number of iterations.

@<Declare Variables@>=
    integer k, limit

@ The initial value of |x| is not quite arbitrary.  This problem has
an $\infty$ of solutions and the initial value can cause convergence
to roots that are not desired.  In the case of |x| $= 1.2$, |x|
will converge to $\pi\over 2$.  We will also need to set |delta_x| to
any nonzero value.  This is needed to allow the iteration loop to
execute at least once.

Of course, in most reasonable programs, this would not be hardcoded,
there would be a dialog to establish these values.  The |limit| of
10 iterations and |delta_x_max|$=0.5$ are based upon knowing the problem
and algorithm's performance.

@<Initialization of Variables@>=
    x = const(1.2)
    delta_x=const(0.000001)
    limit = 10 /* more than safe */

@ The following is the heart of the Newton method.  We will continue
calculating $x_{k+1}$ until $\Delta x$ becomes sufficiently small.
With each iteration of the loop, the number of signicant digits
in |x| will double, approximately.

@<Iteration Loop@>=
    k = 1
    do while ((delta_x .gt. const(0.0)) .and. (k .le. limit))
      @<Calculate Newton change@>
      @<Apply convergence strategies?@>
      k = k + 1
      end do

@ In the vanilla Newton method, |x_0| will receive the value
of |x| from the previous iteration.  To calculate |x|, we use the Newton
formula described in Module 1 by substituting |x_0| for $x_k$.
Now, we will calulate |delta_x| to determine if the iteration loop
should be terminated.

@<Calculate Newton change@>=
       x_0 = x
       delta_x = - (cos(x_0))/(-sin(x_0))
       x = delta_x + x_0

@* Convergence Strategies.  It is well know that if certain conditions
are met, then the convergence of the Newton method is quadratic
and monotone.  This applies within the {\it radius of convergence}.
In some cases, we can apply strategies that will expand the radius of
convergence and still allow rapid convergence in the immediate proximity
of the solution.  The strategies will include constraining the step or
change and noticing when convergence is not monotone.
@^Monotone@>
@<Declare Variables@>=
      floating delta_x_max, delta_x_prev

@  These variables need initial values.

@<Initialization of Variables@>=
    delta_x_prev = const(0.0)
    delta_x_max = const(0.5)

@ When we solve practical problems, we will have some estimate of the
answer, or at least its order of magnitude.  After calculating
|delta_x|, we do not let its magnitude exceed |delta_x_max|.
We don't check for quadratic convergence.@^Quadratic@>
@<Apply convergence strategies?@>=
       if (delta_x > delta_x_max) then
         delta_x = delta_x_max
       else if (delta_x < (-delta_x_max)) then
         delta_x = -delta_x_max
         end if
@^Damping@>

@ This problem is an example of those problems that do not exhibit
monotone convergence.  In this case, the assumptions are not met
because there is an inflection point at the solution.  When we do
not have monotone convergence, we should note it.  We find this
by comparing the signs of consectutive changes.
@^Monotone@>
@<Apply convergence strategies?@>=
      if ((delta_x*delta_x_prev) .lt. const(0.0)) then
        write(*,*) 'Oscillating', k
@.Oscillating@>
        end if
      x = x_0 + delta_x
      delta_x_prev = delta_x

@ Finally, it's time to let the world know the results.
@.The solution ...@>
@<Print Results@>=
       write (*,*) 'The solution to cos(x)=0 is ', x

@* Index.