This notebook contains additional material to the Mathematica GuideBooks  (Springer-Verlag, New York, 2004) http://www.MathematicaGuideBooks.org . It uses the stylesheet of the Mathematica GuideBooks, which is embedded.

Symbolics Volume, Chapter 3

Find a high-order series approximation for the solution of (∑_(k=1)^n)1/k=m for n=n(m) for large m.

Normal[Series[FunctionExpand[HarmonicNumber[n]], 
       {n, Infinity, 0}]] // Solve[# == m, n]&

A first order solution will contain the product log function.

Normal[Series[FunctionExpand[HarmonicNumber[n]], 
       {n, Infinity, 1}]] // Solve[# == m, n]&

To obtain a high order series solution, we look for a Laurent series in ε=exp(m-γ). The series must have the form n∝ε-1/2+(Subscript[c, 1])/ε+(Subscript[c, 2])/(ε^2)+…. It is straightforward to calculate the first 100 (Subscript[c, k]) recursively.

Module[{k = 0, nS = ε - 1/2, δnS, εnS, c, sol, o = 100},
 Do[k++; 
    (* add new series term *)
    δnS = Normal[Series[EulerGamma + 1/n + PolyGamma[0, n], 
                {n, Infinity, 2}]] /.  n -> nS + c[k]/ε^k;
    (* calculate series *)
    εnS = Series[δnS, {ε, Infinity, k + 1}] - (EulerGamma - Log[1/ε]);
    (* cancel series term *)
    coeff = Coefficient[Normal[εnS], ε, -(k + 1)];
    sol = Solve[coeff == 0, c[k]];
    nS = nS + c[k]/ε^k /. sol[[1]], {o}];
    (* resulting series *)
    n𝒮[ε_] = nS];

Here are some of the terms of the resulting series.

n𝒮[ε] // Take[#, -16]&

The left graphic shows the size of the coefficients of the series in ε. The right graphic shows the values of the series for m=1.

maxExp = Max[-Last /@ Cases[n𝒮[ε], ε^_, Infinity]];
n𝒮List[ε_] = MapIndexed[#1 ε^(2 - #2[[1]])&, 
                Reverse[CoefficientList[ε^maxExp #& /@ n𝒮[ε], ε]]];

Off[Graphics::gptn];
Show[GraphicsArray[
Block[{$DisplayFunction = Identity},
      (* show base 10 logarithm *)
      {ListPlot[Log[10, Abs[n𝒮List[ε] /. ε -> 1 // N]]],
       ListPlot[Log[10, Abs[n𝒮List[ε] /. (ε -> Exp[m - EulerGamma] /. 
                                                 m -> N[1, 20])]]]}]]]

The difference (Subscript[H, n(m)])-m shows the quality of the series. Due to the argument exp(m-γ), the quality increases exponentially with m.

δm[m_] := HarmonicNumber[Plus @@ n𝒮List[Exp[m - EulerGamma]]] - m

$MaxExtraPrecision = 100;
ListPlot[Table[{m, Log[10, Abs[δm[m] // N[#, 20]&]]}, {m, 36}]]