Message Boards Message Boards

1
|
6315 Views
|
8 Replies
|
9 Total Likes
View groups...
Share
Share this post:

Record breaking direct summation of MRB constant

POSTED BY: Marvin Ray Burns
8 Replies
POSTED BY: Henrik Schachner

Dear Henrik,

that's a beautiful piece of mathemagic.

Thanks for sharing,

Marco

POSTED BY: Marco Thiel

Dear Marco,

thank you for your nice and encouraging comments! Yes, it is almost magical - and it shows a sort of "beauty of simplicity"! I learned about Shanks transformation when reading - or trying to read - Benders book (which for me is somewhat challenging) and watching his lectures (which I find excellent!). Bender shows lots of most interesting methods I had never heard of before, particularly not in my math courses (about 30 years ago ...).

Best wishes -- Henrik

POSTED BY: Henrik Schachner

Henrik et al,

I have a big job ahead of me if I try to use Shanks to break my record of calculated digits of the MRB constant!! The terms needed to break my record of 3,014,991 digits seems impractical to list! Maybe I just need a new way of computing the Shanks' terms (someway of using known acceleration methods to accumulate them).

POSTED BY: Marvin Ray Burns

Dear Marvin,

as a preliminary remark I need to say that I have not followed this obviously long ongoing discussion closely, because I am not really interested in this kind of calculations. So forgive me please if my idea was posted long time ago or if my post is completely off topic.

It is obvious from the definition that your MRB constant emerges from an alternating series. Any partial sum has an equivalent behavior (in the following I am using not tens of thousands of terms but just the first 100 terms - as a proof of concept):

(* first 100 terms of the series: *)    
m100 = Table[N[(-1)^k (k^(1/k) - 1), 200], {k, 1, 100}];
(* partial sums of the series *)
am100 = Accumulate@m100;
ListPlot[am100, GridLines -> Automatic]

enter image description here

This is a typical situation where a Shanks transformation seems promising:

shanks[ac_List] := 
 Table[(ac[[n + 1]] ac[[n - 1]] - ac[[n]]^2)/(ac[[n + 1]] + ac[[n - 1]] - 2 ac[[n]]), {n, 2, Length[ac] - 1}]

which already leads to a huge improvement:

enter image description here

The nice thing about Shanks transformation is that it can be applied iteratively, leading to further improvements, e.g.:

sac = NestList[shanks, am100, 7];

The difference between the last and the second last term can be regarded as a measure of the error. What we finally get is:

enter image description here

Doing these calculations is just a matter of milliseconds.

Finally I can not resist in citing the great Carl Bender in his very first lecture on mathematical physics (at about min 43):

If you need to find the sum of a series, the worst possible technique is to add the numbers of the series together one at a time ...

Regards -- Henrik

POSTED BY: Henrik Schachner

Thank you Henrik! I think the Shanks transformation will help me. It seems, however, Mathematica has an issue with 10's and 100's of iterations of Shanks.

POSTED BY: Marvin Ray Burns

Let g[n] be the error in using f, as you use 10^n terms, as shown in the original post; using the following piece of code:

m = NSum[(-1)^k (k^(1/k) - 1), {k, 1, Infinity}, 
   WorkingPrecision -> 1000, Method -> "AlternatingSigns"];

N::meprec: Internal precision limit $MaxExtraPrecision = 50.` reached while evaluating -E^(2 I Interval[{0,\[Pi]}])+E^(2 I Interval[{0,\[Pi]}]). >>

f[k_] := Log[(k - Log[k])/k];

ClearAll[g]; 
g[n_] := Module[{}, 
  m - (NSum[(-1)^k (k^(1/k) - 1 + f[k]), {k, 1, 10^n}, 
      WorkingPrecision -> Floor[3 n + 100], 
      Method -> "AlternatingSigns"] - 
     NSum[(-1)^k*f[k], {k, 1, Infinity}, 
      WorkingPrecision -> Floor[3 n + 100], 
      Method -> "AlternatingSigns"])]

. You find there is a pattern to g[n]

For[n = 2, n < 20, 
 Print[RootApproximant[
   MantissaExponent[g[10 n]/g[10 (n + 1)]][[1]]]]; n++]

8/27

27/64

64/125

125/216

216/343

343/512

512/729

729/1000

1000/1331

1331/1728

1728/2197

2197/2744

2744/3375

3375/4096

4096/4913

4913/5832

5832/6859

6859/8000

, and find g[10 n]/g[10 (n + 1)] ->n^3/(n+1)^3 *10 to a factor:

For[n = 2, n < 20, Print[
  RootApproximant[
    MantissaExponent[
      g[10 n]/g[10 (n + 1)]
      ][[1]]
    ] - n^3/(n + 1)^3
  ]; n++]

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0
POSTED BY: Marvin Ray Burns

Again. let g[n] be the error in using f, as you use 10^n terms, as shown in the original post; using the following piece of code, and then it appears to be the case that the g[10 n]'s are rational factors of Log[10]^3, to better and better precision as n gets large:

 m = 
  NSum[(-1)^k (k^(1/k) - 1), {k, 1, Infinity}, 
   WorkingPrecision -> 1000, Method -> "AlternatingSigns"];



 f[k_] := Log[(k - Log[k])/k];

ClearAll[g];
g[n_] := Module[{}, 
  m - (NSum[(-1)^k (k^(1/k) - 1 + f[k]), {k, 1, 10^n}, 
      WorkingPrecision -> Floor[3 n + 100], 
      Method -> "AlternatingSigns"] - 
     NSum[(-1)^k*f[k], {k, 1, Infinity}, 
      WorkingPrecision -> Floor[3 n + 100], 
      Method -> "AlternatingSigns"])]



 Table[
 RootApproximant[MantissaExponent[Log[10]^3/g[10 n]][[1]]], {n, 2, 20}]

(* {3/20, 4/9, 3/16, 24/25, 5/9, 120/343, 15/64, 40/243, 3/25, \ 1200/1331, 25/36, 1200/2197, 150/343, 16/45, 75/256, 1200/4913, \ 50/243, 1200/6859, 3/20} *)

POSTED BY: Marvin Ray Burns
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract