Asymptotic solution of the saddle-point equation

Question

As an example we have the following equation:

$$\sum _{j=1}^{\infty } \frac{r^j}{\left(1-r^j\right)^2}=n$$

 Sum[r^j/(1 - r^j)^2, {j, 1, Infinity}] == n

I'm looking for a solution for unknown r. My question is if it possible to find an asymptotic solution of this equation with Mathematica ?

The first iteration (for j = 1) is

 Solve[r/(1-r)^2 == n, r]

 (* {{r -> -((-1 - 2 n + Sqrt[1 + 4 n])/(2 n))}, {r -> (1 + 2 n + Sqrt[1 + 4 n])/(2 n)}} *)

The first of two solutions gives a minimum r

 (* N[% /. n -> 1000] *)
 (* {{r -> 0.968873}, {r -> 1.03213}} *)

Algebraic expressions for next roots are complicated, but here is a numerical approximation:

 (* N[Solve[r/(1 - r)^2 + r^2/(1 - r^2)^2 == n, r] /. n -> 1000] *)
 (* {{r -> -1.01593}, {r -> -0.984316}, {r -> 0.965265}, {r -> 1.03598}} *)

The root with the smallest absolute value is always at position 3. The following program finds an asymptotics iterations from 2 - 5 terms.

 (* Do[so = Solve[(Sum[r^j/(1 - r^j)^2, {j, 1, terms}]) == n, r, Reals]; Quiet[rasy = Table[Expand[FullSimplify[Normal[Series[r /. so[[root]], {n, Infinity, 1}]], n > 0]], {root, 1, Length[so]}]]; Print[rasy[[3]]];, {terms, 2, 5}] *)

 (* 1 - Sqrt[5]/(2 Sqrt[n]) + 5/(8 n) *)
 (* 1 - 7/(6 Sqrt[n]) + 49/(72 n) *)
 (* 1 - Sqrt[205]/(12 Sqrt[n]) + 205/(288 n) *)
 (* 1 - Sqrt[5269]/(60 Sqrt[n]) + 5269/(7200 n) *)

$$1-\frac{\sqrt{5}}{2 \sqrt{n}}+\frac{5}{8 n}$$

$$1-\frac{7}{6 \sqrt{n}}+\frac{49}{72 n}$$

$$1-\frac{\sqrt{205}}{12 \sqrt{n}}+\frac{205}{288 n}$$

$$1-\frac{\sqrt{5269}}{60 \sqrt{n}}+\frac{5269}{7200 n}$$

The final result (the limit of these iterations) is:

$$r \sim 1-\frac{\pi }{\sqrt{6 n}}+\frac{\pi ^2}{12 n}$$

But how to find it with Mathematica ? Generalized methods are welcome.

The numerical check:

 (* With[{n = 1000000000}, N[Sum[r^j/(-1 + r^j)^2, {j, 1, Infinity}]/n /. r -> 1 - π/Sqrt[6 n] + π^2/(12 n), 10]] *)
 (* 0.9999876725 *)

UPDATE

One method how to find this limit is guess a formula for the coefficients. We have

$$r_k=1-\frac{c_k}{\sqrt{n}}+\frac{c_k^2}{2 n}$$

 Simplify[FindSequenceFunction[{1/2, 5/8, 49/72, 205/288, 5269/7200, 5369/7200, 266681/352800, 1077749/1411200}, k]]
 (* 1/12 (π^2 - 6 PolyGamma[1, 1 + k]) *)

 Limit[%, k -> Infinity]
 (* π^2/12 *)

$$\frac{c_k^2}{2}=\frac{\pi ^2}{12}$$

But such proof is not rigorous...

Answer 1

Mathematica can compute the closed form solution of your sum:

s = Sum[r^k/(1-r^k)^2, {k, Infinity}, Assumptions->0<r<1]

QPolyGamma[1, 1, r]/Log[r]^2

So, the natural thing to do is to find the series expansion of the above expression, but it seems that Mathematica has a Series bug for QPolyGamma:

Series[QPolyGamma[1, 1, r], {r, 1, 0}]

QPolyGamma[1, 1, 1]

$\frac{1}{6} \left(\pi ^2-3\right)+O\left((r-1)^1\right)$

$π^2/6$

Note that the leading order term is off by -1/2. Looking in the literature, we find the following asymptotic expansion for QPolyGamma near $r=1$:

QPolyGamma[1, 1, r] ∼ π^2/6 - Sum[Zeta[-k]/k! BernoulliB[k,1] Log[r]^(k+1), {k, 0, Infinity}]

$\psi _r^{(1)}(1)\sim \frac{\pi ^2}{6}-\sum _{k=0}^{\infty } \frac{B_k(1) \zeta (-k) \log ^{k+1}(r)}{k!}$

It turns out that Zeta[-2 k] is 0 for $k\geq 1$, and BernoulliB[2k+1, 1] is 0 for $k>1$. Hence, only the first 2 terms contribute, and we find that:

qapprox[r_] = π^2/6 - Sum[Zeta[-k]/k! BernoulliB[k,1] Log[r]^(k+1), {k, 0, 2}]

π^2/6 + Log[r]/2 + Log[r]^2/24

Let's compare qapprox with QPolyGamma for $r\lesssim 1$ (the function qapprox is missing a term needed when $r\gtrsim 1$):

Table[N@{qapprox[r], QPolyGamma[1, 1, r]}, {r, 1 - 2^-Range[10]}]

{{1.31838, 1.31838}, {1.50454, 1.50454}, {1.57891, 1.57891}, {1.61284, 1.61284}, {1.6291, 1.6291}, {1.63707, 1.63707}, {1.64102, 1.64102}, {1.64298, 1.64298}, {1.64396, 1.64396}, {1.64445, 1.64445}}

The function qapprox is simple enough that Solve can invert it:

sol = r /. Simplify[Solve[qapprox[r]/Log[r]^2==n && n>10], n>10 && r<1]

{E^((6 - 2 Sqrt[9 + (-1 + 24 n) π^2])/(-1 + 24 n)), E^(( 2 (3 + Sqrt[9 + (-1 + 24 n) π^2]))/(-1 + 24 n))}

Choosing the solution corresponding to $r<1$:

sol = First @ DeleteCases[sol, x_ /; Simplify[x>1, n>10]]

E^((6 - 2 Sqrt[9 + (-1 + 24 n) π^2])/(-1 + 24 n))

If you still want the asymptotic expansion:

Series[sol, {n, Infinity, 2}] //Simplify //TeXForm

$1-\frac{\pi \sqrt{\frac{1}{n}}}{\sqrt{6}}+\frac{3+\pi ^2}{12 n}-\frac{\left(27+39 \pi ^2+4 \pi ^4\right) \left(\frac{1}{n}\right)^{3/2}}{144 \left(\sqrt{6} \pi \right)}+\frac{63+21 \pi ^2+\pi ^4}{864 n^2}+O\left(\left(\frac{1}{n}\right)^{5/2}\right)$

Finally, here's a check:

With[{r = sol /. n -> N@{10, 100, 1000}}, Sum[r^k/(1-r^k)^2,{k,Infinity}]]

{10., 100., 1000.}

Answer 2

You can solve your problem in principal using InverseSeries:

nr = Simplify[ Series[Sum[r^j/(1 - r^j)^2, {j, 1, 5}], {r, 1, 3}]]
InverseSeries[nr, n ]
(* 1 + 1/60 Sqrt[5269] Sqrt[1/n] + 5269/(7200 n) +... *)

but I don't know how to generalize to infinite sum...

appendix

the infinite sum evaluates to

nr\[Infinity]=Simplify[ Series[Sum[r^j/(1 - r^j)^2, {j, 1, Infinity}], {r,1, 1}]]

thereby the summation of the first to parts can be calculated

Sum[1/j^2, {j, 1, Infinity}]
(* Pi^2/6 *)
nr\[Infinity]=Pi^2/6 (1/(r-1)^2+1/(r-1)) +...

Unfortunately the part O[(r-1)^0 ] doesn't converge...

solution Creating a seriesdata object with the previous result the inverse series evaluates to

InverseSeries[SeriesData[r, 1, {Pi^2/6, Pi^2/6, const}, -2, 0, 1], n]
(*SeriesData[n,DirectedInfinity[], {1, 6^Rational[-1, 2] Pi, Rational[1, 12] Pi^2}, 0, 3, 2]*)