From an earlier post, I learned that a closed form approximating expression for the series sum of a harmonic progression can be defined as $$\sum_{k=1}^{f}\frac{1}{1+ak}=\frac{1}{a} (H_{f+ \frac{1}{a} } - H_{\frac{1}{a} } )$$ where $H_{f+ \frac{1}{a}}$ and $H_{\frac{1}{a} }$ are the approximating harmonic numbers and defined as $$H_{\frac{1}{a} } = \gamma +log( \frac{1}{a} )+\frac{1}{2 \frac{1}{a}}-\frac{1}{12( \frac{1}{a}) ^{2}}+\frac{1}{120 ( \frac{1}{a})^{4} }+O \big( \frac{1}{( \frac{1}{a})^{5}} \big)$$ and $$H_{f+ \frac{1}{a} } = \gamma +log(f+ \frac{1}{a} )+\frac{1}{2(f+ \frac{1}{a} )}-\frac{1}{12(f+ \frac{1}{a} )^{2} }+\frac{1}{120(f+ \frac{1}{a} )^{4} }+O \big( \frac{1}{(f+ \frac{1}{a} )^{5}} \big)$$
When I numerically test this approximation over a range of positive values for $f$ and $a$, where $a \ll 1$, the results do not tend toward a progressively more accurate solution as the number of Bernoulli terms increase. Rather, after the first two Bernoulli terms, the solution starts to diverge towards a progressively less accurate result. I don't understand why this is so.
On a whim, I began to test different numerical constants in the Bernoulli terms (i.e., 12, 120, etc. above) and arrived at a seemingly arbitrary number set (i.e., 12,172,2084) that converges to a result that is two orders of magnitude more accurate (i.e. 0.03% vs 6%). The expressions I used to approximate the harmonic numbers are;
$$H_{\frac{1}{a} } = \gamma +log( \frac{1}{a} )+\frac{1}{2 \frac{1}{a}}-\frac{1}{12( \frac{1}{a}) ^{2}}+\frac{1}{172 ( \frac{1}{a})^{4} }-\frac{1}{2084 ( \frac{1}{a})^{6} }$$ and $$H_{f+ \frac{1}{a} } = \gamma +log(f+ \frac{1}{a} )+\frac{1}{2(f+ \frac{1}{a} )}-\frac{1}{12(f+ \frac{1}{a} )^{2} }+\frac{1}{172(f+ \frac{1}{a} )^{4} }-\frac{1}{2084(f+ \frac{1}{a} )^{6} }$$
Can someone please explain to me (1) why the generally accepted approximating expression produces a low accuracy diverging result and (2) what is the significance of my arbitrary solution that performs much better?