ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Wed, 20 May 2020 11:35:19 +0200Problem with implicit_plothttps://ask.sagemath.org/question/51443/problem-with-implicit_plot/I have three huge degree 31 bivariate polynomials (20,000 characters long each) I want to plot, but I keep getting a lot of noise in my plot. I can't upload it, but the point is that in some regions I just get colorful noise. I've tried defining the polynomials over RealField(n) and increasing the number of plot_points, but neither of these approaches work. Any ideas on how to work around this? Thanks.
**Edit:** Tried using sympy's plot_implicit and it's so (SO!) slow. Then used numpy's contour_plot and it's fast, but has the same problem as sage.
Here's the code that produces the polynomials and plot. Be patient as it could be a bit slow (depending on your machine).
plot = Graphics()
m = 32-1
D = [(i,j) for i in range(0,m+1) for j in range(0,m+1) if i+j<m+1]
#Polygon
P = Polyhedron( vertices= [(0, 0), (32, 44), (23, 0), (10, 14), (2, 3)] )
points = P.integral_points()
plot_pts = point(points, rgbcolor=(0, 0, 0), size = 20).plot()
plot_np = P.plot(fill = False, point=False, line='black')
M = matrix(ZZ, len(points), len(D), 0)
for row_num, row in enumerate(points):
for col_num, column in enumerate(D):
i, j = row
a, b= column
#Matrix for interpolation:
M[row_num, col_num] = (i^a)*(j^b)
R = PolynomialRing(QQ, 2, 'xy')
S = PolynomialRing(RealField(500), 2, 'uv')
x, y = R.gens()
u, v = S.gens()
K = M.right_kernel()
Kdim = K.dimension()
print(Kdim)
if Kdim > 0:
for l in range(Kdim):
K_basis = K.basis()[l]
#Writing the interpolating polynomial
f=0
for order, bidegree in enumerate(D):
a, b = bidegree
f += list(K_basis)[order] * u^a * v^b
F = f.factor()
f = F[6][0]
cols = ['red', 'blue', 'green']
interpolation = implicit_plot(f, (u,-1,34), (v,-1,45), plot_points=100, color=cols[l])
plot += interpolation
plot += plot_pts + plot_np
plot.show(figsize=10)
**Edit 2:** Using the mpmath library in Python and with the aid of Sébastien's code below I wrote a routine that allows us to control the root finding method and precision of our computations. I tried several methods, secant (default), newton, hailley, mnewton, etc. without success. Changing the precision and tolerance of the root finding function from mpmath didn't help either. I think this polynomial just behaves too wildly in the region of the plot.
Here's the code and relevant documentation for the "mpmath.findroot" function:
http://mpmath.org/doc/current/calculus/optimization.html :
from sympy import *
import matplotlib.pyplot as plt
import mpmath as mp
mp.dps = 100
stepx = 0.1
stepy = 0.5
xrange = np.arange(7.5,12.5, stepx)
yrange = np.arange(0, 5, stepy)
def plot_roots_of_f(f):
L = []
for u in xrange:
for v in yrange:
Root = mp.findroot(f.eval({x:u}),v, method = 'mnewton', tol = E-60, verbose=False,verify=False)
L.append([u,Root])
return plt.plot(*zip(*L),linestyle='None', marker='.')
plot_roots_of_f(f)
plt.show()Mon, 18 May 2020 09:30:23 +0200https://ask.sagemath.org/question/51443/problem-with-implicit_plot/Comment by tmonteil for <p>I have three huge degree 31 bivariate polynomials (20,000 characters long each) I want to plot, but I keep getting a lot of noise in my plot. I can't upload it, but the point is that in some regions I just get colorful noise. I've tried defining the polynomials over RealField(n) and increasing the number of plot_points, but neither of these approaches work. Any ideas on how to work around this? Thanks.</p>
<p><strong>Edit:</strong> Tried using sympy's plot_implicit and it's so (SO!) slow. Then used numpy's contour_plot and it's fast, but has the same problem as sage.</p>
<p>Here's the code that produces the polynomials and plot. Be patient as it could be a bit slow (depending on your machine).</p>
<pre><code>plot = Graphics()
m = 32-1
D = [(i,j) for i in range(0,m+1) for j in range(0,m+1) if i+j<m+1]
#Polygon
P = Polyhedron( vertices= [(0, 0), (32, 44), (23, 0), (10, 14), (2, 3)] )
points = P.integral_points()
plot_pts = point(points, rgbcolor=(0, 0, 0), size = 20).plot()
plot_np = P.plot(fill = False, point=False, line='black')
M = matrix(ZZ, len(points), len(D), 0)
for row_num, row in enumerate(points):
for col_num, column in enumerate(D):
i, j = row
a, b= column
#Matrix for interpolation:
M[row_num, col_num] = (i^a)*(j^b)
R = PolynomialRing(QQ, 2, 'xy')
S = PolynomialRing(RealField(500), 2, 'uv')
x, y = R.gens()
u, v = S.gens()
K = M.right_kernel()
Kdim = K.dimension()
print(Kdim)
if Kdim > 0:
for l in range(Kdim):
K_basis = K.basis()[l]
#Writing the interpolating polynomial
f=0
for order, bidegree in enumerate(D):
a, b = bidegree
f += list(K_basis)[order] * u^a * v^b
F = f.factor()
f = F[6][0]
cols = ['red', 'blue', 'green']
interpolation = implicit_plot(f, (u,-1,34), (v,-1,45), plot_points=100, color=cols[l])
plot += interpolation
plot += plot_pts + plot_np
plot.show(figsize=10)
</code></pre>
<p><strong>Edit 2:</strong> Using the mpmath library in Python and with the aid of Sébastien's code below I wrote a routine that allows us to control the root finding method and precision of our computations. I tried several methods, secant (default), newton, hailley, mnewton, etc. without success. Changing the precision and tolerance of the root finding function from mpmath didn't help either. I think this polynomial just behaves too wildly in the region of the plot.</p>
<p>Here's the code and relevant documentation for the "mpmath.findroot" function: </p>
<p><a href="http://mpmath.org/doc/current/calculus/optimization.html">http://mpmath.org/doc/current/calculu...</a> :</p>
<pre><code>from sympy import *
import matplotlib.pyplot as plt
import mpmath as mp
mp.dps = 100
stepx = 0.1
stepy = 0.5
xrange = np.arange(7.5,12.5, stepx)
yrange = np.arange(0, 5, stepy)
def plot_roots_of_f(f):
L = []
for u in xrange:
for v in yrange:
Root = mp.findroot(f.eval({x:u}),v, method = 'mnewton', tol = E-60, verbose=False,verify=False)
L.append([u,Root])
return plt.plot(*zip(*L),linestyle='None', marker='.')
plot_roots_of_f(f)
plt.show()
</code></pre>
https://ask.sagemath.org/question/51443/problem-with-implicit_plot/?comment=51449#post-id-51449Could you please try to upload the code somewhere so taht we can give a try ?Mon, 18 May 2020 19:11:28 +0200https://ask.sagemath.org/question/51443/problem-with-implicit_plot/?comment=51449#post-id-51449Comment by jga for <p>I have three huge degree 31 bivariate polynomials (20,000 characters long each) I want to plot, but I keep getting a lot of noise in my plot. I can't upload it, but the point is that in some regions I just get colorful noise. I've tried defining the polynomials over RealField(n) and increasing the number of plot_points, but neither of these approaches work. Any ideas on how to work around this? Thanks.</p>
<p><strong>Edit:</strong> Tried using sympy's plot_implicit and it's so (SO!) slow. Then used numpy's contour_plot and it's fast, but has the same problem as sage.</p>
<p>Here's the code that produces the polynomials and plot. Be patient as it could be a bit slow (depending on your machine).</p>
<pre><code>plot = Graphics()
m = 32-1
D = [(i,j) for i in range(0,m+1) for j in range(0,m+1) if i+j<m+1]
#Polygon
P = Polyhedron( vertices= [(0, 0), (32, 44), (23, 0), (10, 14), (2, 3)] )
points = P.integral_points()
plot_pts = point(points, rgbcolor=(0, 0, 0), size = 20).plot()
plot_np = P.plot(fill = False, point=False, line='black')
M = matrix(ZZ, len(points), len(D), 0)
for row_num, row in enumerate(points):
for col_num, column in enumerate(D):
i, j = row
a, b= column
#Matrix for interpolation:
M[row_num, col_num] = (i^a)*(j^b)
R = PolynomialRing(QQ, 2, 'xy')
S = PolynomialRing(RealField(500), 2, 'uv')
x, y = R.gens()
u, v = S.gens()
K = M.right_kernel()
Kdim = K.dimension()
print(Kdim)
if Kdim > 0:
for l in range(Kdim):
K_basis = K.basis()[l]
#Writing the interpolating polynomial
f=0
for order, bidegree in enumerate(D):
a, b = bidegree
f += list(K_basis)[order] * u^a * v^b
F = f.factor()
f = F[6][0]
cols = ['red', 'blue', 'green']
interpolation = implicit_plot(f, (u,-1,34), (v,-1,45), plot_points=100, color=cols[l])
plot += interpolation
plot += plot_pts + plot_np
plot.show(figsize=10)
</code></pre>
<p><strong>Edit 2:</strong> Using the mpmath library in Python and with the aid of Sébastien's code below I wrote a routine that allows us to control the root finding method and precision of our computations. I tried several methods, secant (default), newton, hailley, mnewton, etc. without success. Changing the precision and tolerance of the root finding function from mpmath didn't help either. I think this polynomial just behaves too wildly in the region of the plot.</p>
<p>Here's the code and relevant documentation for the "mpmath.findroot" function: </p>
<p><a href="http://mpmath.org/doc/current/calculus/optimization.html">http://mpmath.org/doc/current/calculu...</a> :</p>
<pre><code>from sympy import *
import matplotlib.pyplot as plt
import mpmath as mp
mp.dps = 100
stepx = 0.1
stepy = 0.5
xrange = np.arange(7.5,12.5, stepx)
yrange = np.arange(0, 5, stepy)
def plot_roots_of_f(f):
L = []
for u in xrange:
for v in yrange:
Root = mp.findroot(f.eval({x:u}),v, method = 'mnewton', tol = E-60, verbose=False,verify=False)
L.append([u,Root])
return plt.plot(*zip(*L),linestyle='None', marker='.')
plot_roots_of_f(f)
plt.show()
</code></pre>
https://ask.sagemath.org/question/51443/problem-with-implicit_plot/?comment=51455#post-id-51455I added the code, @tmonteilTue, 19 May 2020 00:42:34 +0200https://ask.sagemath.org/question/51443/problem-with-implicit_plot/?comment=51455#post-id-51455Comment by Sébastien for <p>I have three huge degree 31 bivariate polynomials (20,000 characters long each) I want to plot, but I keep getting a lot of noise in my plot. I can't upload it, but the point is that in some regions I just get colorful noise. I've tried defining the polynomials over RealField(n) and increasing the number of plot_points, but neither of these approaches work. Any ideas on how to work around this? Thanks.</p>
<p><strong>Edit:</strong> Tried using sympy's plot_implicit and it's so (SO!) slow. Then used numpy's contour_plot and it's fast, but has the same problem as sage.</p>
<p>Here's the code that produces the polynomials and plot. Be patient as it could be a bit slow (depending on your machine).</p>
<pre><code>plot = Graphics()
m = 32-1
D = [(i,j) for i in range(0,m+1) for j in range(0,m+1) if i+j<m+1]
#Polygon
P = Polyhedron( vertices= [(0, 0), (32, 44), (23, 0), (10, 14), (2, 3)] )
points = P.integral_points()
plot_pts = point(points, rgbcolor=(0, 0, 0), size = 20).plot()
plot_np = P.plot(fill = False, point=False, line='black')
M = matrix(ZZ, len(points), len(D), 0)
for row_num, row in enumerate(points):
for col_num, column in enumerate(D):
i, j = row
a, b= column
#Matrix for interpolation:
M[row_num, col_num] = (i^a)*(j^b)
R = PolynomialRing(QQ, 2, 'xy')
S = PolynomialRing(RealField(500), 2, 'uv')
x, y = R.gens()
u, v = S.gens()
K = M.right_kernel()
Kdim = K.dimension()
print(Kdim)
if Kdim > 0:
for l in range(Kdim):
K_basis = K.basis()[l]
#Writing the interpolating polynomial
f=0
for order, bidegree in enumerate(D):
a, b = bidegree
f += list(K_basis)[order] * u^a * v^b
F = f.factor()
f = F[6][0]
cols = ['red', 'blue', 'green']
interpolation = implicit_plot(f, (u,-1,34), (v,-1,45), plot_points=100, color=cols[l])
plot += interpolation
plot += plot_pts + plot_np
plot.show(figsize=10)
</code></pre>
<p><strong>Edit 2:</strong> Using the mpmath library in Python and with the aid of Sébastien's code below I wrote a routine that allows us to control the root finding method and precision of our computations. I tried several methods, secant (default), newton, hailley, mnewton, etc. without success. Changing the precision and tolerance of the root finding function from mpmath didn't help either. I think this polynomial just behaves too wildly in the region of the plot.</p>
<p>Here's the code and relevant documentation for the "mpmath.findroot" function: </p>
<p><a href="http://mpmath.org/doc/current/calculus/optimization.html">http://mpmath.org/doc/current/calculu...</a> :</p>
<pre><code>from sympy import *
import matplotlib.pyplot as plt
import mpmath as mp
mp.dps = 100
stepx = 0.1
stepy = 0.5
xrange = np.arange(7.5,12.5, stepx)
yrange = np.arange(0, 5, stepy)
def plot_roots_of_f(f):
L = []
for u in xrange:
for v in yrange:
Root = mp.findroot(f.eval({x:u}),v, method = 'mnewton', tol = E-60, verbose=False,verify=False)
L.append([u,Root])
return plt.plot(*zip(*L),linestyle='None', marker='.')
plot_roots_of_f(f)
plt.show()
</code></pre>
https://ask.sagemath.org/question/51443/problem-with-implicit_plot/?comment=51456#post-id-51456The code runs in few seconds on my machine (does not seem computationally heavy?). The output is the overlay the three implicit plot which indeed looks like "*colorful noise*". But what is the expected result? Is the equation `f==0` suppose to pass through all of the integer points inside the polygon?Tue, 19 May 2020 09:41:34 +0200https://ask.sagemath.org/question/51443/problem-with-implicit_plot/?comment=51456#post-id-51456Comment by jga for <p>I have three huge degree 31 bivariate polynomials (20,000 characters long each) I want to plot, but I keep getting a lot of noise in my plot. I can't upload it, but the point is that in some regions I just get colorful noise. I've tried defining the polynomials over RealField(n) and increasing the number of plot_points, but neither of these approaches work. Any ideas on how to work around this? Thanks.</p>
<p><strong>Edit:</strong> Tried using sympy's plot_implicit and it's so (SO!) slow. Then used numpy's contour_plot and it's fast, but has the same problem as sage.</p>
<p>Here's the code that produces the polynomials and plot. Be patient as it could be a bit slow (depending on your machine).</p>
<pre><code>plot = Graphics()
m = 32-1
D = [(i,j) for i in range(0,m+1) for j in range(0,m+1) if i+j<m+1]
#Polygon
P = Polyhedron( vertices= [(0, 0), (32, 44), (23, 0), (10, 14), (2, 3)] )
points = P.integral_points()
plot_pts = point(points, rgbcolor=(0, 0, 0), size = 20).plot()
plot_np = P.plot(fill = False, point=False, line='black')
M = matrix(ZZ, len(points), len(D), 0)
for row_num, row in enumerate(points):
for col_num, column in enumerate(D):
i, j = row
a, b= column
#Matrix for interpolation:
M[row_num, col_num] = (i^a)*(j^b)
R = PolynomialRing(QQ, 2, 'xy')
S = PolynomialRing(RealField(500), 2, 'uv')
x, y = R.gens()
u, v = S.gens()
K = M.right_kernel()
Kdim = K.dimension()
print(Kdim)
if Kdim > 0:
for l in range(Kdim):
K_basis = K.basis()[l]
#Writing the interpolating polynomial
f=0
for order, bidegree in enumerate(D):
a, b = bidegree
f += list(K_basis)[order] * u^a * v^b
F = f.factor()
f = F[6][0]
cols = ['red', 'blue', 'green']
interpolation = implicit_plot(f, (u,-1,34), (v,-1,45), plot_points=100, color=cols[l])
plot += interpolation
plot += plot_pts + plot_np
plot.show(figsize=10)
</code></pre>
<p><strong>Edit 2:</strong> Using the mpmath library in Python and with the aid of Sébastien's code below I wrote a routine that allows us to control the root finding method and precision of our computations. I tried several methods, secant (default), newton, hailley, mnewton, etc. without success. Changing the precision and tolerance of the root finding function from mpmath didn't help either. I think this polynomial just behaves too wildly in the region of the plot.</p>
<p>Here's the code and relevant documentation for the "mpmath.findroot" function: </p>
<p><a href="http://mpmath.org/doc/current/calculus/optimization.html">http://mpmath.org/doc/current/calculu...</a> :</p>
<pre><code>from sympy import *
import matplotlib.pyplot as plt
import mpmath as mp
mp.dps = 100
stepx = 0.1
stepy = 0.5
xrange = np.arange(7.5,12.5, stepx)
yrange = np.arange(0, 5, stepy)
def plot_roots_of_f(f):
L = []
for u in xrange:
for v in yrange:
Root = mp.findroot(f.eval({x:u}),v, method = 'mnewton', tol = E-60, verbose=False,verify=False)
L.append([u,Root])
return plt.plot(*zip(*L),linestyle='None', marker='.')
plot_roots_of_f(f)
plt.show()
</code></pre>
https://ask.sagemath.org/question/51443/problem-with-implicit_plot/?comment=51457#post-id-51457@Sébastien exactly. These three polynomials should define nice graphs passing through said points.This must be the case because for each fixed value of x there can be at most 31 values of y in the plot. These vary continuously as x does, so there's just no way it looks as the plot suggests. In my case this is very clear for x<10, where the plot looks well. The issue happens further away from the origin.Tue, 19 May 2020 10:15:44 +0200https://ask.sagemath.org/question/51443/problem-with-implicit_plot/?comment=51457#post-id-51457Comment by Sébastien for <p>I have three huge degree 31 bivariate polynomials (20,000 characters long each) I want to plot, but I keep getting a lot of noise in my plot. I can't upload it, but the point is that in some regions I just get colorful noise. I've tried defining the polynomials over RealField(n) and increasing the number of plot_points, but neither of these approaches work. Any ideas on how to work around this? Thanks.</p>
<p><strong>Edit:</strong> Tried using sympy's plot_implicit and it's so (SO!) slow. Then used numpy's contour_plot and it's fast, but has the same problem as sage.</p>
<p>Here's the code that produces the polynomials and plot. Be patient as it could be a bit slow (depending on your machine).</p>
<pre><code>plot = Graphics()
m = 32-1
D = [(i,j) for i in range(0,m+1) for j in range(0,m+1) if i+j<m+1]
#Polygon
P = Polyhedron( vertices= [(0, 0), (32, 44), (23, 0), (10, 14), (2, 3)] )
points = P.integral_points()
plot_pts = point(points, rgbcolor=(0, 0, 0), size = 20).plot()
plot_np = P.plot(fill = False, point=False, line='black')
M = matrix(ZZ, len(points), len(D), 0)
for row_num, row in enumerate(points):
for col_num, column in enumerate(D):
i, j = row
a, b= column
#Matrix for interpolation:
M[row_num, col_num] = (i^a)*(j^b)
R = PolynomialRing(QQ, 2, 'xy')
S = PolynomialRing(RealField(500), 2, 'uv')
x, y = R.gens()
u, v = S.gens()
K = M.right_kernel()
Kdim = K.dimension()
print(Kdim)
if Kdim > 0:
for l in range(Kdim):
K_basis = K.basis()[l]
#Writing the interpolating polynomial
f=0
for order, bidegree in enumerate(D):
a, b = bidegree
f += list(K_basis)[order] * u^a * v^b
F = f.factor()
f = F[6][0]
cols = ['red', 'blue', 'green']
interpolation = implicit_plot(f, (u,-1,34), (v,-1,45), plot_points=100, color=cols[l])
plot += interpolation
plot += plot_pts + plot_np
plot.show(figsize=10)
</code></pre>
<p><strong>Edit 2:</strong> Using the mpmath library in Python and with the aid of Sébastien's code below I wrote a routine that allows us to control the root finding method and precision of our computations. I tried several methods, secant (default), newton, hailley, mnewton, etc. without success. Changing the precision and tolerance of the root finding function from mpmath didn't help either. I think this polynomial just behaves too wildly in the region of the plot.</p>
<p>Here's the code and relevant documentation for the "mpmath.findroot" function: </p>
<p><a href="http://mpmath.org/doc/current/calculus/optimization.html">http://mpmath.org/doc/current/calculu...</a> :</p>
<pre><code>from sympy import *
import matplotlib.pyplot as plt
import mpmath as mp
mp.dps = 100
stepx = 0.1
stepy = 0.5
xrange = np.arange(7.5,12.5, stepx)
yrange = np.arange(0, 5, stepy)
def plot_roots_of_f(f):
L = []
for u in xrange:
for v in yrange:
Root = mp.findroot(f.eval({x:u}),v, method = 'mnewton', tol = E-60, verbose=False,verify=False)
L.append([u,Root])
return plt.plot(*zip(*L),linestyle='None', marker='.')
plot_roots_of_f(f)
plt.show()
</code></pre>
https://ask.sagemath.org/question/51443/problem-with-implicit_plot/?comment=51459#post-id-51459I think the generic `implicit_plot` function is no good for your function `f` which involves sums of large numbers (possibly done with using float inputs) that must equal to zero at the end. I would suggest you to write your own drawing function based on factorization of the univariate polynomial `f` after evaluation at specific values of `u`. For example:
sage: f.subs(u=14.7).factor()
(-1.02680916303147e30) * (v - 10.0971622442134) * (v - 2.04438249413817) * (v + 1.88294590341942) * (v + 12.3103608056188) * (v + 29.6136637586803) * (v^2 - 56.4157620800683*v + 804.755073557794) * (v^2 - 50.4692924243472*v + 691.858754703194) * (v^2 - 44.2941922218622*v + 606.707600654288) * (v^2 - 33.9849183617846*v + 477.330795667438) * (v^2 - 21.5731314846998*v + 318.8705...) * (v^2 - 10.1....)Tue, 19 May 2020 13:53:18 +0200https://ask.sagemath.org/question/51443/problem-with-implicit_plot/?comment=51459#post-id-51459Answer by Sébastien for <p>I have three huge degree 31 bivariate polynomials (20,000 characters long each) I want to plot, but I keep getting a lot of noise in my plot. I can't upload it, but the point is that in some regions I just get colorful noise. I've tried defining the polynomials over RealField(n) and increasing the number of plot_points, but neither of these approaches work. Any ideas on how to work around this? Thanks.</p>
<p><strong>Edit:</strong> Tried using sympy's plot_implicit and it's so (SO!) slow. Then used numpy's contour_plot and it's fast, but has the same problem as sage.</p>
<p>Here's the code that produces the polynomials and plot. Be patient as it could be a bit slow (depending on your machine).</p>
<pre><code>plot = Graphics()
m = 32-1
D = [(i,j) for i in range(0,m+1) for j in range(0,m+1) if i+j<m+1]
#Polygon
P = Polyhedron( vertices= [(0, 0), (32, 44), (23, 0), (10, 14), (2, 3)] )
points = P.integral_points()
plot_pts = point(points, rgbcolor=(0, 0, 0), size = 20).plot()
plot_np = P.plot(fill = False, point=False, line='black')
M = matrix(ZZ, len(points), len(D), 0)
for row_num, row in enumerate(points):
for col_num, column in enumerate(D):
i, j = row
a, b= column
#Matrix for interpolation:
M[row_num, col_num] = (i^a)*(j^b)
R = PolynomialRing(QQ, 2, 'xy')
S = PolynomialRing(RealField(500), 2, 'uv')
x, y = R.gens()
u, v = S.gens()
K = M.right_kernel()
Kdim = K.dimension()
print(Kdim)
if Kdim > 0:
for l in range(Kdim):
K_basis = K.basis()[l]
#Writing the interpolating polynomial
f=0
for order, bidegree in enumerate(D):
a, b = bidegree
f += list(K_basis)[order] * u^a * v^b
F = f.factor()
f = F[6][0]
cols = ['red', 'blue', 'green']
interpolation = implicit_plot(f, (u,-1,34), (v,-1,45), plot_points=100, color=cols[l])
plot += interpolation
plot += plot_pts + plot_np
plot.show(figsize=10)
</code></pre>
<p><strong>Edit 2:</strong> Using the mpmath library in Python and with the aid of Sébastien's code below I wrote a routine that allows us to control the root finding method and precision of our computations. I tried several methods, secant (default), newton, hailley, mnewton, etc. without success. Changing the precision and tolerance of the root finding function from mpmath didn't help either. I think this polynomial just behaves too wildly in the region of the plot.</p>
<p>Here's the code and relevant documentation for the "mpmath.findroot" function: </p>
<p><a href="http://mpmath.org/doc/current/calculus/optimization.html">http://mpmath.org/doc/current/calculu...</a> :</p>
<pre><code>from sympy import *
import matplotlib.pyplot as plt
import mpmath as mp
mp.dps = 100
stepx = 0.1
stepy = 0.5
xrange = np.arange(7.5,12.5, stepx)
yrange = np.arange(0, 5, stepy)
def plot_roots_of_f(f):
L = []
for u in xrange:
for v in yrange:
Root = mp.findroot(f.eval({x:u}),v, method = 'mnewton', tol = E-60, verbose=False,verify=False)
L.append([u,Root])
return plt.plot(*zip(*L),linestyle='None', marker='.')
plot_roots_of_f(f)
plt.show()
</code></pre>
https://ask.sagemath.org/question/51443/problem-with-implicit_plot/?answer=51460#post-id-51460It seems that the polynomial f itself does not contain the good information.
sage: f.parent()
Multivariate Polynomial Ring in u, v over Real Field with 500 bits of precision
Because when I draw its roots with:
def roots_of_f(start,stop,step):
u_range = srange(start, stop, step)
L = [(u,v) for u in u_range
for v in f.subs(u=u).univariate_polynomial().roots(multiplicities=False)
if (u,v) in P]
return L
I get (`points` is a list in the code and overwrites the sage function `points`, so I am using `point` below):
sage: point(roots_of_f(1,20,.05)) + plot_np
![image description](/upfiles/15898901611271885.png)
Replacing `u` and `v` by `x` and `y` to use the ring over the rational field seems to give the same thing.
sage: R
Multivariate Polynomial Ring in x, y over Rational Field
Are you sure about that polynomial `f` ? because the following seems weird:
sage: f.subs(u=20)
0Tue, 19 May 2020 14:12:22 +0200https://ask.sagemath.org/question/51443/problem-with-implicit_plot/?answer=51460#post-id-51460Comment by Sébastien for <p>It seems that the polynomial f itself does not contain the good information.</p>
<pre><code>sage: f.parent()
Multivariate Polynomial Ring in u, v over Real Field with 500 bits of precision
</code></pre>
<p>Because when I draw its roots with:</p>
<pre><code>def roots_of_f(start,stop,step):
u_range = srange(start, stop, step)
L = [(u,v) for u in u_range
for v in f.subs(u=u).univariate_polynomial().roots(multiplicities=False)
if (u,v) in P]
return L
</code></pre>
<p>I get (<code>points</code> is a list in the code and overwrites the sage function <code>points</code>, so I am using <code>point</code> below):</p>
<pre><code>sage: point(roots_of_f(1,20,.05)) + plot_np
</code></pre>
<p><img alt="image description" src="/upfiles/15898901611271885.png"></p>
<p>Replacing <code>u</code> and <code>v</code> by <code>x</code> and <code>y</code> to use the ring over the rational field seems to give the same thing.</p>
<pre><code>sage: R
Multivariate Polynomial Ring in x, y over Rational Field
</code></pre>
<p>Are you sure about that polynomial <code>f</code> ? because the following seems weird:</p>
<pre><code>sage: f.subs(u=20)
0
</code></pre>
https://ask.sagemath.org/question/51443/problem-with-implicit_plot/?comment=51468#post-id-51468I replaced the function `draw_roots_of_f` in the answer by a function `roots_of_f` to return a list of roots. This allows me to verify the output of the roots finding algorithm. The following looks ok:
sage: L = roots_of_f(4,5,1)
sage: L
[(4, 5), (4, 4), (4, 3), (4, 2), (4, 1), (4, 0)]
sage: max(abs(f(uv)) for uv in L)
0
But, I agree that the following is weird:
sage: L = roots_of_f(4,5,.5)
sage: max(abs(f(uv)) for uv in L)
8.89546321071247e49
In fact the roots returned by `roots_of_f` are not precise enough:
sage: f(4,2)
0
sage: f(4,2.00000015504052)
8.65952734939081e46Wed, 20 May 2020 09:19:23 +0200https://ask.sagemath.org/question/51443/problem-with-implicit_plot/?comment=51468#post-id-51468Comment by jga for <p>It seems that the polynomial f itself does not contain the good information.</p>
<pre><code>sage: f.parent()
Multivariate Polynomial Ring in u, v over Real Field with 500 bits of precision
</code></pre>
<p>Because when I draw its roots with:</p>
<pre><code>def roots_of_f(start,stop,step):
u_range = srange(start, stop, step)
L = [(u,v) for u in u_range
for v in f.subs(u=u).univariate_polynomial().roots(multiplicities=False)
if (u,v) in P]
return L
</code></pre>
<p>I get (<code>points</code> is a list in the code and overwrites the sage function <code>points</code>, so I am using <code>point</code> below):</p>
<pre><code>sage: point(roots_of_f(1,20,.05)) + plot_np
</code></pre>
<p><img alt="image description" src="/upfiles/15898901611271885.png"></p>
<p>Replacing <code>u</code> and <code>v</code> by <code>x</code> and <code>y</code> to use the ring over the rational field seems to give the same thing.</p>
<pre><code>sage: R
Multivariate Polynomial Ring in x, y over Rational Field
</code></pre>
<p>Are you sure about that polynomial <code>f</code> ? because the following seems weird:</p>
<pre><code>sage: f.subs(u=20)
0
</code></pre>
https://ask.sagemath.org/question/51443/problem-with-implicit_plot/?comment=51466#post-id-51466Thanks for your reply! The polynomials are not irreducible and contain the 6 vertical lines x=19,20,21,22,23,24. I've added two lines in the original code to get rid of them. From what you wrote it seems like the issue goes all the way down to the "roots" algorithm. I gather this method relies on "Brent's method", which requires some kind of regularity from the function. Presumably these polynomials don't satisfy the required hypotheses. I will try to see if manually implementing some other root finding algorithms could help.Tue, 19 May 2020 23:18:13 +0200https://ask.sagemath.org/question/51443/problem-with-implicit_plot/?comment=51466#post-id-51466Comment by jga for <p>It seems that the polynomial f itself does not contain the good information.</p>
<pre><code>sage: f.parent()
Multivariate Polynomial Ring in u, v over Real Field with 500 bits of precision
</code></pre>
<p>Because when I draw its roots with:</p>
<pre><code>def roots_of_f(start,stop,step):
u_range = srange(start, stop, step)
L = [(u,v) for u in u_range
for v in f.subs(u=u).univariate_polynomial().roots(multiplicities=False)
if (u,v) in P]
return L
</code></pre>
<p>I get (<code>points</code> is a list in the code and overwrites the sage function <code>points</code>, so I am using <code>point</code> below):</p>
<pre><code>sage: point(roots_of_f(1,20,.05)) + plot_np
</code></pre>
<p><img alt="image description" src="/upfiles/15898901611271885.png"></p>
<p>Replacing <code>u</code> and <code>v</code> by <code>x</code> and <code>y</code> to use the ring over the rational field seems to give the same thing.</p>
<pre><code>sage: R
Multivariate Polynomial Ring in x, y over Rational Field
</code></pre>
<p>Are you sure about that polynomial <code>f</code> ? because the following seems weird:</p>
<pre><code>sage: f.subs(u=20)
0
</code></pre>
https://ask.sagemath.org/question/51443/problem-with-implicit_plot/?comment=51469#post-id-51469Thank you, @Sébastien. I think what we're seeing is just a more contrived version of this: https://en.wikipedia.org/wiki/Wilkinson%27s_polynomial . Your latest experiments also point in this direction. I also did some experiments on my side, please see the original post for the code and a brief summary. I'm also accepting this answer because it provided a solid explanation as to why every method fails. Specially the last comment showing the instability.Wed, 20 May 2020 11:35:19 +0200https://ask.sagemath.org/question/51443/problem-with-implicit_plot/?comment=51469#post-id-51469