[Python-ideas] PEP 485: A Function for testing approximate equality

[Chris Barker]

On Mon, Jan 26, 2015 at 9:33 AM, Chris Barker <chris.barker at noaa.gov> wrote:

> yes, but in an iterative solution you generally compute a solution, then
> use that to compute a new solution, and you want to know if the new one is
> significantly different than the previous -- so an asymmetric test does
> make some sense. But again either work work, and pretty much the same.
>
> Example forthcoming....
>

OK -- here is a real-life application -- probably a bit too much (and too
specialized) for the PEP, but for the record:

Waves in the ocean have a relationship between the wave length, wave
period, water depth (and the acceleration of gravity), known as the
dispersion relationship. So-called because the wave speed is a function of
the frequency and wave length, and so this relation also predicts the wave
speed -- and waves of different frequencies move at different speeds -- os
"disperse" as the move from the location they were generated -- but I
digress....

The relationship is:

  omega**2 = g*k*tanh(k*h)

where omega is the wave frequence, g is the accelartion of gravity, k
is the wave number (2*pi/wave_length), and h is the water depth.

In the usual case, the frequencey(omega) is known, and you want the
wave number. As the wave number appears in two places, there is no
direct solution, so a numerical method must be used.

The simplist iterative solution is something like this: recast the
equation as:

k_2 = omega**2 / (g * tanh(k_1 * h))

 - guess a k_1
 - compute k_2
 - check if k_2 is close to k_1
 - if not, set k_1 to k_2 repeat.

Here is the code for that:

def dispersion(omega, h, g=9.806):
    "compute the dispersion relation"

    k_1 = 10.0 # initial guess
    while True:
        k_2 = omega**2 / (g * math.tanh(k_1 * h))
        if is_close_to(k_2, k_1, tol=1e-5):
            break
        k_1 = k_2
    return k_1

note that I've provided g to only four significant figures, so I set the
tolerance to 1e-5 -- no need for more precision than that.

Granted, there are better ways to do this that run faster, and I'm sure
there is code out there to do it already (well, I know there is, I wrote
some of it...). In fact, I had code very much like this that I used in grad
school (not Python, though). Since then, I needed this a lot, and took the
time to write up a faster-converging method using Newton's method in C. But
frankly this works just fine, and any of the proposals on the table for
is_closr_to would work fine with it as well.

If I were a student and needed this, and is_close_to ws in the stdlib --
I'd probably use it.

In fact, you write a generic iteration routine:

def iterate(func, x_initial, *args):
    """
    iterate to find a solution to the function passed in

    func should be a function that takes x as a first argument,
    and computes an approximation to x as a result.

    x_initial is an initial guess for the unknown value
    """
    x_1 = x_initial
    while True:
        x_2 = func(x_1, *args)
        if is_close_to(x_2, x_1):
            break
        x_1 = x_2
    return x_2

Where you can pass in the function you want to iterate over. Here it is for
the above:

def disp(k, omega, h, g=9.806):
    """
    the linear  wave dispersion relationship

    k as a function of k, omega, h, g
    """
    return omega**2 / (g * math.tanh(k * h))

and I can then call it like this:

k2 = iterate(disp, 10, omega, h )

Code on gitHub if you care.

https://github.com/PythonCHB/close_pep/blob/master/iteration_example.py

-Chris



-- 

Christopher Barker, Ph.D.
Oceanographer

Emergency Response Division
NOAA/NOS/OR&R            (206) 526-6959   voice
7600 Sand Point Way NE   (206) 526-6329   fax
Seattle, WA  98115       (206) 526-6317   main reception

Chris.Barker at noaa.gov
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/python-ideas/attachments/20150126/41d048a8/attachment-0001.html>