Numeric root-finding in Python

[Dave Angel]

On 02/12/2012 06:05 PM, Steven D'Aprano wrote:
> On Sun, 12 Feb 2012 12:18:15 -0800, Mark Dickinson wrote:
>
>> On Feb 12, 6:41 am, Steven D'Aprano<steve
>> +comp.lang.pyt... at pearwood.info>  wrote:
>>
>>>              err = -a/b  # Estimate of the error in the current w.
>>>              if abs(err)<= 1e-16:
>>>                  break
>> If the result you're expecting is around -1.005, this exit condition is
>> rather optimistic:  the difference between the two Python floats either
>> side of this value is already 2.22e-16, so you're asking for less than
>> half a ulp of error!
> I was gradually coming to the conclusion on my own that I was being
> overly optimistic with my error condition, although I couldn't put it
> into words *why*. Thanks for this Mark, this is exactly the sort of thing
> I need to learn -- as is obvious, I'm no expert on numeric programming.
>
>
me either.  But comments below.
>> As to the rest;  your error estimate simply doesn't have enough
>> precision.  The main problem is in the computation of a, where you're
>> subtracting two almost identical values.<SNIP>
> <SNIP>
>
Two pieces of my history that come to mind.  40+ years ago I got a 
letter from a user of our computer stating that our math seemed to be 
imprecise in certain places.   He was very polte about it, and admitted 
to maybe needing a different algorithm.  The letter was so polite that I 
(as author of the math microcode) worked on his problem, and found the 
difficulty, as well as a solution.

The problem was figuring out the difference in a machining table between 
being level every place on its surface (in which case it would be 
slightly curved to match the earth), or being perfectly flat (in which 
case some parts of the table would be further from the earth's center 
than others)  The table was 200 feet long, and we were talking 
millionths of an inch.  He solved it three ways, and got three different 
answers.  The first two differed in the 3rd place, which he thought far 
too big an error, and the third answer was just about exactly half the 
others.

Well the 2:1 discrepancy just happens when you change your assumption of 
what part of the flat table is level.  If the center is level, then the 
edges are only 100 feet out, while if the edge is level, the other edge 
is 200 feet out.

But the other solution was very interesting.  Turns out he sketched a 
right triangle, with narrow angle at the center of the earth, side 
opposite being 200 feet.  He then calculated the difference between the 
other two sides.  one 8000 miles, and the other 8000 miles plus a few 
microinches.  He got that distance by subtracting the sine from the 
tangent, or something similar to that.  I had microcoded both those 
functions, and was proud of their accuracy.  But if you subtract two 13 
digit numbers that only differ in the last 3, you only get 3 digits 
worth of accuracy, best case.

Solution was to apply some similar triangles, and some trivial 
approximations, and the problem turned out not to need trig at all, and 
accurate to at least 12 places.  if I recall, it was something like  
8000mi is to 200 feet, as 200 feet is to X.  Cross multiply and it's 
just arithmetic.


The other problem was even earlier.  It was high school physics, and the 
challenge was to experimentally determine the index of refraction of air 
to 5 places.  Problem is our measurements can't be that accurate.  So 
this is the same thing in reverse.  Find a way to measure the difference 
of the index of refraction of air and vacuum, to one or two places, and 
add that to 1.00000000

taken together with lots of other experience, i try to avoid commiting 
an algorithm to code before thinking about errors, convergence, and 
exceptional conditions.  I've no experience with Lambert, but I suspect 
it can be attacked similarly.





-- 

DaveA