This notebook contains additional material to the Mathematica GuideBooks  (Springer-Verlag, New York, 2004) http://www.MathematicaGuideBooks.org . It uses the stylesheet of the Mathematica GuideBooks, which is embedded.

Numerics Volume, Chapter 2

Show through numeric integration that the following integral [✶], [✶]

(∫_0^∞)cos(2z)((∏_(k=1)^∞)cos((z/k)))dz

equals π/8 to 42 digits. For the numerical calculation of the infinite product, make use of the expansion

ln(cos(z))=(∑_(j=1)^∞)(((((-1))^j) (2^(2 j-1)) ((4^j)-1) )/(j (2 j)!))(Subscript[B, 2j])(z^(2j))

valid for |z|<π/2. Here the (Subscript[B, k]) are the Bernoulli numbers.

Module[{z = 1, cosT}, 
       cosT = Table[Cos[N[z, 22]/k], {k, 1000}];
       ListPlot[Log[10, Abs[1 - cosT]]]];

To speed up the calculation of the infinite product, we split it into two parts:

(∏_(k=1)^∞)cos((z/k))=((∏_(k=1)^(Subscript[k, 0]))cos((z/k)))((∏_(k=(Subscript[k, 0])+1)^∞)cos((z/k)))=((∏_(k=1)^(Subscript[k, 0]))cos((z/k)))exp((∑_(k=(Subscript[k, 0])+1)^∞)ln(cos((z/k))))

Here (Subscript[k, 0]) has to be chosen in such a way that z/((Subscript[k, 0])+1)<π/2 and the above expansion applies. The first product we calculate directly by carrying out the multiplications. For the tail part, we use the above expansion of ln(cos(x)).

(∑_(k=(Subscript[k, 0])+1)^∞)ln(cos((z/k)))=(∑_(k=(Subscript[k, 0])+1)^∞)(∑_(j=1)^∞)(((((-1))^j) (2^(2 j-1)) ((4^j)-1) )/(j (2 j)!))(Subscript[B, 2j]) ((((z/k)))^(2j))
=(∑_(j=1)^∞)(((((-1))^j) (2^(2 j-1)) ((4^j)-1) )/(j (2 j)!))(Subscript[B, 2j]) ζ(2j,(Subscript[k, 0])+1)(z^(2j))

Here ζ(s,a)=(∑_(k=a)^∞)(k^(-s)). This definition yields immediately the recursion ζ(s,a+1)=ζ(s,a)-(a^(-s)) starting with ζ(2j,1)=| (2^(2 j-1)) (π^(2 j)) (Subscript[B, 2 j])/(2 j)!|.

To calculate the last sum efficiently for many values of x, we write the terms of the sums a s a product of three terms: a (Subscript[k, 0])- and z-independent one, one that contains the (Subscript[k, 0])-dependent term and the one that contains the z-dependent one. Only the last has to be calculated for each value of z. The (Subscript[k, 0])-dependent terms we calculate recursively once. This yields the following implementation of a function makeListDefinitions that calculates kMax globally visible lists L2[k0] and the globally visible list L1 to precision prec.

$MaxExtraPrecision = 10000;

makeListDefinitions[o_, kMax_, prec_] :=
 (Clear[L1, L2, L2E];
  (* k0-independent terms *)
  L1 = Table[((-1)^j 2^(2j - 1) (4^j - 1) BernoulliB[2j])/(j (2j)!), 
             {j, o}] // N[#, prec]&; 
  (* calculate k0-dependent terms recursively *)
  L2E[1] = Table[(-1)^(j - 1) 2^(2j - 1) Pi^(2j) BernoulliB[2j]/(2j)!, {j, o}];      
  L2E[a_] := L2E[a] =  With[{temp = L2E[a - 1] - Table[(a - 1)^(-2j), {j, o}]}, 
                            (* save memory *) L2E[a - 1] =.; temp];  
  L2[k0_] := L2[k0] = L2E[k0 + 1] // N[#, prec]&;
   (* make all L2's in order *)
  Do[L2[k0], {k0, kMax}];)

Using the lists L1 and L2 calculated through a call to makeListDefinitions allows to calculate the infinite product through the function cos∞Product.

cos∞Product[z_?NumberQ, z0_:1] := 
Module[{k0 = Ceiling[Abs[z/z0]], p1, p2, z2 = z^2}, 
       If[z == 0, 1, 
          (* direct product of first cosines *)
          p1 = Product[Cos[z/k], {k, k0}]; 
          (* infinite tail approximation *)
          L3 = NestList[(# z2)&, z2, Length[L1] - 1]; 
          p2 = Exp[(L1 L2[k0]).L3]; 
          (* product of initial and tail terms *) p1 p2]]

Using 50 terms in the tail sum, we can reliably calculate the product for a visualization. The left graphic shows the product, and the right graphic shows the base 10 logarithm of its absolute value. We see an exponential decay for increasing z.

makeListDefinitions[50, 50, 50]

Show[GraphicsArray[
Block[{$DisplayFunction = Identity},
       {Plot[cos∞Product[z], {z, 0, 2Pi}, PlotRange -> All],
        Plot[Log[10, Abs[cos∞Product[z]]], {z, 0, 20 Pi}, 
             PlotRange -> All, PlotPoints -> 100]}]]]

Here is a view on the real part and the argument of the product over the complex plane.

Show[GraphicsArray[
Block[{$DisplayFunction = Identity},
 {(* 3D plot of the real part *)            
  Plot3D[Re[cos∞Product[x + I y]], {x, -2Pi, 2Pi}, {y, -1, 1}, 
         PlotRange -> All, Mesh -> False, PlotPoints -> 200],
  (* contour plot of the argument *)
  ContourPlot[Arg[cos∞Product[x + I y]], {x, -2Pi, 2Pi}, 
              {y, -3/2, 3/2}, PlotRange -> All,
              PlotPoints -> 200, ColorFunction -> (Hue[0.78 #]&)]}]]]

To quantify the decay of the product along the real axis we fit a straight line to the successive local maxima.

Fit[Reverse[Rest[FoldList[{#2[[1]], Max[#2[[2]], #1[[2]]]}&, 
                    {-Infinity, -Infinity}, Rest[Reverse[
       Table[{z, Log[10, Abs[cos∞Product[N[z, 50]]]]},  
             {z, 0, 20 Pi, 20 Pi/111}]]]]]], {z}, z] // N

We see that Π(x)∝(10^(-2x/3)). Using this result, we define a function zMax that for a given precision goal pg for the integral under consideration, gives an upper integration limit so that the precision goal is met.

zMax[pg_] := Ceiling[3/2 pg] + (* to be save *) 10

Now, we must determine how many terms kMax to take into account in the L1, L2 lists. When calculating the actual values of the product, we must make sure to calculate each function values with a sufficient precision. The sum in the exponent becomes less precise for increasing z.

With[{ζ = 50}, ListPlot[Log[10, N @ Abs[L1 L2[ζ] Table[ζ^(2 k), {k, 50}]]]]];

Approximately, we have (Subscript[σ, k])(z)∝(10^(-2k/5)) where (Subscript[σ, k])(z) is the kth term in the sum.

Fit[Log[10, N @ Abs[L1 L2[50] Table[50^(2 k), {k, 50}]]], {k}, k]

makeListDefinitions[75, 75, 100];
Fit[Log[10, N @ Abs[L1 L2[75] Table[75^(2 k), {k, 75}]]], {z}, z]

Given the upper integration limit zMax, we define a function kMax that yields the number of needed terms kMax.

kMax[zMax_] := Ceiling[5/2 zMax] + (* to be save *) 10

Now, we have all ingredients together to define the function integrateCosProduct that carries out the numerical integration and yields a pg-digits approximation of the integral. We use the GaussKronrod method and split the integration interval into subintervals of length δZ.

integrateCosProduct[pg_, wPrec_:Automatic, 
                    δZ_:2Pi, showFunctionEvaluations_:True] := 
Module[{Z, 𝓀, H, g2, sum = 0}, 
 If[wPrec === Automatic, wp = 3 pg, wp = wPrec];
 (* calculate needed z-range and number of terms *) 
 Z = zMax[pg]; 𝓀 = kMax[Z]; H = Ceiling[Z/δZ]; 
 (* prepare fast numerical evaluation of the product *) 
 makeListDefinitions[𝓀, Ceiling[H δZ], wp];
 (* define function to integrate *)
 g2[z_?NumberQ] := Cos[2 z] cos∞Product[z];
 (* integrate g2 in intervals *) 
 Do[int = NIntegrate[g2[z], {z, j δZ, (j + 1) δZ}, 
                      WorkingPrecision -> wp, Method -> GaussKronrod, 
                      PrecisionGoal -> pg, AccuracyGoal -> pg,
                      MaxRecursion -> 8];
    sum = sum + int, {j, 0, H - 1}];
(* return integral *) N[sum, pg]]

integrateCosProduct[50]

Interestingly, the difference to π/8 is only 7.4((…10)^(-43)).

N[Pi/8, 50]