Re: floating point... it still puzzles me :-)
- From: Gordon Sande <g.sande@xxxxxxxxxxxxxxxx>
- Date: Thu, 18 May 2006 16:52:18 GMT
On 2006-05-18 12:55:42 -0300, Bart Vandewoestyne <MyFirstName.MyLastName@xxxxxxxxxx> said:
Hello group,
I think I need some advice... I am experiencing the following in my
numerical computations done in Fortran:
I plot the absolute error of an integration routine for certain
values of N. Since the value of the integral I compute is 1, the
absolute error shown is also the relative error. Now, for some
plots like the following, I see strange things happening at the
end, for large values of N:
http://www.cs.kuleuven.be/~bartv/stuff/convergence_graphs/double_precision_F/order_min3/niederreiter_f_ronald01_1000000p_1d_Laurie.jpg
http://www.cs.kuleuven.be/~bartv/stuff/convergence_graphs/double_precision_F/order_min3/niederreiter_f_ronald01_1000000p_5d_Laurie.jpg
http://www.cs.kuleuven.be/~bartv/stuff/convergence_graphs/double_precision_F/order_min3/niederreiter_f_ronald02_1000000p_1d_Laurie.jpg
http://www.cs.kuleuven.be/~bartv/stuff/convergence_graphs/double_precision_F/order_min3/niederreiter_f_warnock01_1000000p_1d_Laurie.jpg
http://www.cs.kuleuven.be/~bartv/stuff/convergence_graphs/double_precision_F/order_min3/niederreiter_f_warnock02_1000000p_1d_Laurie.jpg
http://www.cs.kuleuven.be/~bartv/stuff/convergence_graphs/double_precision_F/order_min3/sugihara_murota_f_warnock03_1000000p_1d_Laurie.jpg
Notice
how for large values of N, the typical convergence
behavior is broken. This always happens at about absolute errors
(which are in my case also relative errors) of 1e-13.
My compuations are done in double precision. If i switch to
quadruple precision, then this behavior is gone and the
convergence continuous smoothly.
Now my PhD promotor wants me to find out why I get this
behavior... I have a very very strong feeling that what I am
seeing here is related to the limitedness of floating point
arithmetic and double precision, but I don't know exactly how to
explain this phenomenon. I vaguely remember there is something
special about the value 1e-15 or 1e-16 in numerical computing,
and my strange behavior occurs for relative errors of 1e-13, so
I think I am close to an explanation, but i can't seem to get the
puzzle in place...
Anybody?
Best wishes,
Bart
Ultimately you are just evaluating a sum. The pieces come from a fancy
theory but it is still a sum.
So look up the numerical analysis of summation.
There are modern books with titles like "Accuracy of Numerical Compuatations"
and there are clasical books like "The Algebraic Eigenvalye Problem" and
others of its time.
With a bit of cancellation, as all your terms are not exactly 1, and a bit
of error in detrmining the terms you will end up with some error in the
sum. The books will tell you how to bound them. Roundoff (eps) is at
about 10^-16 and the terms may vary by X with N terms so it is plausible
to look at N X eps. The details are left as an exercise of the reader. ;-)
10^-16 x 30 x 30 ~ 10^-13 which is you level under my guesses of the number
of terms and their fluctuations in size.
OH, yes. Numerical analysis has its own newsgroup where you will get
better advice than in a programming newsgroup. It just happens that
Fortran users tend to know a lot about numerical analysis also. And
a few other hard disciplines as well.
.
- References:
- floating point... it still puzzles me :-)
- From: Bart Vandewoestyne
- floating point... it still puzzles me :-)
- Prev by Date: Re: floating point... it still puzzles me :-)
- Next by Date: Re: Proposal: NOEXTERNAL
- Previous by thread: Re: floating point... it still puzzles me :-)
- Index(es):
Relevant Pages
|
