# Is there a numerical algorithm for finding an asymptotic slope?

I have a series of data points $(x_i,y_i)$ which I expect to (approximately) follow a function $y(x)$ that asymptotes to a line at large $x$. Essentially, $f(x) \equiv y(x) - (ax + b)$ approaches zero as $x \to \infty$, and the same can probably be said of all the derivatives $f'(x)$, $f''(x)$, etc. But I don't know what the functional form for $f(x)$ is, if it even has one that can be described in terms of elementary functions.

My goal is to get the best possible estimate of the asymptotic slope $a$. The obvious crude method is to pick out the last few data points and do a linear regression, but of course this will be inaccurate if $f(x)$ does not become "flat enough" within the range of $x$ for which I have data. The obvious less-crude method is to assume that $f(x) \approx \exp(-x)$ (or some other particular functional form) and fit to that using all the data, but the simple functions I've tried like $\exp(-x)$ or $\dfrac1{x}$ don't quite match the data at lower $x$ where $f(x)$ is large. Is there a known algorithm for determining the asymptotic slope that would do better, or that could provide a value for the slope along with a confidence interval, given my lack of knowledge of exactly how the data approach the asymptote?

This sort of task tends to come up frequently in my work with various data sets, so I'm mostly interested in general solutions, but by request I'm linking to the particular data set that prompted this question. As described in comments, the Wynn $\epsilon$ algorithm gives a value that, as far as I can tell, is somewhat off. Here is a plot:

(It does look like there's a slight downward curve at high x values, but the theoretical model for this data predicts that it should be asymptotically linear.)

• This may well be too elementary - or too vague - for this site, but I figured private beta is the time to try such things out. Commented Dec 6, 2011 at 8:53
• No, I think this is a great question. Not everything must be advanced and fancy. Good solutions to simple problems are important. Commented Dec 6, 2011 at 16:42
• @Dan: was replacing $\exp$ really justified? Commented Dec 6, 2011 at 19:14
• Having exps tends to make things harder for me to read, but I'll admit that it was small enough that I shouldn't have done it.
– Dan
Commented Dec 6, 2011 at 19:20
• I really don't care either way, I just figured I might as well approve the edits because, well, why not. You get a couple reputation off it, whatever that's worth. Commented Dec 6, 2011 at 19:33

It's a rather rough algorithm, but I'd use the following procedure for a crude estimate: if, as you say, the purported $f(x)$ that represents your $(x_i,y_i)$ is already almost linear as $x$ increases, what I'd do is to take differences $\dfrac{y_{i+1}-y_i}{x_{i+1}-x_i}$, and then use an extrapolation algorithm like the Shanks transformation to estimate the limit of the differences. The result is hopefully a good estimate of this asymptotic slope.

What follows is a Mathematica demonstration. The Wynn $\epsilon$ algorithm is a convenient implementation of the Shanks transformation, and it is built in as the (hidden) function SequenceLimit[]. We try out the procedure on the function

$$\frac4{x^2+3}+2 x+e^{-4 x}+3$$

xdata = RandomReal[{20, 40}, 25];
ydata = Table[(3 + 13*E^(4*x) + 6*E^(4*x)*x + x^2 + 3*E^(4*x)*x^2 +
2*E^(4*x)*x^3)/(E^(4*x)*(3 + x^2)), {x, xdata}];

SequenceLimit[Differences[ydata]/Differences[xdata],
Method -> {"WynnEpsilon", Degree -> 2}]
1.999998


I might as well show off how simple the algorithm is:

wynnEpsilon[seq_?VectorQ] :=
Module[{n = Length[seq], ep, res, v, w}, res = {};
Do[ep[k] = seq[[k]];
w = 0;
Do[v = w; w = ep[j];
ep[j] =
v + (If[Abs[ep[j + 1] - w] > 10^-(Precision[w]), ep[j + 1] - w,
10^-(Precision[w])])^-1;, {j, k - 1, 1, -1}];
res = {res, ep[If[OddQ[k], 1, 2]]};, {k, n}];
Flatten[res]]

Last[wynnEpsilon[Differences[ydata]/Differences[xdata]]]
1.99966


This implementation is adapted from Weniger's paper.

• Just curious, but why did you the original form of the function, instead of combining all of the terms? Commented Dec 6, 2011 at 13:37
• It was only for demonstration purposes. :) I included the $\LaTeX$-ed expression so you guys know what the expected answer was supposed to be. What I wanted to demonstrate was that you can do this sort of analysis on some complicated-looking function... Commented Dec 6, 2011 at 13:40
• How close to flat do the points have to be for the algorithm to be effective? Commented Dec 6, 2011 at 14:01
• Okay, last question (I swear), can you generate an error bound on the estimate? Commented Dec 6, 2011 at 18:20
• That's a bit trickier. I've seen some suggested methods in a few papers, but I confess to not having done experiments with them. (Maybe I should, one of these days.) The book by Brezinski and Redivo-Zaglia has a few pointers you might want to look into. Commented Dec 6, 2011 at 18:24