Continued fraction expansion related to exponential generating function

Question

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.

Answer 1

Background:

It is evident from the integral form that the function $f(x)$ will have a triple root at $x=0$, that it is monotone increasing, and (due to the integrand's denominator growing exponentially) finitely bounded above (by exactly $\pi^4/15$ as it turns out).

Here are a few fairly precise function values that illustrate this trend:

$$\begin{array}{|r|c|} \hline \text{x} & \text{f(x) (shown rounded)} \\ \hline 1.0 & 0.22480518802594 \\ \hline 2.0 & 1.17634259660700 \\ \hline 3.0 & 2.55221845329080 \\ \hline 4.0 & 3.87705416153119 \\ \hline 5.0 & 4.89989215833058 \\ \hline 6.0 & 5.58585538083094 \\ \hline 7.0 & 6.00316896121307 \\ \hline 8.0 & 6.23962379489192 \\ \hline 9.0 & 6.36657389887547 \\ \hline 10.0 & 6.43192189678183 \\ \hline \end{array}$$

It was also evident in from the Question itself that the function has a power series expansion about the origin, but that the radius of convergence is only $2\pi$ owing to the nearest singularities (poles). Since these singularities "accumulate" at infinity in the complex plane, it cannot be expected that a convergent power series representation can be centered "at infinity", i.e. a Laurent series.

Nevertheless it may prove desirable to investigate that in more detail, after first exploring the possibility of converting the Maclaurin series into a continued fraction. Below the third-order Debye function being considered is graphed alongside its truncated Maclaurin series polynomials of order 7, 15, and 23 for visual consistency:

Fig. 1 Debye function of 3rd order and its Maclaurin polynomials

Notice the diminishing returns near the threshold radius of convergence, that the improvement in approximation going from order 15 to order 23 is slight in comparison with the earlier improvement from order 7 to order 15. As nonconstant polynomials inevitably veer off to plus or minus infinity, it is natural to hope for better convergence from rational expressions (which allow for the horizontal asymptote seen here).

Graphic produced using gnuplot interface from Maxima; I'll post the Maxima code as a GitHub gist once I've polished it...

Answer 2

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.

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.

Answer 3

Your integral term $$f(t)= t^3/(e^t-1)$$ can be expressed as the continued fraction $$f(t) = t^2/(1+t/(2+t/( -{3 \over 2}+t/(-8+t/({5 \over 6}+t/(18+t/(-{7 \over 12} + ... $$

$$f(t) ={ t^2 \over1+ K_{i=2}^\infty{ t \over c_i} }\ where \ c_i = \begin{cases}i\ even \ c_i=(-1)^{(i+2) \over 2}2({i \over 2})^2 \\ i \ odd \ \ \,c_i = (-1)^{(i-1) \over 2}{i \over { {i-1} \over 2}{{i+1} \over 2} } \end{cases} \ $$

Wallis' Fundamental Recurrence Formulas provide an easy way to calculate your optimal rational polynomials. $$ f(t) \approx {h_i \over k_i}$$ $$ h_i = 1 ,\ 0, \ t^2, \ \ \ \ ( 2t^2 ), \ \ (t^3-3t^2), \ \ (-6t^3 +24t^2), (t^4 -8t^2+20t^2),\ ...$$ $$ k_i = 0, \ 1, \ 1, \ \ (t+2), \ \ ({-t \over 2} -3), \ \ \ \ (6t^2 + 24), \ \ \ \ \ \ \ ({t^2 \over 3}+2t+20), \ ...$$ $$ where \ h_i = t*h_{i-2} + c_ih_{i-1} \\ and \ \ \ \ \ k_i = t*k_{i-2} + c_ik_{i-1} \\ i = 2,...,\infty$$