Neat Fortran/Numerical tricks
Author Message Neat Fortran/Numerical tricks

Hello fellow netters:

I'm teaching a numerical methods (more exactly a numerical tricks)
seminar in a couple months and I'm looking for some added input.
Specifically, what are some of the neat little numerical tricks you use
A couple of examples are:
1)      use
pi=4.0*atan(1.0)
because you don't have to remember pi to 15 or 30 decimal places
and it will automatically be defined to the right precision
necessary for the system & compiler your using.

2)      use Horner's method to turn a polynomial in to products and sums
a*x**4 + b*x**3 + c*x**2 + d*x + e =
(((a*x + b)*x + c)*x + d)*x + e
This eliminates the exponentials and speeds the computation.

The seminar isn't too long so I don't think I can get into more
esoteric material, but what I would like are the coding tricks or
design considerations that other people on the net implement to
obtain accurate and fast results.

Thanks,
R.K.

+-----------------------------------+----------------------------------+

| Central Computer Facility         |  FAX:(415)964-1760               |
| NASA-Ames Research Center         |  WK: (415)604-2677               |
+-----------------------------------+----------------------------------+

Sun, 06 Nov 1994 07:39:27 GMT  Neat Fortran/Numerical tricks

Quote:
(R.K.Owen) writes:
>Hello fellow netters:
>  I'm teaching a numerical methods (more exactly a numerical tricks)
>seminar in a couple months and I'm looking for some added input.
>Specifically, what are some of the neat little numerical tricks you use
>A couple of examples are:
>1)     use
>               pi=4.0*atan(1.0)
>       because you don't have to remember pi to 15 or 30 decimal places
>       and it will automatically be defined to the right precision
>       necessary for the system & compiler your using.
>2)     use Horner's method to turn a polynomial in to products and sums
>               a*x**4 + b*x**3 + c*x**2 + d*x + e =
>               (((a*x + b)*x + c)*x + d)*x + e
>       This eliminates the exponentials and speeds the computation.
>The seminar isn't too long so I don't think I can get into more
>esoteric material, but what I would like are the coding tricks or
>design considerations that other people on the net implement to
>obtain accurate and fast results.
>       Thanks,
>       R.K.
>+-----------------------------------+----------------------------------+

>| Central Computer Facility         |  FAX:(415)964-1760               |
>| NASA-Ames Research Center         |  WK: (415)604-2677               |
>+-----------------------------------+----------------------------------+

Horner's Method (which I have not seen called that) I teach in my Math
Methods of Physics class (to Juniors) both as a Numerical and Analytical
Technique - namely this is the same as what is often called "Synthetic
Division" for evaluating/factoring polynomials.

that even simple operations such as incrementing a time t(i) = t(i-1) + dt
should be thought out carefully:
t(i) = t(i-1) + dt  is in general faster than  t(i) = t(0) + i*dt  (at least
on most machines where mult is slower than addition) but the round off
error accumulates.  Similarly I at least mention things like hard-coding
loops for speed but with trade off of larger code size.

The algorithm for finding roots to quadratic eqn is often very unstable -
see Numerical Recipes by Press et al for a better algorithm

If you do root finding, I have several "fun" examples which I do in class (I
bring a computer in and project the results with an LCD panel)
For Newton's Method
Sin(x)  If x0 <~ 1 will converge to x=0
If x0 = 1.4 will converge to x=pi
If x0 = 1.5 will converge to x=-4pi
ArcTan(x)  If x0 = 0.1 will converge to x=0
If x0 = 0.2 it will blow up
Sign(sqrt(x),x)  Will oscillate for any X value
I teach them to always use (and have them write as an assignment) a combined
bisection/newton or bisection/secant method so they get guaranteed
convergence and often the speed of the better algorithms.

-------------------

St. Peter, MN 56082 USA        Phone:     (507)933-7036
-------------------

St. Peter, MN 56082 USA        Phone:     (507)933-7036

Sun, 06 Nov 1994 23:15:49 GMT  Neat Fortran/Numerical tricks

Quote:

>[for finding real zeros of a scalar function]
>I teach them to always use (and have them write as an assignment) a combined
>bisection/newton or bisection/secant method so they get guaranteed
>convergence and often the speed of the better algorithms.

I would suggest teaching them *about* the (bisection/secant or
bisection/Newton) method, but having them program using a existing
high quality "black box" code which implements the method carefully.
The books by Forsythe, Malcolm, and Moler, and more recently by
Kahaner, Moler, and Nash, exemplify this philosophy.

The key point here is that your students will probably continue to
use the programs from this course for years to come, in a wide range
of problems.  It's thus important for them to become familiar with
*high* *quality* codes.  These need not be any harder to use than
ones of lesser quality, and fortunately are often available (free!)

For the (Van Wijngaarden-Dekker-Brent bisection/secant) method under
discussion, the best implementation seems to be the standard  ZEROIN
code available from netlib.  (Forsythe, Malcolm, and Moler give both
a detailed exposition and a complete source listing of this code.)
Note that the version of this code in Numerical Recipes has been
modified from the original, and can suffer from spurious floating
point underflow.

University of Texas at Austin / Physics Dept / Center for Relativity
and (for a few more months) University of British Columbia / Astronomy

Mon, 07 Nov 1994 01:57:37 GMT

 Page 1 of 1 [ 4 post ]

Relevant Pages