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.

It takes a long time to compute the summation below, and I'd like to know if there are alternative ways
to compute things faster. When replacing $15$ by $\infty$, then I should get $3^{1/3}$. I need Mathematica solution
not mathematical solutions (since here I can handle things on my own).

Clear[b, k]
b[0] := 1;
b[k_] := b[k - 1] + 1/(b[k - 1]^2 + b[k - 1] + 1)
N[1/(15)^(1/3) Sum[1/(b[k]^2 + b[k] + 1), {k, 0, 15}]]
share|improve this question
    
I'd like to see what happens when I set $1000$ instead of $15$. –  Chris's sis Oct 6 '13 at 18:57

3 Answers 3

up vote 4 down vote accepted

EDIT

Clear[b, k]
b[0] := 1;
b[k_] := ( b[k] = N[b[k - 1] + 1/(b[k - 1]^2 + b[k - 1] + 1), 100] )
1/(1000)^(1/3) Sum[1/(b[k]^2 + b[k] + 1), {k, 0, 1000}]

(*1.28966188938645818444514220128278717806515296115261423714827563492762\
  4979228386408865750383532338923*)

very close to 3^(1/3) = 1.4422...

The convergence is slow. So you should do extrapolation.

share|improve this answer
    
Thank you for your answer (+1). It seems that $N$ in front caused some trouble. –  Chris's sis Oct 6 '13 at 19:19
    
yeah, the limit is precisely $3^{1/3}$. –  Chris's sis Oct 6 '13 at 19:32
1  
I actually put N in the definition of b[k], since its exact values are huge rational numbers. –  Vahagn Poghosyan Oct 6 '13 at 19:32
    
Note that here k is local variable, so it is not necessary to Clear it at the beginning of the code. –  Vahagn Poghosyan Oct 6 '13 at 19:39
    
I got the point. Thanks! –  Chris's sis Oct 6 '13 at 19:46
b[0] := 1;
b[k_] := b[k] = N[b[k - 1] + 1/(b[k - 1]^2 + b[k - 1] + 1)]
f[x_] := f[x] = f[x - 1] + (1/(b[x]^2 + b[x] + 1))
f[0] := 0

Manipulate[
 ListPlot[Table[1/x^(1/3) f[x], {x, 1, a}], MaxPlotPoints -> 100, 
  Epilog -> {Green, Line[{{1, 3^(1/3)}, {a, 3^(1/3)}}]}, 
  PlotRange -> {0, 4}], {a, 10, 100000, 1}]

Mathematica graphics

share|improve this answer

Here is an alternative,

Nest[# + 1/(1 + # + #^2) &, 1., n + 1] - 1

that will compute the sum over 0 <= k <= n. For 100-digit precision use

Nest[# + 1/(1 + # + #^2) &, 1.`100, n + 1] - 1

Examples

With[{n = 1000}, 
  1/n^(1/3) (Nest[# + 1/(1 + # + #^2) &, 1.`, n + 1] - 
     1)] // AbsoluteTiming
(* {0.000359, 1.28966} *)

With[{n = 1000}, 
  1/n^(1/3) (Nest[# + 1/(1 + # + #^2) &, 1.`100, n + 1] - 1)] // AbsoluteTiming
(* {0.005293, 
    1.28966188938645818444514220128278717806515296115261423714827563492762\
      4979228386408865750383532338923} *)

With[{n = 1000000}, 
  1/n^(1/3) (Nest[# + 1/(1 + # + #^2) &, 1.`100, n + 1] - 1)] // AbsoluteTiming
(* {5.204264,
    1.42720128317421809656439478563254431280673078216998116597726119064353\
      220814168138689903801038204268} *)

It avoids the recursion and memory limits that the OP's b runs into. (With n = 1000000, my kernel crashed using b and Sum.) And you can use NestList to get all the sums up to n.

share|improve this answer

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.