<div dir="ltr"><br><div class="gmail_extra"><br><div class="gmail_quote">On Mon, Jan 26, 2015 at 9:33 AM, Chris Barker <span dir="ltr"><<a href="mailto:chris.barker@noaa.gov" target="_blank">chris.barker@noaa.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><div>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.</div><div><br></div><div>Example forthcoming....</div></div></div></div></blockquote><div><br></div><div>OK -- here is a real-life application -- probably a bit too much (and too specialized) for the PEP, but for the record:</div><div><br></div><div>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....</div><div><br></div><div>The relationship is:</div><div><br></div><div><div> omega**2 = g*k*tanh(k*h)</div><div><br></div><div>where omega is the wave frequence, g is the accelartion of gravity, k</div><div>is the wave number (2*pi/wave_length), and h is the water depth.</div><div><br></div><div>In the usual case, the frequencey(omega) is known, and you want the</div><div>wave number. As the wave number appears in two places, there is no</div><div>direct solution, so a numerical method must be used.</div><div><br></div><div>The simplist iterative solution is something like this: recast the</div><div>equation as:</div><div><br></div><div>k_2 = omega**2 / (g * tanh(k_1 * h))</div><div><br></div><div> - guess a k_1</div><div> - compute k_2</div><div> - check if k_2 is close to k_1</div><div> - if not, set k_1 to k_2 repeat.</div></div><div><br></div><div>Here is the code for that:</div><div><br></div><div><div>def dispersion(omega, h, g=9.806):</div><div> "compute the dispersion relation"</div><div><br></div><div> k_1 = 10.0 # initial guess</div><div> while True:</div><div> k_2 = omega**2 / (g * math.tanh(k_1 * h))</div><div> if is_close_to(k_2, k_1, tol=1e-5):</div><div> break</div><div> k_1 = k_2</div><div> return k_1</div></div><div><br></div><div>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.</div><div><br></div><div>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.</div><div><br></div><div>If I were a student and needed this, and is_close_to ws in the stdlib -- I'd probably use it.</div><div><br></div><div>In fact, you write a generic iteration routine:</div><div><br></div><div><div>def iterate(func, x_initial, *args):</div><div> """</div><div> iterate to find a solution to the function passed in</div><div> </div><div> func should be a function that takes x as a first argument,</div><div> and computes an approximation to x as a result.</div><div> </div><div> x_initial is an initial guess for the unknown value</div><div> """</div><div> x_1 = x_initial</div><div> while True:</div><div> x_2 = func(x_1, *args)</div><div> if is_close_to(x_2, x_1):</div><div> break</div><div> x_1 = x_2</div><div> return x_2</div></div><div><br></div><div>Where you can pass in the function you want to iterate over. Here it is for the above:</div><div><br></div><div><div>def disp(k, omega, h, g=9.806):</div><div> """</div><div> the linear wave dispersion relationship</div><div><br></div><div> k as a function of k, omega, h, g</div><div> """</div><div> return omega**2 / (g * math.tanh(k * h))</div></div><div><br></div><div>and I can then call it like this:</div><div><br></div><div>k2 = iterate(disp, 10, omega, h )<br></div><div><br></div><div>Code on gitHub if you care.</div><div><br></div><div><a href="https://github.com/PythonCHB/close_pep/blob/master/iteration_example.py">https://github.com/PythonCHB/close_pep/blob/master/iteration_example.py</a><br></div><div><br></div><div>-Chris</div><div><br></div><div><br></div><div> </div></div>-- <br><div class="gmail_signature"><br>Christopher Barker, Ph.D.<br>Oceanographer<br><br>Emergency Response Division<br>NOAA/NOS/OR&R (206) 526-6959 voice<br>7600 Sand Point Way NE (206) 526-6329 fax<br>Seattle, WA 98115 (206) 526-6317 main reception<br><br><a href="mailto:Chris.Barker@noaa.gov" target="_blank">Chris.Barker@noaa.gov</a></div>
</div></div>