Accuracy and error analysis (was Re: [Maxima] primes)

--- Richard Fateman  wrote:
> ok, long explanation.
> Consider the maxima definition
> 
> f[i,n]:=f[i-1,n]-(f[i-1,n]^2-n)/(2*f[i-1,n]);
> 
> This is a newton iteration for square root of n
> Compute the square root of 1+d this way..
> f[0,3]:1;
> f[9,3];  --> 9 iterations towards sqrt(3).

OK.

> next.. sqrt(1+d)
> 
> f[0,1+d]:1   /* initial guess */
> 
> f[2,1+d];   /*after 2 iterations */
> 
> (d/2) - (( - d + ((d/2) + 1)^2 - 1)/(2 * ((d/2) + 1))) + 1
> 
> which you can understand by
> 
> taylor(%,d,0,3)  --> 1 +d/2 -d^2/8 + ....
> 
> so it is getting what I said it should get.

OK.

> But look again at the formula. It takes advantage of
> the fact that each occurrence of "d" is the same "d", and
> thus  d-d/2 is  d/2.  and other formulas might
> make use of simplifications like (1-d)*(1+d) = 1-d^2.

But d as an error shouldn't be used with the Newton iteration - that
does not accurately represent the uncertainty involved.  I'm afraid now
I'm even more confused, because I don't see how you can Newton iterate
towards an inexact number.  There is a d associated with each operation
in the Newton equation resulting in an inexact answer:

f[i-1,n]^2[d1]
(f[i-1,n]^2-n)[d2]
(2*f[i-1,n])[d3]
((f[i-1,n]^2-n)/(2*f[i-1,n]))[d4]
(f[i-1,n]-(f[i-1,n]^2-n)/(2*f[i-1,n]))[d5]

but this doesn't impact the final answer as such since the next
iteration will use the exact return as the next seed.  It could equally
well use either return+dtotal or return-dtotal, since either one will
also iterate toward the same final answer.  It does, however matter
when the propagated error of one Newton calculation (d5) is larger than
the difference between your previous value and your current one.  At
that point, Maxima has reached the limit of its ability to accurately
calculate the sqrt with this method and current ffprec settings. 
Presumably with more precision in floating point, each individual d
would be smaller, resulting in a smaller d5 and the possibility of more
precise Newton calculations.

The "relevant d" to include with the final answer here would be (at a
guess) (either .5 or 1)*10^(power of final unchanging digit).  This or
something like it is the actual uncertainty of the final answer.  Then,
if an a+d type number was supplied, the d would be delt with once the
Newton iteration had been completed (not sure how to combine offhand -
possibly Newton(a)+-(abs(Newton(a-d)-(relevant d of a-d
calc)-(Newton(a+d)+(relevant d from a+d calc)))) would be workable, if
a bit expensive for nontrivial cases.
 
> In general, if maxima is doing
> interval arithmetic with uncertainties, there is not just
> one d.  Each d is different.  e.g. if you are talking about
> one unit in the last place, then d might be a shorthand for something
> like any number in [-1.0e-18, +1.0e-18]

Right.  I associate 1+d with an original measurement uncertainty
entered by the user, not the propagated error resulting from the
numerical limits of Maxima calculations.

> In such a case  d[1] -d[2]/2  is  not d[3]/2  but a number in
> an interval like  [d/2, 3d/2].
> 
> and (1-d)^2  is   (1-d[1])*(1-d[2]).  Instead of getting
> 1-2d+d^2,  you get something like +-(1+2|d| +d^2), which is
> cruder than you'd like.

Yes, that would be an overestimation.

> I hope I've conveyed a little of the problem. 

Of using Newton to propagate the error, yes.  I'm not convinced that's
relevant to the error of interest, but I readily concede I'm not
qualified to have an opinion.

> Oh, and solving the newton iteration problem
> so that you just solve that
> one example does not automate all of error analysis...

Of course not.  But if nothing else trying to understand the problems
involved seems to be a useful way to learn about error analysis in
general.  It remains an issue that uncertainties introduced by CAS
calculations will impact invisibly answers to certain classes of
problems.  If a practical solution (even a difficult to impliment one)
exists it would be useful to know about, in case it ever becomes worth
taking a serious run at.  Or, if it is definitely impossible (in which
case I question a human beings ability to do it as well) then that's
also good to know.

Thanks for explaining these things, Dr. Fateman.  I really appreciate
your taking the time.
CY


		
Yahoo! Mail
Stay connected, organized, and protected. Take the tour:
http://tour.mail.yahoo.com/mailtour.html