I’m trying to approximate Julia sets using roots of polynomials in Python. In particular I want to find the roots of the nth iterate of the polynomial q(z) = z^2-0.5. In other words I want to find the roots of q(q(q(..))) composed n times. In order to get many sample points I would need to compute the roots of polynomials of degree >1000.
I’ve tried solving this problem using both the built in polynomial class of numpy which has a root function and also the function solver of sympy. In the first case precision is lost when I choose degrees larger than 100. The sympy computation simply takes to long time. Here is my code:
p = P([-0.5,0,1])
for k in range(9):
p = p**2-0.5
roots = p.roots()
plt.plot([np.real(r) for r in roots], [np.imag(r) for r in roots],'x')
plt.show()
abs_vector = [np.abs(p(r)) for r in roots]
max = 0
for a in abs_vector:
if a > max:
max = a
print(max)
The max value above gives the largest value of p at a supposed root. However running this code gives me 7.881370400084486e+296 which is very large.
How would one go about computing roots of high degree polynomials with good accuracy in a short amount of time?
Advertisement
Answer
For the n-times composition of a polynomial q you can reconstruct the roots iteratively
q = [1,0,-0.5]
n = 9
def q_preimage(w):
c = q.copy()
c[-1] -= w
return np.roots(c)
rts = [0]
for k in range(n):
rts = np.concatenate([q_preimage(w) for w in rts])
which returns
array([ 1.36444432e+00+0.00095319j, -1.36444432e+00-0.00095319j,
1.40104860e-03-0.92828301j, -1.40104860e-03+0.92828301j,
8.82183775e-01-0.52384727j, -8.82183775e-01+0.52384727j,
8.78972436e-01+0.52576116j, -8.78972436e-01-0.52576116j,
1.19545693e+00-0.21647154j, -1.19545693e+00+0.21647154j,
3.61362916e-01+0.71612883j, -3.61362916e-01-0.71612883j,
1.19225541e+00+0.21925381j, -1.19225541e+00-0.21925381j,
3.66786415e-01-0.71269419j, -3.66786415e-01+0.71269419j,
...
or plotted
plt.plot(rts.real, rts.imag,'ob', ms=2); plt.grid(); plt.show()
