I want to compute numerically the Taylor series coefficients of some (real analytic) function that appears in a problem I am looking at, let's say $ f(x) = \sum_{n\ge0}f_n x^n $.
I can do two things:
I have a way to compute $f(x_i)$ for any positive real $x_i<1$. This is a relatively quick computation and it can be done with arbitrary precision. I can do that for a "large" number of points (say something like 1000).
I have a way to compute directly any "f_n" with arbitrary precision, but this becomes very slow very quickly as $n$ increases. So realistically I can only compute a handful of $f_n$'s (say 10).
An obvious method to do that, is to calculate the first 10 coefficients and then get points of the form $p_i=(x_i,f(x_i))$ and fit the truncated Taylor series for the "rest".
Is there any better method than this obvious one? Is there any rule of thumb about the precision I need for $p_i$'s if I want a given precision of $f_n$'s? How many $f_n$'s is it reasonable to expect to calculate accurately (for some given precision)?
Rather than an answer this is just to make my comment more understandable, and I need some lines for that.
I entangled in a similar problem and had to make clear to myself that:
- the (truncated) Taylor series is a polynomial "osculating" (that has the highest contact order) the function at $x=0$, or any other fixed $x_0$;
- the interpolating polynomial exactly matches the function at various different values of $x$ at discrete distance from each other.
That is, in the Taylor polynomial the interpolation points are infinitely close to each other, and, for a finite degree, it will never interpolate exactly other points at finite distance from $x_0$.
In the reverse, it also means that an interpolating polynomial will never provide you the exact Taylor coefficients, except for the initial term $f(x_0)$.
The bounds on error will be given by the various Remainder formulations, and thus depend on the actual $f(x)$ being considered.
A classical example is $e^x$.