This worksheet demonstrates the error in evaluating a polynomial by directly computing the powers vs. evaluating a nested version. 

>	with(plots):

 

>

 

>	a:=-6.1; b:=3.2; c:=1.5;

 

>	f:=x->x^3+a*x^2+b*x+c;

 

>

f(4.71);

 

Above are the function and x-value from the book. 

Below, we use evalf[j] to simulate the computation of f to j digits of accuracy. Rounding takes place after every computation. Check that the parentheses match up. 

>	froundpower:=(j,x)->evalf[j](evalf[j](evalf[j](evalf[j](evalf[j](x*x)*x) + evalf[j](a*evalf[j](x*x)))+evalf[j](b*x))+c);

 

>	froundpower(3,4.71);

 

With three significant digits (notice: not three decimal places), we have a significant error using the direct method. 

Below, we compute the same polynomial using the nesting method. Observe that there are only five computations, as opposed to the eight required with the direct method. 

>	froundnest:=(j,x)->evalf[j](c+evalf[j](x*evalf[j](b+evalf[j](x*evalf[j](a+x)))));

 

>	froundnest(3,4.71);

 

This is much closer to the correct result of -14.258.  

The next command plots values of x between 0 and N, at intervals of 0.001, and the corresponding function values, as computed by the direct method. 

>	powereval:=(j,N)->pointplot([seq([0.001*k,froundpower(j,0.001*k)], k=0..N*1000)], symbol=POINT, color=blue);

 

The next command does the same thing with the nesting method. 

>	nesteval:=(j,N)->pointplot([seq([0.001*k,froundnest(j,0.001*k)], k=0..N*1000)], symbol=POINT,color=green);

 

The next command is the ordinary plot of the function to high accuracy.  

>	pureeval:=N->plot(x^3+a*x^2+b*x+c,x=0..N, thickness=2, color=red);

 

>	computationdisplay:=(j,N)->display(powereval(j,N),nesteval(j,N),pureeval(N));

 

Here we display the "pure" plot in red, the direct computation in blue, and the nesting computation in green. We model a three-digit decimal approximation on the interval [0,6]. Observe the approximation methods seem to have the most trouble near a turning point of the function.  

>	computationdisplay(3,6);

 

To better see what's going on, we plot the differences between the computation with roundoff and the true result. First the direct method. 

>	powererror:=(j,N)->pointplot([seq([0.001*k,froundpower(j,0.001*k)-f(0.001*k)], k=0..N*1000)], symbol=POINT, color=blue);

 

Then the nesting method. 

>	nesterror:=(j,N)->pointplot([seq([0.001*k,froundnest(j,0.001*k)-f(0.001*k)], k=0..N*1000)], symbol=POINT,color=green);

 

>	errordisplay:=(j,N)->display(powererror(j,N), nesterror(j,N));

 

We display the error in the direct method (in blue) and in the nesting method (in green) for a three-decimal roundoff on the interval [0,6]. Observe the nesting method for this polynomial is virtually always better by a significant margin. It also seems to do much better at the turning point. 

>	errordisplay(3,6);

 

This is not a totally fair comparison; if one changes the parameters slightly, one gets a very different picture. 

>	errordisplay(2,10);

 

We see that sometimes (for example around x=8) the direct method does much better than the nesting method, and one is not necessarily better than the other. Nevertheless, as a general rule it is often better to use nesting. This is partly for accuracy, but mostly to save time in computation. Here nesting takes five computations while direct evaluation takes eight; as the size of the polynomial grows, this difference becomes more drastic. If you're evaluating the function many times, that will add up in the long run.  

Also, while it is not always the case, typically one expects each new computation to introduce a slight error. In the second error comparison, the green points are mostly closer to zero than the blue points due to error buildup. Around x=8, the multiple errors in the direct method happen purely by coincidence to cancel out, while they happen to add up in the worst way for the nesting method. Despite the possibility of such things, it's almost always better to use nesting.