Numerical evaluation of monomial divided differences

Question

Suppose $f(x)=x^{n+1}$ for some $n\in\mathbb{N}$, and define the divided difference $$f[a,b]=\frac{a^{n+1}-b^{n+1}}{a-b}.$$

I am wondering about the best way to numerically evaluate $f[a,b]$ to high precision, when

i) $n$ may be large

ii) $a,b\in(-1,1)$ or $a,b\in (0,1)$

iii) $a$ may be close to or even equal to $b$, with $f[a,b]$ replaced with $f'(b)=(n+1)b^n$ at $a=b$.

I know

$$f[a,b]=\sum_{k=0}^{n}a^{k}b^{n-k}$$

and this seems like it should be numerically stable but slow for large $n$. For $a$ moderately but not extremely close to $b$, it seems like replacing $f[a,b]$ with $f'(b)$ might lead to a large truncation error. I have also thought about using a Taylor expansion for the numerator of $f[a,b]$ and writing $f(a)=\sum_{k=0}^{m}f^{(k)}(b)(a-b)^{k}/k!+f^{(m+1)}(\xi)(a-b)^{m+1}/(m+1)!$ which would allow cancelling the $a-b$ in the denominator to get

$$f[a,b]=\sum_{k=1}^{m}f^{(k)}(b)(a-b)^{k-1}+\dbinom{n+1}{m+1}(\xi)^{n-m}(a-b)^{m}$$

but I'm not sure how to get a usable bound for the error term. If this problem has been studied before, a reference would also be appreciated.

Answer 1

Define $M=\pmatrix{a&1\\0&b}$.

In general $f[a,b]$ is the top-right entry of $f(M)$ where $f$ is the function $f$ extended to matrices. (See Higham's book for info on functions of matrices.)

In your case you want the top right element of $M^{n+1}$ which you can get by multiplication and repeated squaring in $O(\log n)$ steps. There are no subtractions involved so it should be numerically reasonable.

You can make it more efficient - no need to compute and use that zero in the bottom left. This is analogous to automatic differentiation - which it becomes when $a=b$.

Also see Kahan and Fateman.

Answer 2

Just key Idea not an answer to your question:

Given $f(x)=x^{n+1}$, the divided difference $f[a,b]=\frac{a^{n+1}-b^{n+1}}{a-b}$ can be written as $f[a,b]=\sum_{k=0}^{n}a^{k}b^{n-k}$. For $a=b$, use $f[a,b]=f'(b)=(n+1)b^n$. When $a$ is close to $b$, use the Taylor expansion: $f[a,b]=\sum_{k=1}^{n+1} \binom{n+1}{k} b^{n+1-k} (a-b)^{k-1}$. To ensure numerical stability and efficiency for large $n$, compute $f[a,b]$ using Horner's method: $f[a,b] = b^n + a(b^{n-1} + a(b^{n-2} + \cdots + a(b + a)\cdots))$. For reference i may suggest to you Numerical Analysis, 9th by Richard L. Burden and J. Douglas