(* This file is a short supplement to "Discriminants of some Painleve polynomials". Depending on your Mathematica system, you can either 1) read it in by <<(whatever you call this file) or 2) just cut and paste the file in with the mouse. It might be best to keep two windows open: one to look at this file, as there are comments interspersed throughout; one is your Mathematica session. Part I of this file implements the polynomials that the paper is about, namely p[m], h[m,n], and q[m,n]. Part II discusses discriminants, resultants and complex roots of these polynomials. Part II also compares the behavior of Painleve polynomials with "random polynomials" Part III is somewhat longer and I don't expect this Part to be of too much interest to others. After writing the paper, I retyped in the main formulas to check for typos. Part III shows these checks. *) (* Part I. Implementation of the polynomials. *) (* In all three cases we make use of the abbreviation b (which does not appear in the paper) *) b[j_,f_] := D[f,{x,2}] f - D[f,x]^2 + j f^2 (* For m=0,1,2,... the Yablonsky-Vorobiev polynomial P_m is given as p[m]. Note: here and later, I've kept the code short by not making any attempt to catch errors. So if you type p[-1], you'll go into an infinite recursion which might crash your session. Watch out! *) p[0] = 1 p[1] = x p[m_] := p[m] = -b[-x,p[m-1]]/p[m-2] //Factor (* For m,n=0,1,2,... the biHermite polynomial H_m,n is given as h[m,n] *) h[0,0]=1 h[1,0]=1 h[0,1]=1 h[1,1]= x h[m_,0] := h[m,0] = b[m-1,h[m-1,0]]/(h[m-2,0]*(m-1)) //Factor h[m_,1] := h[m,1] = b[m-1,h[m-1,1]]/(h[m-2,1]*(m-1)) //Factor h[m_,n_] := h[m,n] = b[1-n,h[m,n-1]]/(h[m,n-2]*(1-n)) //Factor (* For m,n=0,1,2,... the Okamoto polynomial Q_m,n is given as q1[m,n]. For m,n= 0,-1,-2,-3,..., the Okamoto polynomial Q_m,n is given as q2[m,n]. Both are then also called q[m,n]. *) q1[0,0] = 1 q1[1,0] = 1 q1[0,1] = 1 q1[1,1] = x q1[m_,0] := q1[m,0] = b[x^2+2 m - 3,q1[m-1,0]]/q1[m-2,0] //Factor q1[m_,1] := q1[m,1] = b[x^2 + 2 m-2,q1[m-1,1]]/q1[m-2,1] //Factor q1[m_,n_] := q1[m,n] = b[x^2 + 3 - m - 2 n,q1[m,n-1]]/q1[m,n-2] //Factor q2[0,0] = 1 q2[1,0] = 1 q2[0,1] = 1 q2[1,1] = x q2[m_,0] := q2[m,0] = b[x^2+2 m +1,q2[m+1,0]]/q2[m+2,0] //Factor q2[m_,1] := q2[m,1] = b[x^2 + 2 m+2,q2[m+1,1]]/q2[m+2,1] //Factor q2[m_,n_] := q2[m,n] = b[x^2 -1 - m - 2 n,q2[m,n+1]]/q2[m,n+2] //Factor q[m_,n_] := Which[ m>=0 && n >= 0,q1[m,n], m<=0 && n<=0,q2[m,n]] (* Part II. Discriminants, resultants and complex roots. Comparison of Painleve polynomials and random polynomials. *) (* randpol[d,c] gives a monic degree d polynomial with remaining coefficients integers randomly chosen from {-c,-c+1,...,c-1,c} *) randpol[d_,c_] := x^d + Sum[Random[Integer,{-c,c}] x^k,{k,0,d-1}] (* If f is a nonzero polynomial in x then disc[f] is its discriminant. disc[randpol[10,10]] //FactorInteger should illustrate what happens with a "randomly chosen" polynomial: at least one large prime divides the discriminant. On the other hand, disc[p[10]] //FactorInteger illustrates what happens with the Painleve polynomials here: all primes dividingthe discriminants are small. *) disc[f_] := If[Exponent[f,x]==0,1, If[Mod[Exponent[f, x], 4] <= 1, Resultant[f, D[f, x], x], -Resultant[f, D[f, x], x]]] (* If f and g are polynomials in x then res[f,g] gives their resultant. res[randpol[10,10],randpol[10,10]] should illustrate that typically at least one large prime divides resultants. However, the resultants treated in the paper (and a few others in the biHermite and Okamoto cases) are unusual in that they are products of small primes only . Examples like res[p[4],p[8]] represent intermediate cases. They factor into small primes raised to large powers times large primes raised to small powers. *) res[f_,g_] := Resultant[f,g,x] (* If f is a polynomial in x then rootplot[f] graphs its roots in the complex plane. rootplot is designed to work at least somewhat on systems with "terminal graphics". So matter what your system, you should see a marked difference between e.g. rootplot[h[3,4]] and rootplot[randpol[12]] . To get pictures like Figure 6.1 of the paper, one has to do a bit more formatting. *) rootplot[f_] := Module[{lis}, lis = x /. NSolve[f==0,x]; ListPlot[Transpose[{Re[lis], Im[lis]}], Axes -> False, Frame -> True, FrameTicks->False]] (* Part III. Here are the checks. Any command that ends D (for "difference") should return 0 for any allowable input. Commands that end Q (for "quotient") should return 1 or -1 for allowable inputs. *) (* First we check that the rational functions defined in Section 1 really do satisfy Painleve 2 or Painleve 4 *) scalef[m_] := scalef[m] = (D[Log[p[m+1]/p[m]],x] /. x->2^(-2/3) x) 2^(-2/3) //Factor pdiffD[m_] := D[scalef[m],{x,2}] - (2 scalef[m]^3 + x scalef[m] -m-1) //Factor scaleh[m_,n_] := (D[Log[h[m,n+1]/h[m,n]],x] /. x->Sqrt[-2] x) * Sqrt[-2] //Factor hdiffD[m_,n_] := D[scaleh[m,n],{x,2}] - (D[scaleh[m,n],x]^2/(2 scaleh[m,n]) + 3 scaleh[m,n]^3/2 + 4 x scaleh[m,n]^2 + 2(x^2-(m+2n+1)) scaleh[m,n] -2m^2/scaleh[m,n]) //Factor scaleg[m_,n_] := (D[Log[q[m,n+1]/q[m,n]],x] + w x /. x->Sqrt[-2/3] x) * Sqrt[-2/3] //Factor qdiffD[m_,n_] := D[scaleg[m,n],{x,2}] - (D[scaleg[m,n],x]^2/(2 scaleg[m,n]) + 3 scaleg[m,n]^3/2 + 4 x scaleg[m,n]^2 + 2(x^2-(m+2n)) scaleg[m,n] - 2/9 (3m-1)^2 /scaleg[m,n]) //Factor (* Next we check that there are no typos in the quadratic relations appearing in the paper. *) pquad1D[m_] := p[m+1] p[m-1] - (D[p[m],x]^2 - D[p[m],{x,2}] p[m] + x p[m]^2) //Expand pquad2D[m_] := p[m+2] p[m-2] - ((4 - (2 m + 1)^2) D[p[m],x]^2 + (x^4 + (2m+1)^2 x) p[m]^2 - 4 x^2 p[m] D[p[m],x] + (2m+1)^2 p[m] D[p[m],{x,2}]) //Expand pquad3D[m_] := p[m+2] p[m-1] - (-(2m+3) D[p[m],x] p[m+1] + (2m+1) D[p[m+1],x] p[m] + x^2 p[m] p[m+1]) //Expand hquad1D[m_,n_] := n h[m,n+1] h[m+1,n-1] - (h[m,n] D[h[m+1,n],x] - h[m+1,n] D[h[m,n],x]) //Expand hquad2D[m_,n_] := m h[m+1,n] h[m-1,n+1] - (h[m,n] D[h[m,n+1],x] - h[m,n+1] D[h[m,n],x]) //Expand hquad3D[m_,n_] := D[h[m,n],x]^2 - h[m,n] D[h[m,n],{x,2}] - (-m(h[m+1,n] h[m-1,n] - h[m,n]^2)) //Expand hquad4D[m_,n_] := D[h[m,n],x]^2 - h[m,n] D[h[m,n],{x,2}] - (n(h[m,n+1] h[m,n-1] - h[m,n]^2)) //Expand qquad1D[m_,n_] := q[m+1,n] q[m-1,n] - ((x^2 + 2 m + n - 1) q[m,n]^2 + D[q[m,n],{x,2}] q[m,n] - D[q[m,n],x]^2) //Factor qquad2D[m_,n_] := D[q[m,n],x] q[m+1,n] - (q[m,n] D[q[m+1,n],x] - x q[m,n] q[m+1,n] + q[m+1,n-1] q[m,n+1]) //Factor qquad3D[m_,n_] := (3(1-m-n)-1) q[m+1,n-1] q[m-1,n+1] + (3m-1) q[m,n+1] q[m,n-1] + (3n-1) q[m-1,n] q[m+1,n] //Factor (* Finally, we check that there are no typos in the discriminant and resultant formulas in the main theorems*) pdiscQ[m_] := disc[p[m]]/Product[j^(j(2m+1-j)^2/4),{j,3,2m-1,2}] pres1Q[m_] := Resultant[p[m],p[m-1],x]/ Product[j^(j(2m+1-j)(2m-1-j)/4),{j,3,2m-3,2}] pres2Q[m_] := Resultant[p[m],p[m-2],x]/ ((2m-1)^((m-1)(m-2)/2) Product[j^(j(2m+1-j)(2m-3-j)/4),{j,3,2m-5,2}]) star[-1] = star[0] = starstar[-1] = starstar[0] = 1 star[n_] := star[n] = Product[j^j,{j,1,n}] starstar[n_] := starstar[n] = Product[star[m],{m,1,n}] e[m_,n_,j_] := Which[ j <= Min[m,n], j^2 - 2(m-j)(n-j), Min[m,n] <= j <= Max[m,n],Min[m,n]^2, Max[m,n] <= j,(m+n-j)^2] hdisc[m_,n_] := Product[j^(j e[m,n,j]),{j,2,m+n-1}] hdiscQ[m_,n_] := disc[h[m,n]]/hdisc[m,n] hres1Q[m_,n_] := res[h[m,n],h[m+1,n]]/ (hdisc[m,n] starstar[m+n-1]/(starstar[m-1] star[m]^n starstar[n-1])) hres2Q[m_,n_] := res[h[m,n],h[m,n+1]]/ (hdisc[m,n] starstar[m+n-1]/(starstar[m-1] star[n]^m starstar[n-1])) hres3Q[m_,n_] := res[h[m,n],h[m+1,n-1]]/ (hdisc[m,n] star[n-1]^m/star[m]^n) f[x_,y_] := Which[ x >= 2 && y >= 2,Product[j^(j (x-(j+1)/3)(y - (j+1)/3)),{j,2,3 Min[x,y]-4,3}], x <=-2 && y <= -2, Product[j^(j(-x-(j-1)/3)(-y-(j-1)/3)),{j,4,3 Min[-x,-y]-2,3}], True,1] qdiscQ[m_,n_] := disc[q[m,n]]/(f[m,m] f[n,n] f[1-m-n,1-m-n]) qresQ[m_,n_] := res[q[m,n],q[m+1,n-1]]/(f[m,m+1] f[n,n-1] f[1-m-n,1-m-n])