# Long waiting time for computing a summation

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}]]

-
I'd like to see what happens when I set $1000$ instead of $15$. – A_math_ninja Oct 6 '13 at 18:57

## 3 Answers

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.

-
Thank you for your answer (+1). It seems that $N$ in front caused some trouble. – A_math_ninja Oct 6 '13 at 19:19
yeah, the limit is precisely $3^{1/3}$. – A_math_ninja Oct 6 '13 at 19:32
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! – A_math_ninja 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}]


-

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.

-