Illustration of Zeros of Polynomials Orthogonal on the Unit Circle with Respect to Weights Having Jump Discontinuities in Higher Derivatives

Some Basic Definitions

The following command defines ell equally spaced angles of discontinuity in the interval (-π,π).  In order of increasing angle the angles are indexed by j varying between 1 and ell.  The set of angles is symmetric about zero.

thetajump[ell_Integer, j_Integer] = 
        Pi*(2*j - 1)/ell - Pi;

The following command defines the weight function ϕ(θ) as a function of the angle theta in the interval ending with thetajump[ell,j].  The vector w of length ell contains the parameters for the weight function, and the integer k determines the order of the derivative in which the weight will have jump discontinuities at the angles thetajump[ell,j].

phi[ell_Integer, w_, k_Integer, j_Integer, theta_] := 
        1 + Exp[w[[j]]]*
            Abs[Sin[ell*(theta - thetajump[ell, j])/2]]^k;

The following command defines the Fourier coefficient of the weight ϕ(θ) with wavenumber q.  Here q is allowed to be real so that a limit can be taken to integer values to correctly handle resonances.

FourierCoeff[ell_Integer, w_, k_Integer, q_]  :=  
        Sum[Integrate[
              phi[ell, w, k, j, theta]*Exp[-I*q*theta], {theta, 
                thetajump[ell, j] - 2*Pi/ell, thetajump[ell, j]}], {j, 1, 
              ell}]/(2*Pi);

The following command defines the Toeplitz matrix of dimension n+1 corresponding to the weight function ϕ(θ) specified by the number of discontinuities ell, the weight parameter vector w, and the order of the discontinuity k.

ToeplitzMatrix[ell_Integer, w_, k_Integer, n_Integer] := 
        Module[{FC = FourierCoeff[ell, w, k, q], fcs}, 
          fcs = Table[Limit[FC, q -> m], {m, -n, n}]; 
          Table[fcs[[j1 - j2 + n + 1]], {j1, 0, n}, {j2, 0, 
              n}]];

The following command defines the field V(θ)=-log(ϕ(θ)).

v[ell_Integer, w_, k_Integer, j_Integer, 
          theta_] := -Log[phi[ell, w, k, j, theta]];

The following command defines the mean value of V(θ) to be the result of a numerical integration to be carried out when the procedure is called.

v0[ell_Integer, w_, k_Integer] := 
        Sum[NIntegrate[
              v[ell, w, k, j, theta], {theta, thetajump[ell, j] - 2*Pi/ell, 
                thetajump[ell, j]}, WorkingPrecision -> 50], {j, 1, 
              ell}]/(2*Pi);

The following command returns a vector of the limiting values of the kth order derivative of V(θ) evaluated as the points of discontinuity are approached from below.

vderiv[ell_Integer, w_, k_Integer] := 
        Table[Limit[D[v[ell, w, k, j, theta], {theta, k}], 
            theta -> thetajump[ell, j]], {j, 1, ell}];

The following command uses the previous command to compute the jumps in the kth derivative of V(θ) at the points of discontinuity.

delta[ell_Integer, w_, 
          k_Integer] := ((-1)^k)*RotateLeft[vderiv[ell, w, k]] - 
          vderiv[ell, w, k];

The following command defines a vector of values of Ω(θ) (the Cauchy transform of V(θ)) at the points thetajump[ell,j] of discontinuity of the kth derivative of V(θ).  The values are obtained by numerical integration when the procedure is called.

omega[ell_Integer, w_, k_Integer] :=  
        Table[2*Im[-Sum[
                    NIntegrate[
                      I*Exp[I*
                            theta]*(v[ell, w, k, j2, theta] - 
                              v[ell, w, k, j2, thetajump[ell, j1]])/(Exp[
                                I*theta] - 
                              Exp[I*thetajump[ell, j1]]), {theta, 
                        thetajump[ell, j2] - 2*Pi/ell, 
                        thetajump[ell, j2]}], {j2, 1, ell}]/(2*Pi*
                    I)], {j1, 1, ell}];

The following command defines a polynomial proportional to the orthogonal polynomial of degree n with respect to the weight ϕ(θ) specified by the data ell (the number of points of discontinuity), w (the vector of weight parameters), and k (the order of the derivative having discontinuities).  The returned polynomial is a function of the variable z.  The coefficients are minors of a Toeplitz matrix, which are computed by first scaling the Toeplitz matrix elements by Exp[v0[ell,w,k]], which stabilizes the calculation with respect to numerical overflow or underflow according to the strong Szego limit theorem for Toeplitz determinants.  The calculation of the matrix elements is done to 50 digits of accuracy.

scaledpi[ell_Integer, w_, k_Integer, n_Integer, z_] := 
        Module[{t = ToeplitzMatrix[ell, w, k, n]}, 
          Sum[((-1)^(n - p))*(z^p)*
              Det[N[Exp[v0[ell, w, k]]*
                    t[[Table[m, {m, 1, n}], 
                        Join[Table[m, {m, 1, p}], 
                          Table[m, {m, p + 2, n + 1}]]]], 50]], {p, 0, 
              n}]];

The following command defines a constant multiple of the rational function whose zeros inside the unit circle are possible points of attraction for spurious zeros of the orthogonal polynomial of degree n.

scaledf[ell_Integer, w_, k_Integer, n_Integer, z_] := 
        Module[{deltavec = delta[ell, w, k], omegavec = omega[ell, w, k]}, 
          Sum[
deltavec[[j]]*
              Exp[I*(n + 1)*thetajump[ell, j]]*
              Exp[I*omegavec[[j]]]/(Exp[I*thetajump[ell, j]] - 
                    z), {j, 1, ell}]];

The Plot Generating Program

This code fragment takes as input the number ell of discontinuities of the derivative of order k of the weight, which is parameterized by a real vector w of length ell, and also the degree n of the polynomial of interest.

This module computes the orthogonal polynomial of degree n with respect to the weight ϕ(θ) on the unit circle specified by the remaining parameters.  The roots of this polynomial are then found numerically, and a graphics object called pidots is created to store them drawn as small black dots for subsequent display.  A graphics object called unitcircle is created storing the unit circle (drawn in black) with exterior tick marks indicating the angles thetajump[ell,j] of discontinuity of the kth derivative of the weight.  A graphics object called limitcircle contains the circle (drawn in green) to which the majority of the zeros are supposed to be attracted as the degree grows.  A graphics object called zerofree contains the circle (drawn in red) outside of which there should be no zeros of the polynomial (asymptotically as the degree grows).  The roots of the rational function whose roots are supposed to characterize the spurious zeros are computed and stored as lavender dots in a graphics object called fdots.  Finally, the graphics objects are all rendered together to give a complete picture of the zeros of the polynomial and their relation to various asymptotic quantities of interest.  This image is the output of the routine.

zeroplot[ell_, w_, k_, n_] := 
        Module[{polyroots = 
              z /. NSolve[scaledpi[ell, w, k, n, z] == 0, z], 
            polyrootsre, polyrootsim, 
            spuriousroots = 
              z /. NSolve[scaledf[ell, w, k, n, z] == 0, z], 
            spuriousrootsre, spuriousrootsim, pidots, unitcircle, tickmarks, 
            limitcircle, fdots}, 
polyrootsre = 
            Re[polyroots]; 
polyrootsim = 
            Im[polyroots]; 
spuriousrootsre = 
            Re[spuriousroots]; 
spuriousrootsim = 
            Im[spuriousroots]; 
pidots = 
            ListPlot[Transpose[{polyrootsre, polyrootsim}], 
              AspectRatio -> 1, DisplayFunction -> Identity, 
              PlotStyle -> 
                PointSize[1/68]]; 
unitcircle = 
            ParametricPlot[{Cos[theta], Sin[theta]}, {theta, 0, 2*Pi}, 
              DisplayFunction -> 
                Identity]; 
tickmarks = 
            Table[ParametricPlot[{r*Cos[thetajump[ell, j]], 
                  r*Sin[thetajump[ell, j]]}, {r, 1, 1.05}, 
                DisplayFunction -> Identity], {j, 1, 
                ell}]; 
limitcircle  =  
            ParametricPlot[{(1 - (k + 1)*Log[n]/n)*
                  Cos[theta], (1 - (k + 1)*Log[n]/n)*
                  Sin[theta]}, {theta, 0, 2*Pi}, PlotStyle -> Hue[1/3], 
              DisplayFunction -> 
                Identity]; 
zerofree = 
            ParametricPlot[{(1 - k*Log[n]/n)*
                  Cos[theta], (1 - k*Log[n]/n)*Sin[theta]}, {theta, 0, 
                2*Pi}, PlotStyle -> Hue[0], 
              DisplayFunction -> Identity]; 
fdots = 
            ListPlot[Transpose[{spuriousrootsre, spuriousrootsim}], 
              PlotStyle -> {Hue[2/3, 1/3, 1], PointSize[0.04]}, 
              DisplayFunction -> Identity];  
          Show[unitcircle, tickmarks, limitcircle, zerofree, fdots, pidots, 
            AspectRatio -> 1, 
            PlotRange -> {{-1.1, 1.1}, {-1.1, 1.1}}, 
            Axes -> False, Frame -> False, 
            DisplayFunction -> $DisplayFunction]];

Here is an example of its use.  Note that the previous definitions must be executed before zeroplot[] is called.

zeroplot[3, {1/2, -1/2, 0}, 2, 51]

A Version Without The Zeros of the Polynomial

This is the same as zeroplot[], but it does not calculate the orthogonal polynomial or its zeros.  This routine is useful for locating the possible spurious zeros, as it is significantly faster than zeroplot[] when n is large.

spuriousplot[ell_Integer, w_, k_Integer, 
          n_Integer] := 
Module[{spuriousroots = 
              z /. NSolve[scaledf[ell, w, k, n, z]], spuriousrootsre, 
            spuriousrootsim, unitcircle, tickmarks, limitcircle, fdots},  
          spuriousrootsre = Re[spuriousroots]; 
          spuriousrootsim = Im[spuriousroots];  
          unitcircle = 
            ParametricPlot[{Cos[theta], Sin[theta]}, {theta, 0, 2*Pi}, 
              DisplayFunction -> 
                Identity]; 
tickmarks = 
            Table[ParametricPlot[{r*Cos[thetajump[ell, j]], 
                  r*Sin[thetajump[ell, j]]}, {r, 1, 1.05}, 
                DisplayFunction -> Identity], {j, 1, 
                ell}]; 
limitcircle  =  
            ParametricPlot[{(1 - (k + 1)*Log[n]/n)*
                  Cos[theta], (1 - (k + 1)*Log[n]/n)*
                  Sin[theta]}, {theta, 0, 2*Pi}, PlotStyle -> Hue[1/3], 
              DisplayFunction -> 
                Identity]; 
zerofree = 
            ParametricPlot[{(1 - k*Log[n]/n)*
                  Cos[theta], (1 - k*Log[n]/n)*Sin[theta]}, {theta, 0, 
                2*Pi}, PlotStyle -> Hue[0], 
              DisplayFunction -> Identity]; 
fdots = 
            ListPlot[Transpose[{spuriousrootsre, spuriousrootsim}], 
              PlotStyle -> {Hue[2/3, 1/3, 1], PointSize[0.04]}, 
              DisplayFunction -> Identity];  
          Show[unitcircle, tickmarks, limitcircle, zerofree, fdots, 
            AspectRatio -> 1, 
            PlotRange -> {{-1.1, 1.1}, {-1.1, 1.1}}, 
            Axes -> False, Frame -> False, 
            DisplayFunction -> $DisplayFunction]];

Here is an example of its use.  The earlier definitions, with the exception of FourierCoeff[], ToeplitzMatrix[], and scaledpi[], must be executed before spuriousplot[] is called.

spuriousplot[3, {1/2, -1/2, 0}, 2, 51]