Problem3

from sage.plot.histogram import Histogram set_random_seed(635004) U = RealDistribution('uniform',[0,1]) T = RealDistribution('uniform',[0,2*pi]) W = RealDistribution('uniform',[1,2]) G = RealDistribution('gaussian', 1) C = ComplexField(1000) P.<z> = C['z'] 

# Compute and return array of roots def getRoot(poly): rootarray = [] for root in poly.roots(): multiplicity=root[1] for mult in range (multiplicity): rootarray.append(root[0]) return rootarray # Compute and return array of critical points def getCP(poly): derpoly = poly.derivative(z) cparray = [] for cp in derpoly.roots(): j = cp[1] for i in range(j): cparray.append(cp[0]) return cparray # Extract the real and imaginary parts of given array def getRealImag(array): real = [entry.real() for entry in array] imag = [entry.imag() for entry in array] return real, imag 

results1 = [] roots1 = [] cps1 = [] kabluchko1 = [] root_distribution1 = [] cp_distribution1 = [] # Take our f function to be an indicator function that returns 1 if the root's/critical point's imaginary part lies between 0 and 0.5 def f1(input): if 0 <= input.real() <= 0.5: return 1 else: return 0 # Because f ranges over 1/4 of a unit square expectation1 = 0.5 poly = C(1) for n in range(1, 101): poly *= (z - C(U.get_random_element()) - C(U.get_random_element())*C(I)) roots1 = getRoot(poly) cps1 = getCP(poly) # Compute sum of f(x_i) and sum of f(omega_i) roots_sum = 0 for root in roots1: roots_sum += f1(root) cps_sum = 0 for cp in cps1: cps_sum += f1(cp) # Show Kabluchko's theorem kabluchko = 1/n*roots_sum - 1/n*cps_sum kabluchko1.append(kabluchko) # Show our final result root_distribution = 1/sqrt(n)*(roots_sum-n*expectation1) root_distribution1.append(root_distribution) cp_distribution = 1/sqrt(n)*(cps_sum-(n-1)*expectation1) cp_distribution1.append(cp_distribution) result = root_distribution - cp_distribution results1.append(result) 

# Graph roots and critical points to make sure they look okay Gout1 = Graphics() Gout1 += points(roots1,color='white', markeredgecolor = 'red', pointsize=20,zorder = 1, marker = "o") Gout1 += points(cps1,color='white', markeredgecolor = 'blue', pointsize=20,zorder = 1, marker = "o") show(Gout1,aspect_ratio=1) 

# Graph Kabluchko's result Gout1 = list_plot([(i + 1, kabluchko1[i]) for i in range(len(kabluchko1))], color="green") show(Gout1) 

# Graph the difference for each order of polynomial Gout1 = list_plot([(i + 1, results1[i]) for i in range(len(results1))], color="green") show(Gout1) 

real_root_distribution1, imag_root_distribution1 = getRealImag(root_distribution1) real_cp_distribution1, imag_cp_distribution1 = getRealImag(cp_distribution1) 

# Plot a histogram to see distribution of critical points; may be trivial if n is small # We want to see a histogram that resembles a bell curve (Gaussian distribution!) histogram(imag_root_distribution1) 

histogram(real_cp_distribution1) 

results2 = [] roots2 = [] cps2 = [] kabluchko2 = [] # Take f to be an indicator function over the top half of the unit circle, i.e, when the imaginary part of an input is non-negative def f2(input): if input.imag() >= 0: return 1 else: return 0 expectation2 = 0.5 poly = C(1) for n in range(1, 201): t = T.get_random_element() poly *= (z - C(cos(t)) - C(sin(t))*C(I)) roots2 = getRoot(poly) cps2 = getCP(poly) # Compute sum of f(x_i) and sum of f(omega_i) roots_sum = 0 for root in roots2: roots_sum += f2(root) cps_sum = 0 for cp in cps2: cps_sum += f2(cp) # Show Kabluchko's theorem kabluchko = 1/n*roots_sum - 1/n*cps_sum kabluchko2.append(kabluchko) # Show our final result result = 1/sqrt(n)*(roots_sum-n*expectation2) - 1/sqrt(n)*(cps_sum-(n-1)*expectation2) results2.append(result) 

Gout2 = Graphics() Gout2 += points(roots2,color='white', markeredgecolor = 'red', pointsize=20,zorder = 1, marker = "o") Gout2 += points(cps2,color='white', markeredgecolor = 'blue', pointsize=20,zorder = 1, marker = "o") show(Gout2,aspect_ratio=1) 

# Graph Kabluchko's result Gout2 = list_plot([(i + 1, kabluchko2[i]) for i in range(len(kabluchko2))], color="green") show(Gout2) 

Gout2 = list_plot([(i + 1, results2[i]) for i in range(len(results2))], color="green") show(Gout2) 

results3 = [] roots3 = [] cps3 = [] kabluchko3 = [] root_distribution3 = [] cp_distribution3 = [] # Take our f function to be an indicator function that returns 1 if the root's/critical point's imaginary part is between 1 and 1.2 def f3(input): if 1 <= input.imag() <= 1.2: return 1 else: return 0 expectation3 = 0.2 poly = C(1) for n in range(1, 201): poly *= (z - C(W.get_random_element()) - C(W.get_random_element())*C(I)) roots3 = getRoot(poly) cps3 = getCP(poly) # Compute sum of f(x_i) and sum of f(omega_i) roots_sum = 0 for root in roots3: roots_sum += f3(root) cps_sum = 0 for cp in cps3: cps_sum += f3(cp) # Show Kabluchko's theorem kabluchko = 1/n*roots_sum - 1/n*cps_sum kabluchko3.append(kabluchko) # Show our final result root_distribution = 1/sqrt(n)*(roots_sum-n*expectation3) root_distribution3.append(root_distribution) cp_distribution = 1/sqrt(n)*(cps_sum-(n-1)*expectation3) cp_distribution3.append(cp_distribution) result = root_distribution - cp_distribution results3.append(result) 

Gout3 = Graphics() Gout3 += points(roots3,color='white', markeredgecolor = 'red', pointsize=20,zorder = 1, marker = "o") Gout3 += points(cps3,color='white', markeredgecolor = 'blue', pointsize=20,zorder = 1, marker = "o") show(Gout3,aspect_ratio=1) 

Gout3 = list_plot([(i + 1, kabluchko3[i]) for i in range(len(kabluchko3))], color="green") show(Gout3) 

Gout3 = list_plot([(i + 1, results3[i]) for i in range(len(results3))], color="green") show(Gout3) 

Gout3 = list_plot([(i + 1, kabluchko3[i]) for i in range(len(kabluchko3))], color='red', pointsize=15,zorder = 1, marker = "o", alpha=.7, legend_label="Kabluchko") Gout3 += list_plot([(i + 1, results3[i]) for i in range(len(results3))], color='blue', pointsize=15,zorder = 1, marker = "o", alpha=.7, legend_label="Problem 3", axes_labels=["n",""], fontsize=8) show(Gout3) 

real_root_distribution3, imag_root_distribution3 = getRealImag(root_distribution3) real_cp_distribution3, imag_cp_distribution3 = getRealImag(cp_distribution3) 

histogram(real_cp_distribution3) 

histogram(imag_cp_distribution3) 

results4 = [] roots4 = [] cps4 = [] kabluchko4 = [] root_distribution4 = [] cp_distribution4 = [] # Take our f function to be the identity function def f4(input): return input expectation4 = 0 poly = C(1) for n in range(1, 201): poly *= (z - C(G.get_random_element()) - C(G.get_random_element())*C(I)) roots4 = getRoot(poly) cps4 = getCP(poly) # Compute sum of f(x_i) and sum of f(omega_i) roots_sum = 0 for root in roots4: roots_sum += f4(root) cps_sum = 0 for cp in cps4: cps_sum += f4(cp) # Show Kabluchko's theorem kabluchko = 1/n*roots_sum - 1/n*cps_sum kabluchko4.append(kabluchko) # Show our final result root_distribution = 1/sqrt(n)*(roots_sum-n*expectation4) root_distribution4.append(root_distribution) cp_distribution = 1/sqrt(n)*(cps_sum-(n-1)*expectation4) cp_distribution4.append(cp_distribution) result = root_distribution - cp_distribution results4.append(result) 

Gout4 = Graphics() Gout4 += points(roots4,color='white', markeredgecolor = 'red', pointsize=20,zorder = 1, marker = "o", legend_label="Roots") Gout4 += points(cps4,color='white', markeredgecolor = 'blue', pointsize=20,zorder = 1, marker = "o", axes_labels=["Re","Im"], fontsize=8, legend_label="Critical points") show(Gout4,aspect_ratio=1) 

Gout4 = list_plot([(i + 1, kabluchko4[i]) for i in range(len(kabluchko4))], color="green") show(Gout4) 

Gout4 = list_plot([(i + 1, results4[i]) for i in range(len(results4))], color="green") show(Gout4) 

Gout4 = list_plot([(i + 1, kabluchko4[i]) for i in range(len(kabluchko4))], color='white', markeredgecolor = 'red', pointsize=15,zorder = 1, marker = "o", alpha=.7, legend_label="Kabluchko") Gout4 += list_plot([(i + 1, results4[i]) for i in range(len(results4))], color='white', markeredgecolor = 'blue', pointsize=15,zorder = 1, marker = "o", alpha=.7, legend_label="Problem 3", axes_labels=["Re","Im"], fontsize=8) show(Gout4) 

real_root_distribution4, imag_root_distribution4 = getRealImag(root_distribution4) real_cp_distribution4, imag_cp_distribution4 = getRealImag(cp_distribution4) 

histogram(real_cp_distribution4) 

histogram(imag_cp_distribution4)