# Approximation:- Algorithmic considerations

Hello

I want to approximate a function $f$ on $(a,b)$. The function is singular at the points $a$ and $b$, however I have asymptotic expansions at these points. I can also construct Taylor polynomials for any point of expansion $x_i \in (a,b)$ of finite order $N_i$,

$$\sum_{n=0}^{N_{i}}\frac{f^{(n)}(x_{i})}{n!}(x-x_{i})^{n}.$$

I'm not a 100% on how to bound the remainder, the taylor coefficients are given recursively and hard to examine. However I know for certain that $f$ is a strictly increasing function and I "think" its derivatives behave much like the derivatives of $tan(x)$. So for a closed subinterval of $(a,b)$, the nth order derivative $f^{(n)}(x)$ attains its maximum at the endpoints. Based on this conjecture one could put a bound on Lagrange's remainder.

So with what I have I am interested in building an approximation algorithm, and was hoping to collect other peoples thoughts ... what would be a reasonable or rather a good way to go about this?

As a first step I was thinking maybe, for a given tolerance $\epsilon$, find the intervals $(a,a_0]$ and $[b_0,b)$ on which the asymptotic expansions are valid.

Then on the interval $[a_0,b_0]$ perhaps use one or more Taylor series expansions for the approximation. I say more than one because its expensive to compute the coefficients and speed is a concern for me. I'm happy enough with this basic algorithm but is there a better way. Perhaps choosing the points at which the Taylor polynomials are constructed and the degree of the polynomials in an optimal or efficient way.

Then as an alternative or extension, I was thinking maybe to combine the information from the Taylor Polynomials at the points $a_0=x_0 < x_1 < \cdots < x_M=b_0$ to construct a multipoint Pade approximation, but i'm not sure what the optimal way to do this is. How could the error be controlled?

Then of course one could economize the approximations etc. There's a host of possibilities. And I'm hoping to get some feedback from an experienced numerical analyst to save my self some time and exploring dead ends. What type of algorithms would you devise given what I have?

If it matters, I can also compute $f(x)$ to within a desired tolerance using a slow numerical scheme. But I dont mind doing this if it helps me control the error.

• Just a thought - instead of splitting your interval to use your asymptotic expansion, can't you interpolate f - f_a - f_b, where f_a and f_b are your expansions at the endpoints? That way you get a non-singular problem. In the end I guess it depends what your end goal is - do you want to integrate f? Solve a BVP that depends on f ? If your f is simple enough, I'd suggest a trivial linear interpolation, based on evaluation of f at each point of a grid. About your last remark, presumably, it isn't cheaper to compute $f^(n)(x_i)$ than $f(x_i)$, is it? – Antoine Levitt Dec 5 '11 at 22:46
• Antoine Levitt ... Sounds interesting but i'm not quite sure what you mean by interpolate $g(x)=f(x)−f_a(x)−f_b(x)$. I'm assuming $f_a$ and $f_b$ are the first few terms of the asymptotic expansions (which are divergent) restricted to a particular domain?? Unfortunately tho f is continuously differentiable, from my experiments interpolation schemes fail, near the singular points. Regarding the last question. its cheaper to compute $f^{(n)}(x_i)$ for small n. Say up to 30. then I guess it would be cheaper to compute $f(x_i)$. Why do you ask? – mark Dec 5 '11 at 23:45
• I didn't understand what you wrote about the singularity. Assume $a = 0$, $f(x) = 1/x + x + o(x)$ at 0. Then, define $g(x) = f(x) - 1/x$. $g$ isn't singular anymore, so you can interpolate it, and reconstruct f from it. I ask about the cost because I'm not convinced it's worth the bother to go higher up in the approximation scheme. Have you tried a simple linear interpolation? Why is it unsuitable? You haven't answered the question of your end goal. – Antoine Levitt Dec 6 '11 at 6:36
• I see what you mean now. This is a good idea but unfortunately the singularity cannot be factored out so easily. Yes I have tried linear interpolation, it is slow and fails near the singular points. The end goal is to approximate $f$ for use in a simulation – mark Dec 6 '11 at 16:54
• Why don't you just tell what your $f$ is? It'll allow us to see your particular problem instead of giving general advices that may easily miss the point. – fedja Dec 7 '11 at 1:05

## 1 Answer

First of all, if you need a polynomial approximation, I'd rather rescale $f$ to put it on $(-1,1)$ and interpolate it at the Chebyshev nodes http://en.wikipedia.org/wiki/Chebyshev_nodes. Furthermore, instead of using the monomial basis, as for Taylor polynomials, it is better to expand the function in a series of Chebyshev polynomials, and truncate it at a certain value. Next, multipoint Padé approximants (rational interpolants) work well at the Chebyshev nodes too. Finally, if you can allow for some algebraic functions as interpolants, you can try to approximate first a function of the form $h= f g$, where $g$ is a known function such that $h$ has no singularities at the end points.

Btw, you might want to check out the Chebfun project http://www2.maths.ox.ac.uk/chebfun/ where several of the ideas mentioned here are nicely implemented.

• I dont know the explicit form of $f$ so I'm not sure if I can generate the coefficients in the Chebyshev expansions. Do you know of any references describing the construction of multipoint pade approximations? – mark Dec 6 '11 at 16:47
• Actually, Chebyshev-Pade approximation could do better in your context (these are rational functions $r$ whose Chebyshev series should match that of $f$ as far as possible). I think that Ricardo Pachon has implemented many of these algorithms in Matlab, you might want to check his webpage, maths.ox.ac.uk/node/10864 Take a look especially at his recent publications listed there. – Andrei MF Dec 7 '11 at 0:41
• OP says that he has the inverse error function; since it's quite easy to generate series coefficients for that (via, say, Lagrangian inversion), he should have no trouble constructing Padé approximants. – J. M. isn't a mathematician Dec 8 '11 at 2:25
• @Andrei MF. Thanks. This is a learning exercise for me so I want to do as much of the coding myself.To construct the coefficients of the Chebyshev series expansions and in turn the Chebyshev-Pade rational approximation I must use FFT which requires me to evaluate $f$ at the chebyshev nodes. I can evaluate $f$ at these points through slow root finding. So I was wondering if I can avoid this and instead economize the Taylor series to obtain an expansion in chebyshev polynomials (of course like FFT the coefficients will only be approximate). Is this a viable approach? or a stupid idea? – mark Dec 8 '11 at 7:09
• Sorry, I probably don't get your question right, but Chebyshev nodes are explicit, see e.g. en.wikipedia.org/wiki/Chebyshev_nodes So, don't see why you need the root finding. – Andrei MF Dec 8 '11 at 10:49