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

In exercise 35 of Chapter 1 of the Symbolics volume, the average distance of two randomly chosen points from within the unit square was calculated. Now, calculate the probability distribution 𝓅(d) for the occurring lengths 0≤d≤(2^(1/2 )).

cf = Compile[{{n, _Integer}, {bins, _Integer}},
            Module[{tab = Table[0, {bins}], λ, b},
                   Do[λ = Sqrt[(Random[] - Random[])^2 + 
                               (Random[] - Random[])^2];
                      b = Ceiling[bins λ/Sqrt[2]];
                      tab[[b]] = tab[[b]] + 1, {n}];
                    tab/n]];

The resulting distribution is relatively broad with a maximum around d≈0.5.

lcp = With[{bins = 100, n = 10^8}, 
ListPlot[MapIndexed[{Sqrt[2] #2[[1]]/bins, bins/Sqrt[2] #1}&, 
                    cf[n, bins]], PlotStyle -> PointSize[0.006]]];

The exact distribution is given by the following integral

𝓅(λ)=(∫_0^1)(∫_0^1)(∫_0^1)(∫_( 0)^1)δ(λ-(Sqrt[((((Subscript[x, 2])-(Subscript[x, 1])))^2)+((((Subscript[y, 2])-(Subscript[y, 1])))^2)])) d(Subscript[y, 2])d(Subscript[y, 1])d(Subscript[x, 2])d(Subscript[x, 1])

Mathematica cannot carry out even the innermost integration symbolically.

𝒹[{x1_, y1_}, {x2_, y2_}] := Sqrt[(x1 - x2)^2 + (y1 - y2)^2]

Integrate[DiracDelta[λ - 𝒹[{x1, y1}, {x2, y2}]], {y1, 0, 1 - x1}]

So we will use the well known formula (∫_ℐ)δ(y-f(x))dx=(∫_ℐ)(∑_k)δ(y-(Subscript[x, k]))/|f'((Subscript[x, k]))|dx where the (Subscript[x, k]) are the zeros of f(x) in the interval ℐ. For the integrand under consideration, we have up to two zeros.

sol = Solve[λ - 𝒹[{x1, y1}, {x2, y2}] == 0, y2] // Simplify

The denominator absolute value is identical for the two roots.

{int1y2, int2y2} = (1/D[λ - 𝒹[{x1, y1}, {x2, y2}], y2] /. # // 
                                  Simplify[#, λ > 0]&)& /@ sol

inty2 = int1y2;

For carrying out the various integration, we subdivide the (Subscript[x, 1])-(Subscript[x, 2])-(Subscript[y, 1]) space into regions of identical sign of the root expression (this automatically ensures the existence of the roots).

cond1 = 0 < x1 < 1 && 0 < x2 < 1 && 0 < y1 < 1 &&
        0 < sol[[1, 1, 2]] < 1 && λ > 0;

cad1 = CylindricalDecomposition[cond1, {λ, x1, x2, y1}];

We group the resulting regions with respect to the λ-domains.

expandCAD[cad_] := 
{#[[1, 1]], Rest /@ #}& /@ Split[List @@ #, #1[[1]] == #2[[1]]&]& @ 
  DeleteCases[cad  //. a_ && (b_ || c_) :> (a && b) || (a && c) /. 
                                                 LessEqual -> Less,
              _?(MemberQ[#, _Equal, Infinity]&)]

Here is the subdivision corresponding to the first root.

cad1R = expandCAD[cad1]

And here is the subdivision corresponding to the second root.

cond2 = 0 < x1 < 1 && 0 < x2 < 1 && 0 < y1 < 1 &&
        0 < sol[[2, 1, 2]] < 1 && λ > 0;

cad2 = CylindricalDecomposition[cond2, {λ, x1, x2, y1}];
cad2R = expandCAD[cad2]

We have four λ-intervals.

{First /@ cad1R, First /@ cad2R}

λIntervals = %[[1]];

Now, we must carry out the remaining three integrations for the four λ-domains. The function x1x2y1Integration iteratively forms the indefinite integral and then calculates the values at the upper and lower limit to obtain the value of the definite integral. We do the definite integration manually to avoid that Mathematica spends time on checking the integration range for singularities and discontinuities; because the integration domains were constructed to ensure a smooth dependence in their inside, this is not needed.

x1x2y1Integration[cad_, startIntegrand_] :=
Module[{𝒟ℒ, 𝓋2, x1L, x1U, 𝓋3, x2L, x2U, 𝓋4, y1L, y1U},
 (* form difference of integral at upper limit minus
    integral at lower limit *)
 𝒟ℒ = (Limit[#1, #2[[1]] -> #2[[3]], Direction -> +1] - 
       Limit[#1, #2[[1]] -> #2[[2]], Direction -> -1])&;     
 (* integration regions *)
 {{𝓋2, x1L, x1U}, {𝓋3, x2L, x2U}, {𝓋4, y1L, y1U}} = 
            cad /. Inequality[a_, Less, b_, Less, c_] :>
                                 {b, a, c} /. And -> List;
 (* carry out the three integrations *)
 Fold[𝒟ℒ[Integrate[#1, #2[[1]]], #2]&, startIntegrand, 
      {{𝓋4, y1L, y1U}, {𝓋3, x2L, x2U}, {𝓋2, x1L, x1U}}]]

Summing up the contributions for all (Subscript[x, 1])-(Subscript[x, 2])-(Subscript[y, 1])regions belonging to a given λ-domain gives the probability distribution we are looking for.

sumUp[λRange_] := 
  Simplify[Plus @@ #, λRange]& @ (x1x2y1Integration[#, inty2]& /@
   Flatten[Last /@ Cases[{cad1R, cad2R}, {λRange, _}, Infinity]])

The resulting distribution has different functional form in the two intervals [0,1] and [1,(2^(1/2))].

{𝓅[{0, 1/2}], 𝓅[{1/2, 1}], 𝓅[{1, Sqrt[5]/2}], 𝓅[{1, Sqrt[5]/2}]} = 
                                           sumUp /@ λIntervals

Here is a graphic of the resulting distribution. The earlier calculated numerical values agree excellently with the theoretical curve. The resulting curve is continuous differentiable at λ=1 (the second derivative is discontinuous).

Plot[Which[0 < λ < 1, 𝓅[{0, 1/2}], 1 < λ < Sqrt[2], 𝓅[{1, Sqrt[5]/2}]], 
     {λ, 0, Sqrt[2]}, Prolog -> {PointSize[0.009], Hue[0],
                                 Cases[lcp, _Point, Infinity]},
     PlotStyle -> {{Thickness[0.002], GrayLevel[0.8]}}];

𝓅D[λ_] := UnitStep[λ] UnitStep[1 - λ] 2 λ (Pi + λ (λ - 4)) + 
          UnitStep[λ - 1] UnitStep[Sqrt[2] - λ]*
            (2 λ (Pi - 2 - λ^2 + 4 Sqrt[λ^2 - 1] - 4 ArcTan[Sqrt[λ^2 - 1]]))

By construction, the distribution is normalized to 1.

Integrate[𝓅[{0, 1/2}], {λ, 0, 1}] + Integrate[𝓅[{1, Sqrt[5]/2}], {λ, 1, Sqrt[2]}]

The maximum occurs at (Subscript[d, max])≈0.48. Interestingly, the maximum of the distribution has nearly exactly the same distance from 0.5 as the expectation value (the maximum is slightly less than 0.5 and the maximum is slightly larger than 0.5.)

{#, N @ #}& @ Min[λ /. Solve[D[𝓅[{0, 1/2}], λ] == 0, λ]]

{(2 + Sqrt[2] + 5 ArcSinh[1])/15 - 1/2, 
 1/2 -  (4 - Sqrt[16 - 3 Pi])/3} // N

We end by calculating once more the expectation value. Using the obtained form of the distribution directly results in a complicated expression containing hypergeometric functions for the expectation value.

av1 = Integrate[λ S[{0, 1/2}], {λ, 0, 1}]

av2 = Integrate[λ 𝓅[{1, Sqrt[5]/2}], {λ, 1, Sqrt[2]}]

Calculating the last integral through the indefinite integral yields a result in elementary functions.

(* form indefinite integral *)
intλ = Integrate[λ 𝓅[{1, Sqrt[5]/2}], λ]

Forming the indefinite integral yields once more the result from solution 35a.

av2b = Limit[intλ, λ -> Sqrt[2], Direction -> 1] - 
       Limit[intλ, λ -> 1, Direction -> -1];

av1 + av2b // FullSimplify