Application Center - Maplesoft

App Preview:

An Example Of Round-Off Error

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application




Bruno Guerrieri

Florida A&M University

The following is in no way an original idea.  I heard about it, many years ago, from my major professor, Dr. Christopher Hunter (Florida State University), and, at a later date, have come across essentially the same idea in "Computer Methods for Mathematical Computations" by G.E. Forsythe, M.A. Malcolm, C.B. Moller, Prentice-Hall (1977), page 16.  The results you will be getting on your machine will, perhaps, be different from ours, due to different machine accuracy, but we hope that our presentation will remain consistent.

>    restart:

The quantity "alpha" coming up is an INTEGER.  We have only used it with the values 25, 30, and 40.  You are welcome to experiment but at your risks and perils!

>    alpha := 40:

Consider the following integral

>    Quest := Int(x^alpha/(x+8),x=0..1);

Quest := Int(x^40/(x+8),x = 0 .. 1)

 Before we evaluate this integral, let us consider a string of inequalities that could be unnerving but the exercise has some value:

abs(Int(x^alpha/(x+8),x = 0 .. 1)) <= Int(abs(x^alpha/(x+8)),x = 0 .. 1)  which can be shown to be less than 1/7 using a variation of the triangle inequality such abs(abs(x)-abs(8)) <= abs(x-(-8))   and using the fact that x is between 0 and 1.  This will help us to realize that some upcoming answers cannot be correct.

Maple probably thinks nothing of evaluating the integral!

>    Quest = int(x^alpha/(x+8),x=0..1);

Int(x^40/(x+8),x = 0 .. 1) = 2658455991569831745807614120560689152*ln(3)-3987683987354747618711421180841033728*ln(2)-209123016346760841446125740863301305248944644727499/1335732864265800

>    Quest = evalf(int(x^alpha/(x+8),x=0..1));

Int(x^40/(x+8),x = 0 .. 1) = .16e28

Sometime, depending on the value of alpha, the above answer comes up negative.  A little bit worrisome in that case since the integral should come out as a positive number!  

Plotting the (positive for x between 0 and 1) integrand would confirm that.  It is because of these negative answers and also with the upper bound seen earlier that we are very leary of an answer such as the one above.

>    plot(x^alpha/(x+8),x=0..1);

[Maple Plot]

What about finding the antiderivative and then applying the Fundamental Theorem of Calculus?  Would we get the same answer?

>    antideriv := int(x^alpha/(x+8),x);

antideriv := -166153499473114484112975882535043072*x+1/40*x^40-8589934592/29*x^29-549755813888/27*x^27+17179869184/7*x^28-32768/35*x^35-2097152/33*x^33+131072/17*x^34+524288*x^32+536870912/15*x^30-1342...
antideriv := -166153499473114484112975882535043072*x+1/40*x^40-8589934592/29*x^29-549755813888/27*x^27+17179869184/7*x^28-32768/35*x^35-2097152/33*x^33+131072/17*x^34+524288*x^32+536870912/15*x^30-1342...
antideriv := -166153499473114484112975882535043072*x+1/40*x^40-8589934592/29*x^29-549755813888/27*x^27+17179869184/7*x^28-32768/35*x^35-2097152/33*x^33+131072/17*x^34+524288*x^32+536870912/15*x^30-1342...
antideriv := -166153499473114484112975882535043072*x+1/40*x^40-8589934592/29*x^29-549755813888/27*x^27+17179869184/7*x^28-32768/35*x^35-2097152/33*x^33+131072/17*x^34+524288*x^32+536870912/15*x^30-1342...
antideriv := -166153499473114484112975882535043072*x+1/40*x^40-8589934592/29*x^29-549755813888/27*x^27+17179869184/7*x^28-32768/35*x^35-2097152/33*x^33+131072/17*x^34+524288*x^32+536870912/15*x^30-1342...
antideriv := -166153499473114484112975882535043072*x+1/40*x^40-8589934592/29*x^29-549755813888/27*x^27+17179869184/7*x^28-32768/35*x^35-2097152/33*x^33+131072/17*x^34+524288*x^32+536870912/15*x^30-1342...

 We are not sure we expected the answer to be what it is, so let us verify it by calculating the derivative, in the hope of recovering the original integrand.

>    simplify(diff(antideriv,x));


>    evalf(subs(x=1,antideriv)) - evalf(subs(x=0,antideriv));


Should we check using numerical integration?

>    evalf(Quest);


 Since we have, so far, three distinct different answers.  We could forge ahead with a different method, a longer and less intuitive one, but one that we can keep an eye on as we proceed.   Let us define the function   A(n)  = Int(x^n/(x+8),x = 0 .. 1)   and  note that   A(n)+8*A(n-1)  =   Int(x^n/(x-8),x = 0 .. 1)+8*Int(x^(n-1)/(x-8),x = 0 .. 1)  = Int(x^(n-1)*(x+8)/(x+8),x = 0 .. 1)   = Int(x^(n-1),x = 0 .. 1)   = 1/n   

  Since A(0)  = Int(1/(x+8),x = 0 .. 1)   =  ln(9/8), we could "recursively" apply the formula A(n) = -8*A(n-1)+1/n   up to n = alpha.

Let us solve the recursion relation symbolically, using the power of Maple (cf. Heck's "Introduction to Maple",Springer-Verlag, (2nd Ed.))

>    Recurrence := rsolve({A(n) = -8*A(n-1) + 1/n,A(0) = ln(9/8)},A(n));

Recurrence := (2*ln(3)-3*ln(2))*(-8)^n+Sum((-8)^(n-n0)/n0,n0 = 1 .. n)

>    evalf(subs(n=alpha,Recurrence));


 We obtain an answer we already came across before but this is not to be taken as a confirmation.  Let us, for verification purposes, also evaluate the recursion formula numerically

>    Ans := ln(9/8);

>    for n from 1 to alpha do Ans := evalf(-8*Ans + 1/n); od:

>    Ans;

Ans := ln(9/8)


 Well, enough wild-goose chases! We are not quite sure what Maple is doing with respect to the integration process we carried out at the beginning but, as far as the recursion process is concerned, what could possibly have gone wrong?  Well, you guessed it.  Round-off errors swallowed us!  But how? Well, you recall the initial value A(0) = ln(9/8)?  This number was not quite the number that was stored in the variable "Ans", a few lines above, since we are not working with infinite precision here.  In fact, let us run the following procedure (borrowed from Marron's "Numerical Analysis", MacMillan (2nd Ed.)) which will tell us the level of precision of our (and your) machine.

>    Machine_precision := proc()

>        local k,u,test;
    u := 1; test := 1 + u;

>        for k from 1 to 100 while test-1 <> 0 do

>           u := u/2;

>           test := evalf(1 + u);

>        od;

>    end:

>    Machine_precision();


 The above calculation (if you dissect it) was adding to the number 1 the quantity 2^(-k)  and, by the time we were adding 2^(-32)  for our machine, we could not tell the difference (of course, the number may be different for your machine and everything we are saying from now on will have to take this into account).  Another way to put it is to say that the quantity that our computer stored for ln(9/8) was roughly off by 2^(-32)  .  After that, this rather minute error was AMPLIFIED repetitively, everytime we multiplied the "A(n) at bat" by a factor of 8.  A quick calculation, trying to measure the "evolution" of that error (we are multiplying by 8, alpha many time) tells us the following:

>    Error := evalf(2^(-32)*8^(alpha));

Error := .3094850098e27

 Here, let us make a feeble attempt at saying that the value for "Ans" and "Error" are within the same ballpark (in absolute value and as we could expect) but you probably will not be convinced.  And anyway, the more important question is whether we will ever obtain the right answer, or, at least, a satisfactory one!  Well, here is the grand finale.  We are actually going to make the roundoff error work for us(!!!).  Let us revisit the recursion formula seen earlier and solve it for A(n-1) .  We obtain   A(n-1) = -A(n)/8+1/(8*n)   .  All we have to do now is start at, say A(2*alpha) and work our way DOWN to A(alpha).  Think of alpha as being 30 when you are reading this.

 Excuse me, but what is A(2*alpha) equal to, may we ask?  Here is the beauty of the whole thing.  In a way, the value that A(2*alpha) assumes almost does not matter because we will be dividing the possible error of a wild guess by a factor of 8 at each round!  So, for simplicity, let us assume that  A(2*alpha) = 1000.0.  Here we go!  if you wish to start further up the ladder for as better clean-up, then change the value of "start" to that "higher up" choice.

>    start := 2*alpha;

start := 80

>    Descending_Iterate[start] := 1000.0:

>    for n from 1 to start-alpha+1 do

>       Descending_Iterate[start-n] := - Descending_Iterate[start-n+1]/8 + 1/(8*(start-n+1)):

>    od:

>    NewTable := [seq([start-n+1,Descending_Iterate[start-n+1]],n=1..start-alpha+1)]:

>    array(1..start-alpha+1,1..2,NewTable);

matrix([[80, 1000.0], [79, -124.9984375], [78, 15.62638697], [77, -1.951695807], [76, .2455853525], [75, -.2905343222e-1], [74, .5298345695e-2], [73, .1026895977e-2], [72, .1583966770e-2], [71, .153811...

 Compare that with an earlier answer:

>    evalf(Quest);





 So, now we know what the exact answer to 10s ( maybe, depending on the parameter alpha you use) is!