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.

Symbolics Volume, Chapter 1

Calculate the value of the integral

(Subscript[ℐ, n])=(∫_0^1)(∫_0^1)…(∫_0^1)(∫_0^1)(e^(|1-(Subscript[x, n])|))(e^(|(Subscript[x, n])-(Subscript[x, n-1])|))… (e^(|(Subscript[x, 3])-(Subscript[x, 3])|)) (e^(|(Subscript[x, 2])-(Subscript[x, 1])|))(e^(|(Subscript[x, 1])-0|))d(Subscript[x, n])d(Subscript[x, n-1])… d(Subscript[x, 2])d(Subscript[x, 1])

for 2≤n≤10. Is it possible to calculate (Subscript[ℐ, 100])?

integrand[n_] := Product[Exp[Abs[x[k] - x[k - 1]]], {k, 1, n + 1}] /.
                                         {x[0] -> 0, x[n + 1] -> 1};
                                         
iterators[n_] := Table[{x[k], 0, 1}, {k, n}];

Here is the resulting integral for n=4.

integrand[4]

(HoldForm[Integrate[##]]& @@ {integrand[6], Sequence @@ iterators[6]}) /.
                               x[k_] :> Subscript[x, k] // TraditionalForm

The built-in Integrate can do such integrals, but they are quite slow.

ReleaseHold[%] // Timing

The time required to calculate these integrals grows quite fast with n, so obtaining a value for n=10 seems not possible through a direct call to Integrate. So, we will implement our own integration routine for this special integral. We will split each absolute value into two parts that do not contain absolutes values anymore (the first and last exponents need no splitting, they have constant sign over the integration region [0,1]). This gives us (2^(n-1)) n-dimensional integrals.

(* thread over sums *)
removeAbs[p_Plus, k_] := removeAbs[#, k]& /@ p

(* split integration region into two for each variable *)
removeAbs[integrate[body_, iterators_], k_] :=
With[{(* first and last exponent need only one region *)
      body1 = body /. {Abs[x[1]] -> x[1], 
              Abs[1 - x[Length[iterators]]] -> 1 - x[Length[iterators]]}},
 integrate[body1 /. Abs[x[k] - x[k - 1]] :> x[k - 1] - x[k], 
           MapAt[{x[k], 0, x[k - 1]}&, iterators, k]] +
 integrate[body1 /. Abs[x[k] - x[k - 1]] :> x[k] - x[k - 1], 
           MapAt[{x[k], x[k - 1], 1}&, iterators, k]] ]

splitIntegrals[n_] := Fold[removeAbs, integrate[integrand[n], iterators[n]], 
                           Range[2, n]] // ExpandAll;

Here are the eight abs-free integrals for n=4.

splitIntegrals[4]

We convert the multiple integrals into explicit iterated integrals.

splitIntegrals[6][[4]] /. integrate[integrand_, iterators_] :> 
                          FoldList[integrateInd, integrand, Reverse[iterators]]

To carry out a single integral of the multivariate integral, we define a function integrateInd. Because the integrands are smooth and free of singularities, we can just evaluate the indefinite integral (which is fast) and substitute limits.

integrateInd[f_, {x_, x0_, x1_}] := 
     If[(* integrand does not contain integration variable? *)
        FreeQ[f, x, Infinity], (x1 - x0) f, 
        (* avoid size swell by applying Together *)
        Together[((# /. x -> x1) - (# /. x -> x0))&[Integrate[f, x]]]]

integrateInd[f_, {x_, x0_, x1_}] := 
     If[(* integrand does not contain integration variable? *)
        FreeQ[f, x, Infinity], Together[(x1 - x0) f], 
        (* avoid size swell by applying Together *)
        Together[((# /. x -> x1) - (# /. x -> x0))&[Integrate[f, x]]]]

Now we have all ingredients together to define the function absMultiIntegral that calculates the original integral in dimension n.

absMultiIntegral[n_] := 
 splitIntegrals[n] /. integrate[integrand_, iterators_] :> 
     Fold[integrateInd, integrand, Reverse[iterators]] // Together

Now we can carry out the integrals much faster (although still with exponentially in n growing times).

Table[{n, absMultiIntegral[n] // Timing}, {n, 2, 10}]

We can achive a further speed-up by recognizing that we can write (Subscript[ℐ, n])=(Subscript[(OverTilde[ℐ]), n])(1), where (Subscript[(OverTilde[ℐ]), n])(1) obeys the recursion relation (Subscript[(OverTilde[ℐ]), n])(x) =(∫_0^1)exp(|x-y|)(Subscript[(OverTilde[ℐ]), n-1])(y) dy. This allows to do only two integrals per integration variable. The next function intR implements this recursive approach. At each stage, we keep (Subscript[(OverTilde[ℐ]), n])(x) fully expanded.

(* save dynamic values for intR[n][x_] *)
Off[RuleDelayed::rhs];

intR[n_][x_] := intR[n][x_] = (* thread integration over sums and expand *)
ExpandAll[Exp[ x] ((# /. y -> x) - (# /. y -> 0))&[
             Integrate[# Exp[-y], y]& /@  intR[n - 1][y]] +
          Exp[-x] ((# /. y -> 1) - (# /. y -> x))&[
             Integrate[# Exp[+y], y]& /@  intR[n - 1][y]]]

(* initial value *)
intR[1][x_] := Exp[2 - x]/2 - Exp[x]/2 + x Exp[x]

The recursive approach allows to calculate (Subscript[ℐ, 10]) in about one second.

(* do recursion step by step *)
(Do[intR[n][x], {n, 2, 10}]; Together[intR[10][1]]) // Timing

And in a approximately an hour it is now possible to calculate (Subscript[ℐ, 100]). Assumping that the calculation time t depends on n in the form c (n^α), we can estimate the time needed to calculate (Subscript[ℐ, 100])by extrapolating from the quickly computable first 30 values of (Subscript[ℐ, n]).

(* calculation times as a function of n *)
timeData = Table[{n, Timing[intR[n][x]][[1, 1]]}, {n, 11, 30}];
(* parameter fitting an assumed timing of the form t[α] = Exp[c] n^α *)
{c, α} = CoefficientList[Fit[Log @ timeData, {1, n}, n], n];
Sum[Exp[c] n^α, {n, 100}] Second

This means that in approximately an hour it is now possible to calculate (Subscript[ℐ, 100]).

Do[(* to see some progress *)
   Print[{k, Timing[Length[intR[k][x]]]}], {k, 31, 100}] // Timing;
   
intR[100][1]

Here are the leading terms in e of the result for (Subscript[ℐ, 100])=2.2825… (10^16):
 1/2^100 E^101 + 97/2^100 E^99+ 3533/(2^98 3) E^97 + 103745/(2^96 3^2) E^95 +
 240415611/(2^100 5 7) E^93 + …