Floating point arithmetic: $(x-2)^9$

Question

This is taken from Trefethen and Bau, 13.3.

Why is there a difference in accuracy between evaluating near 2 the expression $(x-2)^9$ and this expression:

$$x^9 - 18x^8 + 144x^7 -672x^6 + 2016x^5 - 4032x^4 + 5376x^3 - 4608x^2 + 2304x - 512 $$

Where exactly is the problem?

Thanks.

Answer 1

$(x-2)$ is small by definition of $x$, $(x-2)^9$ is even much smaller but can be computed with small relative error. The single terms of the expanded polynomial are much larger and therefore you will suffer from catastrophic cancellation (see e.g. Wiki or do a web search). Example: $$x=2 + 10^{-2} \Longrightarrow (x-2)^9 = 10^{-18}.$$ This is below the machine epsilon for double! Now with your expanded polynom you have to compute $-512 + 23.04 \pm \dots$.

And here is the actual computation with double and evaluation of the polynomial using Horner for $x=2.01:$

(x-2)^9 = 9.99999999999808E-019      poly(x) =  -3.75166564481333E-012

As expected the polynomial result is completely wrong (in terms of relative error), but note that even the relative error for the first is about $2\cdot 10^{-13},$ which is a result from the fact that $0.01$ cannot represented exactly as double.

Answer 2

I have tried to plot this 2 expresions near x=2.0. Look at the image. The behaviour ( shape) seems very different .

Here is the Maxima CAS code :

draw2d(
  file_name = "d",
  color = blue,
  key = "second function",
  explicit(x^9 - 18*x^8 + 144*x^7 -672*x^6 + 2016*x^5 - 4032*x^4 + 5376*x^3 - 4608*x^2 + 2304*x - 512 ,x, 2-0.0000005,2.0000005),

  color = red,
  key = " y = (x-2)^9",
  explicit((x-2)^9, x, 2-0.0000005,2.0000005),

  terminal = 'svg) $   

Note that above 2 expressions are the same !

expand((x-2)^9)

So problem is in computer arithmethic not in math ( informally speaking) More : stable evaluation of polynomials . See also discussion