On Mon, Jan 26, 2015 at 9:33 AM, Chris Barker <chris.barker@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@noaa.gov