Conversion from Taylor- to continued fraction expansion

This routine converts a vector  of Taylor coefficients into an (equally long) vector of continued fraction coefficients. A function is then given, which evaluates this expansion. The numerical result is compared to the exact values in case of a Stieltjes' function.

Enter a test case (Stieltjes' function)

Remove["Global`*"]

\(TraditionalForm
    \`Enter\ the\ input\ vector\ of\ Taylor\ coefficients.\ Here\ we\ have\ 
        choosen\ a\_i = \ \((\(-1\))\)\^i\ \(i!\), \ 
    corresponding\ to\ the\ Stieltjes'\ function\ \ 
      \(∫\_0\%∞\( e\^\(-\ t\)\/\(1 + z\ t\)\) 
          dt\ \ 
          \(\((not\ strictly\ a\ Taylor\ series, \ 
              since\ it\ has\ zero\ radius\ of\ convergence)\) : \)\)\)

a = {1, -1, 2, -6, 24, -120, 720}

Algorithm to get continued fraction coefficients

Determine the highest power series coefficient that is entered.

n = Length[a] - 1

Next,  
    we need to initialize the first columns (superscript zero) for
       the variables (Subscript[A, i]^((0))) and   (Subscript[B, i]^((0))) 
      (in the notation of B& O;  
        in our notation  A[i, 0] and B[i, 0] resp.)

Do [  A[i, 0] = a[[i + 1]];  B[i, 0] = 0, {i, 0, n}]; 
    B[0, 0] = 1;

Here comes the main loop; equations (8.4.4a) and (8.4.4b) in B&O:

\(Do[\ \ \ \ \ 
      Do[
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        A[i, j] = B[i + 1, j - 1] - A[i + 1, j - 1]/A[0, j - 1]; 

        			\ \ \ \ \ \ \ \ \ \ \ 
        B[i, j] = \ A[i, j - 1]/A[0, j - 1]\ \ \ \ , {i, 0, n - j}]\ \ 

      \ \ \ \ \ , {j, 1, n}]\)

We now pick off the continued fraction coefficients from the top row in the A-array

c = Table[A[0, j], {j, 0, n}]

Function that evaluates a continued fraction at a point

The input arguments to the function cfval[c_,z_] are
c   - a vector with continued fraction coefficients, for example as generated above,
z   - independent variable

\(cfval[c_, z_] := 
	\ \ 
      \((len = Length[c] - 1; r[\(-1\)] = 0; \ s[\(-1\)] = 1; 
        r[0] = c[\([1]\)]; s[0] = 1; 
		\ 
        Do\ [\ 	r[k] = r[k - 1] + c[\([k + 1]\)]*z*r[k - 2]\ ; \ 

          			\ \ \ \ \ \ \ \ \ \ 
          s[k] = s[k - 1] + c[\([k + 1]\)]*z*s[k - 2]\ , 
          \ \ \ \ \ {k, 1, len}]; 
	\ \ r[len]/s[len]\ )\)\)

Test of evaluation function

We now plot the imaginary part of the Stieltjes' function, and compare it with what we get from the continued fraction approximation

z = x + I y ;

Stieltjes' function:

Plot3D[
      Im[NIntegrate[Exp[-t]/(1 + z t), {t, 0, ∞}]], {x, 
        -1, 1}, {y, -1, 1}, PlotPoints -> 20]

Continued fraction approximation:

Plot3D[Im[cfval[c, z]], {x, -1, 1}, {y, -1, 1}, 
      PlotPoints -> 20]

We can also compare along the positive real axis. Direct power series expansion (dotted) and continued fraction expansion (dashed) are compared tothe exact result (solid curve):

\(taylor = 
      Plot[Sum[a[\([k + 1]\)]\ z^k, {k, 0, Length[a] - 1}], {z, 0, 3}, 

        			\ \ \ PlotRange -> {0, 2}, PlotStyle -> Dashing[{0.01, 0.01}], 
        DisplayFunction -> Identity]; 

    cfract = Plot[cfval[c, z], {z, 0, 3}, PlotRange -> {0, 2}, 
        PlotStyle -> Dashing[{0.08, 0.08}], DisplayFunction -> Identity]; 

    exact = Plot[
        NIntegrate[Exp[\(-t\)]/\((1 + z*t)\), {t, 0, Infinity}], {z, 0, 3}, 
        PlotRange -> {0, 2}, DisplayFunction -> Identity]; 

    Show[{taylor, cfract, exact}, DisplayFunction -> $DisplayFunction]\)

We should finally note that increasing the number of power series coefficients that we started with will rapidly make the direct evaluation result worse, but the continued fraction result even better.