4
$\begingroup$

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}]]
$\endgroup$
1
  • $\begingroup$ I'd like to see what happens when I set $1000$ instead of $15$. $\endgroup$ Oct 6 '13 at 18:57
4
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ Thank you for your answer (+1). It seems that $N$ in front caused some trouble. $\endgroup$ Oct 6 '13 at 19:19
  • $\begingroup$ yeah, the limit is precisely $3^{1/3}$. $\endgroup$ Oct 6 '13 at 19:32
  • 1
    $\begingroup$ I actually put N in the definition of b[k], since its exact values are huge rational numbers. $\endgroup$ Oct 6 '13 at 19:32
  • $\begingroup$ Note that here k is local variable, so it is not necessary to Clear it at the beginning of the code. $\endgroup$ Oct 6 '13 at 19:39
  • $\begingroup$ I got the point. Thanks! $\endgroup$ Oct 6 '13 at 19:46
5
$\begingroup$
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

$\endgroup$
3
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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