[Math] Continued fraction expansion related to exponential generating function

algorithmsbernoulli numberscontinued-fractionspower seriesspecial functions

A recent SciComp.SE Question motivates us to ask for a nice continued fraction expansion of the following Maclaurin series:

$$ f(x) = \sum_{n=0}^\infty \frac{B_n\; x^{n+3}}{n! (n+3)} = \int_0^x \frac{t^3}{e^t – 1} dt $$

where $B_n$ are the (first) Bernoulli numbers.

Although the positive power of $t$ in the numerator of the integrand removes the singularity at the origin, the exponential function has complex period $2\pi i$. Therefore the radius of convergence for the above series is $2\pi$, the distance from the origin to the nearest singularities.

Where power series convergence is always limited by the proximity of singularities, continued fraction expansions typically have larger domains of convergence, perhaps including the whole real line. Such a circumstance weighs in favor of rational over polynomial approximations over a widened real interval, and nice continued fractions often give apt starting values for a minimax rational approximation.

Best Answer

The method of Viskovatov

While smooth real functions can have at most one convergent power series expansion in the neighborhood of a point (Taylor series), expanding the same function as a continued fraction presents many options. A basic technique (Viskovatov, 1806) operates on the ratio of a pair of power series (or polynomials) to produce a Thiele-type continued fraction or with slight variation, a corresponding-type continued fraction:

$$ p(x)/q(x) = c_0 + \cfrac{x^{\alpha_1}}{c_1 + \cfrac{x^{\alpha_2}}{c_2 + \cfrac{x^{\alpha_3}}{c_3 + \cdots}}} $$

where it is assumed $q(0)$ is nonzero.

Specifically, after subtracting off the constant term $c_0 = p(0)/q(0)$:

$$ p(x)/q(x) = c_0 + \frac{p(x) - c_0 q(x)}{q(x)} $$

If the latter numerator is not identically zero, then $p(x) - c_0 q(x) = x^{\alpha_1} r(x)$ since it has a root at $x=0$ with multiplicity $\alpha_1 \ge 1$. Thus in that case we expand:

$$ p(x)/q(x) = c_0 + \frac{x^{\alpha_1}}{q(x)/r(x)} $$

By construction $r(0)$ is nonzero, so the method applies to $q(x)/r(x)$ and may be continued until potentially an identically zero numerator is reached at some stage. Terminating fractions are then essentially rational expressions in $x$.

Truncating such an expansion, omitting the ellipsis at some stage, gives rational expressions analogous to nested multiplication forms of polynomials in that divisions replace the multiplies, especially in the simple case where exponents $\alpha_i$ are always one:

c0 + x/(c1 + x/(c2 + x/(c3 + ...)))

Some regularity can be seen in the Thiele fraction expansion around $x=0$ of the function $f(x)/x^3 = \sum_{k=0}^{\infty} \frac{B_k x^k}{k!(k+3)}$:

$$ \frac{f(x)}{x^3} = \frac{1}{3} + \cfrac{x}{-8 + \cfrac{x}{-\frac{15}{16} + \cfrac{x}{8 + \cfrac{x}{\frac{7}{5} + \cfrac{x}{-8 + \cfrac{x}{-\frac{27}{16} + \cfrac{x}{8 + \cfrac{x}{\frac{275}{161} + \cdots}}}}}}}} $$

The curious reader will find a few more of the constants:

$$ c_9 = -8, c_{10} = -\frac{144417}{97264} $$ $$ c_{11} = 8, c_{12} = \frac{36954241}{25202043} $$ $$ c_{13} = -8, c_{14} = -\frac{1763915461119}{1034996631248} $$ $$ c_{15} = 8, c_{16} = \frac{53785750138088275}{33572832373135209} $$ $$ c_{17} = -8, c_{18} = -\frac{559934841935054771682651}{394380619248992328976208} $$ $$ c_{19} = 8, c_{20} = \frac{123408586919942234000562495331463}{74004329695604288642520453854805} $$

seem to dispose of easy conjectures about the even index constants $c_{2k}$ while also confirming the alternating appearances of $\pm 8$ in odd indices. For the profoundly curious reader I have an Appendix below that describes my implementation in Maxima CAS of the method of Viskovatov.

For the sake of comparison with the Maclaurin approximations given in the Background section, I've chosen Thiele-type expansions of orders 4, 12, and 20 (as multiplied by removed factor $x^3$) to match up with polynomial degrees 7, 15, and 23 previously shown. In one sense this is a conservative choice, as the higher odd index Bernoulli coefficients being zero would allow those same Maclaurin polynomials to be construed as approximations of degree 8, 16, and 24.

3rd order Debye function and Thiele fractions

Fig. 2 Debye function of 3rd order and its Thiele-type fractions

The plot demonstrates that increasing the order of expansion does not "hit a wall" in the way the power series did, that the domain of approximation appears to extend smoothly with increasing orders.

Appendix: Maxima CAS implementation

I always anticipate barking my shins at every turn when trying to perform symbolic calculations with a CAS, and this experience with Maxima was in line with that. Maxima has no datatype for continued fractions, but it does provide Taylor series constructions and a function to compute Pade approximations. The latter was useful as a check on the results of my Viskovatov method.

This is a user-defined Maxima function that returns a rational expression:

/*
  Maxima function that performs Thiele fraction expansion by
  applying Viskovatov's method to the ratio of two power series.

  viskovatovRat(Variable,PowerSeries0,PowerSeries1,Depth)

    returning a rational expression

  Expands about Variable=0; assumes Depth is nonnegative integer.
  Assumes value of PowerSeries1 at Variable=0 is nonzero.
*/

viskovatovRat(x,p0,p1,n) :=
block ( [a,c,p2],
    c : ev(p0,[x=0])/ev(p1,[x=0]),
    if n <= 0 
       then return(c),
    p2 : ratsimp(p0 - c*p1),
    a : 0,
    if p2 = 0
       then return(c),
    while ev(p2,[x=0]) = 0
    do (
          a : a+1,
          p2 : ratsimp(p2/x)
    ),
    c + (x^a/viskovatovRat(x,p1,p2,n-1))
)$

I also created a variant that returns a list of constants plus powers of $x$ which was convenient for compactness of display.

Related Question