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.

Finding curcuit resistance

One of the advanced physics problems which had made an early impression on me, was to find a resistance of an n-steps ladder made out of identical resistors, and to solve this problem for an infinite steps ladder.
Solving of this problem in Mathematica makes it easy to deeper explore involved mathematical concepts, and these include continued fractions, algebraic numbers, recurrence equations to name a few. For instance, the resistance of an n-steps ladder is
Despite appearance of algebraic numbers in the answer, resistance is a rational number for every natural integer n. First few of these read 3 R, 11 R/4, 41 R/15 and so on.

The solution can be downloaded as a Mathematica Player notebook. One would need freely available Mathematica Player to view it and evaluate it. Enjoy!

Divisibility by 7

Some divisibility criteria are being taught in elementary school, some in middle school. The following is the summary of what a high schooler ought to know.

1. A number is divisible by 2, 5, or 10 if its last digit is divisible by, correspondingly, 2, 5, or 10.
2. A number is divisible by 3, or 9 if the total of its digits is divisible by, correspondingly, 3, or 9.
3. A number is divisible by 4 if its last 2 digits are divisible by 4.
4. A number is divisible by 8 if its last 3 digits are divisible by 8.
   

One might wonder when is the number divisible by 7 ? To answer this question it is useful to remember some properties of modular arithmetic. These are

1. (a*b) mod n = (a mod n)*(b mod n) mod n
2. (a+b) mod n = (a mod n) + (b mod n) mod n

Thus for n = d[0] + d[1]*10^1 + d[2]*10^2 + ... + d[k]*10^k,

n mod 7 = sum( d[k] (10^k mod 7), over k )

Quick check in Mathematica shows that 10^k mod 7 is periodic with period 6:

In[30]:= Table[Mod[PowerMod[10, k, 7], 7, -3], {k, 0, 12}]

Out[30]= {1, 3, 2, -1, -3, -2, 1, 3, 2, -1, -3, -2, 1}

This gives the divisibility by 7 test. The number is divisible by 7 if the cyclic sum of its digits (starting with rightmost one) with {1,3,2,-1,-3,-2} is divisible by 7. The following is a Mathematica code to compute this cyclic sum:

DivisibilityBySevenTest[n_Integer] :=
With[{dig = Mod[Reverse[IntegerDigits[n]],7]},
   dig.PadRight[{},Length[dig], {1,3,2,-1,-3,-2}]
]

Consider an easy example. Let n = 28. Since 8 mod 7 = 1, we have (8 mod 7)*1+2*3 = 7, and, of course, 28 = 7*4.

In[47]:= DivisibilityBySevenTest[28]

Out[47]= 7

Now take n=2009. We do (9 mod 7)*1 + 0*3+0*2+2*(-1) = 2-2 = 0.

In[48]:= DivisibilityBySevenTest[2009]

Out[48]= 0

Indeed:

In[46]:= 2009/7

Out[46]= 287

When manually applying the divisibility test for large numbers, it might be easier to use the fact that 1000^k mod 7 is (-1)^(k-1). Thus, we might first form an alternating sum of triples of digits, and test the result afterwards. For example, let n = 123,456,789. It is divisible by 7 if 789 - 456+123 = 333+123 = 456 is divisible by 7. To check this, we compute 6*1+5*3+4*2 = 6+15+8=29 = 1 mod 7. Thus the number n is not divisible by 7. Indeed,

In[50]:= QuotientRemainder[123456789, 7]

Out[50]= {17636684, 1}

Applying the divisibility by 7 directly we get

In[53]:= DivisibilityBySevenTest[123456789]

Out[53]= -13

Wolfram Demonstration project

The ease with which Mathematica language allows to express ideas is most impressively seen from the success of the Demonstration project. It features over 2000, and growing, freely available demonstrations. All it takes to explore this right away, is free Mathematica Player.

One cool demonstration I just saw was a demonstration of Einstein's proof of Pythagorean theorem by John Kiehl. The heart of the proof is in quadratic scaling of triangle's area with the linear size the triangle.

It is very easy to become the publisher too. You need ideas, Mathematica 6 and an account.

Brackets for zeros of just plotted function

Given a complicated function one often waits significant amount of time for the function to plot. Once plotted you might want to determine its zeros, or other constant values saving oneself a call to

FindRoot

. The following function extract zero brackets from the given plot:

AllZeroInThePlot[gr_Graphics, val_:0] :=
Catch[Module[{a, pos, x, k},
  a = First[gr]; pos = Position[a, Line[_List]];
  If[pos==={}, Throw[$Failed]];
  a = Extract[a, pos]; k = 0;
  k = (Reap[Scan[Function[ln,
     x = ln /. Line -> Identity;
     x = Split[x, Sign[Last[#1]-val] == Sign[Last[#2]-val] &];
     If[Length[x] > 1,
         x = ((Sequence @@ Through[{First, Last}[#]]) & /@ x);
         x = Flatten[First /@ x];
         x = Prepend[x, First[x]];
         x = Append[x, Last[x]];
         x = Rest[Most[Partition[x, 2]]];
     Sow[x, k]; k++;];],a]]//Last);       
 (Flatten[#, 1] & /@ k) /. {x1_, x2_} /; x1 == x2 :> {x1}
]]

Then

In[64]:=
gr=Plot[{AiryAi[-x],AiryBi[-x]},{x,0,5}];

In[65]:=
zeros =AllZeroInThePlot[gr]

Out[65]=
{ { {2.28795,2.50089}, {3.98017,4.19035}},
 { {1.04426,1.24898}, {3.12924,3.33915}, {4.79223,5.}}
}

Or, use it to find where function attains a particular value:

In[67]:=
gr2=Plot[Sin[x]/x,{x,0,5}];

In[68]:=
AllZeroInThePlot[gr2,0.5]

Out[68]=
{{{1.88289,2.09816}}}

ExcudedForms for Simplify

We shall define

Simplify2[e_, opts__]

admitting an option

ExcludedForms

, just like

FullSimplify

does.

The idea behind, is to wrap parts matching excluded patterns around with

Hold

, do the simplification and restore held parts back to normal.

In[38]:=
Options[Simplify2] =
  Union[Options[Simplify], {ExcludedForms :> {}}];

In[39]:=
Simplify2[e_, (opts___)?OptionQ] :=
Catch[Module[{popts, tag, ef, p},
 popts = FilterOptions[Simplify, opts];
 ef = ExcludedForms/.{opts}/.Options[Simplify2];
 If[ef === {}, Throw[Simplify[e, popts]]];
 ef = Alternatives @@ ef;
 Simplify[
    e /. RuleDelayed @@ {p:ef, Hold[tag[p]]}
 , popts] /. {Hold[tag[p_]] :> p}
]]

In[40]:=
Simplify2[Sin[x]^4 + Cos[x]^4 + Sin[x]^2 + Cos[x]^2]

Out[40]=
(1/4)*(7 + Cos[4*x])

In[41]:=
Simplify2[Sin[x]^4 + Cos[x]^4 + Sin[x]^2 + Cos[x]^2,
  ExcludedForms -> {(_Sin | _Cos)^4}]

Out[41]=
1 + Cos[x]^4 + Sin[x]^4

Enjoy!

Magic Squares

Magic square is a square of integer numbers such that sum
of numbers in each row, each column and two diagonals
are the same.

Generating magic squares is easy with Mathematica. For instance,
solving for 4 by 4 matrices

vars = Variables[matrix = Array[a, {4, 4}]];

eqs = Equal @@ Join[  
 Plus @@@ matrix,                      (* sum of rows *)
 Plus @@@ Transpose[matrix],           (* of columns *)
 {Tr[matrix ], Tr[Reverse /@ matrix] } (* of diagonals *)
];

sol =
matrix /. ToRules[And @@ 
     Cases[Reduce[eqs , vars, Integers], _Equal]]

we get an answer in terms of free 8 integers

{{n1, n2, n3, n4}, 
 {n5, n6, n7, n1 + n2 + n3 + n4 - n5 - n6 - n7}, 
 {n8, n1 - n4 + n5 - n7 + n8, 
  n2 + n3 + 2 n4 - n5 - n6 - n8, n6 + n7 - n8}, 
 {n2 + n3 + n4 - n5 - n8, n3 + 2 n4 - n5 - n6 + 
   n7 - n8, n1 - n3 - n4 + n5 + n6 - n7 + n8, 
  -n4 + n5 + n8}
}