# Approximating power series coefficients — Why does a clearly illegitimate method (sometimes) work so well?

For reasons that don't matter here, I want to estimate the power series coefficients $$t_{ij}$$ for the rational function $$T(x,y)= {(1+x)(1+y)\over 1- x y(2+x+y+x y)}=\sum_{i,j} t_{ij}x^iy^j$$

Using a method that I cannot justify, I get highly accurate estimates when $$i=j$$ and highly inaccurate estimates when $$|i-j|$$ strays at all far from zero.

My questions are:

Q1) Why does my apparently illegitimate method work so well when $$i=j$$?
Q2) Why does the answer to Q1) not apply when $$i\neq j$$ ?
(Of course, once the answer to Q1) is known, the answer to Q2) might be self-evident.)

I'll first present the method, then explain why I think it shouldn't work, then present the evidence that it works anyway when $$i=j$$, and then present the evidence that it rapidly goes haywire when $$i\neq j$$.

The Apparently Illegitimate Method:

Note that $$t_{ij}=t_{ji}$$, so we can limit ourselves to estimating $$t_{j+k,j}$$ for $$k\ge 0$$.

I) Define $$T_k(y)=\sum_jt_{k+j,j}y^j$$ For example, a residue calculation gives

$$T_0(y)= {1-y-\sqrt{1-4y+2y^2+y^4}\over y\sqrt{1-4y+2y^2+y^4}}$$

It turns out that all of the $$T_k$$ share a branch point at $$\zeta\approx .2956$$ and are analytic in the disc $$r<\zeta$$.

II) Write $$L_k=\lim_{y\mapsto \zeta} T_k(y)\sqrt{y-\zeta}$$. Discover that $$L_0\approx 1.44641$$ and $$L_k=L_0/\zeta^{k/2}$$.

III) Approximate $$T_k(y)\approx L_k/\sqrt{y-\zeta}$$

IV) Expand the right hand side in a power series around $$y=0$$ and equate coefficients to get $$t_{ij}\approx \pm{L_0\over\sqrt{\zeta}}\pmatrix{-1/2\cr j\cr}\zeta^{-(i+j)/2} \approx \pm 2.66036 \pmatrix{-1/2\cr j\cr}\zeta^{-(i+j)/2}\qquad(E1)$$

Remarks:

1. Obviously one could try to improve this approximation at Step III by using more terms in the power series for $$T_k$$ at $$y=\zeta$$. This doesn't seem to help, except when $$k=0$$, in which case the original approximation is already quite good.

2. For $$k\ge 2$$, $$T_k(y)$$ has a zero of order $$k-1$$ at the origin. Thus one could modify this method by approximating $$T_k(y)/(y^{k-1})$$ instead of $$T_k(y)$$ This yields $$t_{ij}\approx \pm{2.66036}\pmatrix{-1/2\cr 1-i+2j}\zeta^{-(i+j)/2}\qquad(E2)$$ (E2) is (much) better than (E1) in the range $$i\ge 2j+1$$, where it gets exactly the correct value, namely zero. Otherwise, it seems neither systematically better nor worse.

Why Nothing Like This Should Work: The expansion of $$T_k(y)$$ at $$\zeta$$ contains nonzero terms of the form $$A_{i,j}(\zeta-y)^j$$ for all positive integers $$j$$. (I'm writing $$i=j+k$$ to match up with the earlier indexing.) The truncation at Step III throws all these terms away. Therefore the expansion around the origin in Step IV ignores (among other things) the contribution of $$A_{ij}$$ to the estimate for $$t_{ij}$$. So unless we can control the sizes of the $$A_{ij}$$, we have absolutely no control over the quality of the estimate.

And in fact, even when $$k=0$$, the $$A_{j,j}$$ are not small. For example, $$t_{8,8}=8323$$ and my estimate for $$t_{8,8}$$ is a respectable $$8962.52$$. But $$A_{8,8}$$, which should have contributed to that estimate and got truncated away, is equal to $$58035$$. It seems remarkable that I can throw away multiple terms of that size and have the effects nearly cancel. I'd like a conceptual explanation for this.

But When $$i=j$$, It Works Anyway: and these get even better if you truncate just slightly farther out.

Why any explanation can't be too general: The paper 'A New Method for Computing Asymptotics of Diagonal Coefficients of Multivariate Generating Functions' by A. Raichev and M. Wilson has the precise machinery that can solve this problem. Get a copy and these brief notes correspond to their symbols, for the diagonal case $$f_{n,n} = [x^n \, y^n ] \frac{I(x)}{J(x)} = [x^n \, y^n ] \frac{(1+x)(1+y)}{1-xy(2+x+y+xy)}.$$ Solve the simultaneous system for zeros of the denominator $$J$$; in Mathematica, Solve[ { x D[J,x] == y D[J,y], J==0 },{x,y} ]. The proper solution must have both $$x$$ and $$y$$ positive. That set is $$\mathbf c=(\rho,\rho),\,\rho=(\tau-2/\tau-1)/3, \, \tau=(17+3\sqrt{33})^{1/3} \approx 0.543689.$$ This solution set, with identical $$c_1 = c_2$$, falls under the purview of a simplified calculation, in which it can be shown $$f_{n,n} \sim \rho^{-2n} \frac{1.5009481}{\sqrt{n}}.$$ (The true amplitude can be written in terms of $$\rho$$, and as many decimal places as wanted are possible, but I'm not going to bother to typeset it.)

For comparison:

1. $$n=40$$, $$\text{true}=3.4601\times 10^{20}$$, $$\text{asym} = 3.5261\times 10^{20}$$ , $$\text{absolute % err} = 1.91\%$$.

2. $$n=200$$, $$\text{true}=7.6554\times 10^{104}$$, $$\text{asym} = 7.6847\times 10^{104}$$ , $$\text{absolute % err} = 0.38\%$$.

For the non-diagonal case, you will be looking at $$f_{an,bn}$$. The machinery should work, though it is more complicated. You'll get roots that depend on $$(a,b)$$ and have to solve a complicated determinant to get the amplitude, also dependent on $$(a,b)$$. The question is, why does the non-diagonal case deviate rapidly from the diagonal? Is there a way to understand this qualitatively? I believe the answer lies in the modified form $$f_{an,bn} \sim c_1(a,b)^{-a n} c_2(a,b)^{-b n} \cdot \operatorname{amp}(a,b)/\sqrt{n}.$$ The amplitude will vary only like a polynomial upon changing $$(a,b)$$, but the first two factors have an exponential dependence.

• MathJax likes to typeset code in backticks , so I changed from the TeXed version of your Mathematica code to backticked code. – LSpice Jul 1 at 15:10
• @LSpice Thanks for the edits. I've been wanting to know how to embed links and your edit showed me how. – skbmoore Jul 1 at 15:57
• Glad it helped! I like to edit by hand, but the GUI also has some nice tools. (The link button looks vaguely like two links of chain, if you squint.) There's also some editing help, which is accessible by clicking on the '?' icon when you post. (Personally I disagree with the advice that [link] followed by : url is easier to read than [link](url)`, but to each their own. Posts aren't meant to be read in source form anyway.) – LSpice Jul 1 at 17:06
• For n=40 and n=200, my method gives better approximations than yours does. My first question is why my method does so well. I don't think you've addressed that. – Steven Landsburg Aug 11 at 5:38
• @StevenLandsburg True, I haven't answered 'why my method does so well,' or that for certain cases yours might be better. I can guess that for certain $n$ a perturbed $\rho$ could give a better result? By the way, my $\rho,$ given as an explicit algebraic number, appears to be the square of your $\zeta,$ given only numerically. You also have some binomials, which may reduce to that factor $1/\sqrt{n}.$ The answer I presentedgotten in my answer. – skbmoore Aug 12 at 14:46