Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I've completed a problem that involves approximating $e$ by a continued fraction:

$$\frac{N_1}{D_1+\frac{N_2}{D_2+\frac{N_3}{\ddots+\frac{N_k}{D_k}}}}$$

with the $N_i$ being the list {1,1,1, …} and the $D_i$ being the list {1,2,1,1,4,1,1,6,1,1,8, …}. (Note: The continued fraction actually approximates $e-2$).

The overarching structure of computing the continued fraction is a straightforward application of Fold with the following function defined:

f[x_, {m_, d_}] := m/(d+x);
eContFracApprox[n_] := Fold[f, 0, (*properly constructed list of pairs*)]

which shifts the crux of the problem to generating the lists for $N_i$ and $D_i$ and combining them properly.

I have a solution, but I'm not happy with it:

Reverse[
MapThread[
  List, {Table[1, {n}], 
   Take[Flatten[{Join[{1, 2}, Table[{1, 1, i}, {i, 4, n, 2}]]}],n]}
]
]

I think it's ugly and it doesn't work for values of $n < 3$. Perhaps someone can suggest a different/better approach? Specifically, I'd like to know if we can do away completely with MapThread and generate the list of pairs directly from an application of Table.


Edit 1:

Just to be clear. I'm looking for an improvement to generating the lists of the $N_i$ and $D_i$. Although, alternate ways of generating continued fractions are also appreciated.

I'll also add my complete function for reference:

eContinuedFractionApprox[n_Integer] := Module[{f},
  f[x_, {m_, d_}] := m/(d + x);
N[Fold[f, 0,
 Reverse[
  MapThread[List,
   {Table[1, {n}], Take[Flatten[{Join[{1, 2}, Table[{1, 1, i},
         {i, 4, n, 2}]]}],
     n]
    }
   ]
  ]
 ],
n
] + 2
 ]
share|improve this question
add comment

5 Answers

up vote 13 down vote accepted

As you say, this is a straightforward application of Fold, which is also perhaps the cleanest solution you can get. I'm guessing that you're seeking Table based approaches since you didn't want to deal with having to "fold properly". I'll complete the Fold application here so that you can see how it is applied:

n = {n1, n2, n3, n4, n5};
d = {d1, d2, d3, d4, d5};

f[x_, {m_, d_}] := m/(d + x);
Fold[f, Last@n/Last@d, Reverse@Most@Transpose@{n, d}]


Now coming to generating a continued fraction approximation for $e-2$, observe that if you prepend a $1$ to the $D_i$ list you showed, then we can get a nice sequence function for it (I have replaced the formal symbols in the output with normal ones for clarity):

func = FindSequenceFunction[{1, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8}]
(* DifferenceRoot[Function[{y, n}, {-2 - y[n] - y[1 + n] - y[2 + n] + y[3 + n] + 
        y[4 + n] + y[5 + n] == 0, y[1] == 1, y[2] == 1, y[3] == 2, y[4] == 1, y[5] == 1}]] *)

So if you want to use $N$ terms in your approximation, just use the above function to obtain $N+1$ entries and discard the first one.

The Fold operation now can be written easily as (for 20 terms):

eMinus2Approx[n_Integer] := With[{d = N@Rest@func@Range[n + 1], f = 1/Plus[##] &},
    Fold[f, 1/Last@d, Reverse@Most@d]
]

eMinus2Approx[20]
(* 0.718282 *)
share|improve this answer
    
The fold procedure isn't the issue. It's the list generation that I'm concerned with. I'll edit my original post to reflect that. –  Haer'Dalis Mar 24 '13 at 18:22
    
@Haer'Dalis Please see my update for a way to generate the list. –  rm -rf Mar 24 '13 at 18:30
    
That is a badass function (FindSequenceFunction). I'm continually amazed at what's been built into MMA. Thanks for that @rm-rf. –  Haer'Dalis Mar 24 '13 at 18:35
add comment

Well, the smartass answer would be

ContinuedFraction[E - 2, n]

So your approximation function (if I understood correctly) would be

FromContinuedFraction[ContinuedFraction[E - 2, n]]

But the way I would generate this particular sequence manually is something along these lines:

ls[n_] := Module[{s}, s = 2 Range[n];
                 {1, Riffle[s, {{1, 1}}]} // Flatten]
share|improve this answer
    
was not aware of ContinuedFraction either. –  Haer'Dalis Mar 24 '13 at 18:42
1  
yea usually it's a question of what Mathematica doesn't have. though i don't think i would have anticipated continued fraction functionality myself, and i don't remember how i happened upon it –  amr Mar 24 '13 at 18:46
add comment

If you are interested in the simple continued fraction representation for any irrational, not just $e-2$, then the following code generates the required list of coefficients, which are the $D_i$ in your notation. The number of coefficients requested is $n$ and the input irrational is $r$.

Map[-#[[1, 1]]&,Rest[NestList[1/(#-Floor[#])&,r,n]]]

For example, $r=\pi$ and $n=10$ returns {3, 7, 15, 1, 292, 1, 1, 1, 2, 1}.

share|improve this answer
add comment

The built-in function ContinuedFractionK[] seems to have not yet been mentioned, so:

n = {n1, n2, n3, n4, n5};
d = {d1, d2, d3, d4, d5};
ContinuedFractionK[n[[k]], d[[k]], {k, 5}]
   (d5*(d4*(d2*d3*n1 + n1*n3) + d2*n1*n4) + (d2*d3*n1 + n1*n3)*n5)/
 (d5*(d4*(d3*(d1*d2 + n2) + d1*n3) + (d1*d2 + n2)*n4) + (d3*(d1*d2 + n2) + d1*n3)*n5)

A Part::pspec error is thrown in this case even when ContinuedFractionK[] possesses the HoldAll attribute; I am not quite sure why. Nevertheless, it works, and you can use Quiet[] if you find the messages somewhat annoying.


As already mentioned, Fold[] is indeed the function you can use for evaluating continued fractions had the CF-related functions not been built-in. rm has detailed the conventional "backward evaluation" method; I will present the forward recursion scheme (formulae 16-19 here, which I also used in this answer):

Divide @@ Last[Fold[{{0, 1}, #2}.#1 &, IdentityMatrix[2], Transpose[{n, d}]]]

Usually, this method is not used in inexact arithmetic since the recursions tend to produce results that overflow, even when their ratio is perfectly representable. In Mathematica, which allows for symbolics and exact arithmetic, one has no such fear. The methods I featured in my other answer can also be adapted accordingly (e.g. Lentz-Thompson-Barnett) to this more general form of the continued fraction.


For the specific case of $e-2$:

Divide @@ Last[Fold[{{0, 1}, #2}.#1 &, IdentityMatrix[2], 
    Transpose[{ConstantArray[1, 18],
               {1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10, 1, 1, 12, 1}}]]] // N
   0.718281828459028
share|improve this answer
    
Welcome back. Added obligatory link to docs. –  rcollyer Apr 2 '13 at 14:42
    
Well, back-ish... –  J. M. Apr 2 '13 at 14:43
1  
ish, works. No big deal. missed your input, though. –  rcollyer Apr 2 '13 at 14:44
add comment

Here's another approach to manually construct the list given by ContinuedFraction[E-2, n]:

cfList[n_] := {0} ~Join~ (Flatten @ Riffle[Table[{1, 1}, {n}], 2 Range @ n] ~Take~ {2, n})

Checking the result:

FromContinuedFraction[cfList @ 10] // N
E - 2.
(*
0.718284
0.718282
*)
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.