Mathematica programming

Sunday, January 20, 2008

Selberg integral

Selberg integral is a multi-dimensional integral involving . A nice review can be found in ArXiv article by P. Forrester and S. Warnaar.

An example of Selberg integral was considered by Mehta:

Attempting to verify this beautiful formula using Mathematica for small values of n:
In[1] := Table[Timing[Integrate[
Product[((x[i]-x[j])^2),{i,n},{j,i+1,n}]*
Exp[-n Sum[x[i]^2,{i,n}]],
Sequence @@ Table[{x[i], -Infinity, Infinity}, {i, n}]]
], {n, 2, 5}]

Out[1] = {{0.86, Pi/4}, {5.703, Pi^(3/2)/(54 Sqrt[3])},
{20.875, (9 Pi^2)/131072},
{147.078, (27 Pi^(5/2))/(195312500 Sqrt[5])}}

In[2]:= Table[(2Pi)^(n/2)/(2n)^(n^2/2) Product[k!,{k,n}],{n,2,5}]

Out[2]= {Pi/4, Pi^(3/2)/(54 Sqrt[3]), (9 Pi^2)/131072,
(27 Pi^(5/2))/(195312500 Sqrt[5])}
In order to achieve more efficient evaluation, and check the formula for higher values of n, we note that Integrate[x^(m-1) Exp[-n x^2],{x,-Infinity,Infinity}]==Sin[Pi m/2]^2 Gamma[m/2]/n^(m/2). We can then integrate variable one by one, carefully expanding only when needed with Expand[expr,patt]:
In[1] := Table[Timing[(Pi/n)^(n/2)*
Fold[Expand[#1, #2] //. {#2^(m_) :>
Cos[Pi*m/2]^2*Pochhammer[1/2, m/2]/n^(m/2)} /. {#2->0} &,
Product[(x[i]-x[j])^2, {i,n}, {j,i+1,n}],
Table[x[k], {k, n}]]], {n, 2, 7}]

Out[1] = {{0., Pi/4},
{0., Pi^(3/2)/(54 Sqrt[3])},
{0.015, (9 Pi^2)/131072},
{0.125, (27 Pi^(5/2))/(195312500 Sqrt[5])},
{2.016, (25 Pi^3)/3343537668096},
{34.594, (273375 Pi^(7/2))/(875799914882589322976 Sqrt[7])}}

In[2]:= %[[All, 2]] ==
Table[(2 Pi)^(n/2)/(2 n)^(n^2/2) *
Array[Factorial, n, 1, Times], {n,2, 7}]

Out[2]= True
My 2GB RAM machine runs out of memory for n=8 for both methods, but
the speed-up of the latter method is clear.

0 Comments:

Post a Comment

<< Home