Hi all,
I have this function:
f(x, y) = 100*abs(x/(245*x^4/(20*x^2 - I*x*(20*y + 1) - 20) - 100*x^2 + I*x + 100))
and need to create a plot showing where its derivative with respect to x is negative:
region_plot(diff(f, x) < 0, (x, 0.7, 1.7), (y, 0, 1))
Unfortunately, region_plot fails with a TypeError: unable to coerce to a real number.
Apparently, rounding errors cause the derivative to have a small complex part. Here is some sample code to demonstrate the effect:
var('x, y')
assume(x, 'real')
assume(y, 'real')
assume(x > 0)
assume(y > 0)
f(x, y) = 100*abs(x/(245*x^4/(20*x^2 - I*x*(20*y + 1) - 20) - 100*x^2 + I*x + 100))
g(x, y) = diff(f, x)
print f
print
print g
print
xv = 1.37
yv = 0.28
print "rational arguments: f(", QQ(xv), ",", QQ(yv), ") =", N(f(QQ(xv), QQ(yv)))
print "rational arguments: g(", QQ(xv), ",", QQ(yv), ") =", N(g(QQ(xv), QQ(yv)))
print
print "float arguments: f(", xv, ",", yv, ") =", f(xv, yv)
print "float arguments: g(", xv, ",", yv, ") =", g(xv, yv)
print
print "float arguments: real(f(", xv, ",", yv, ")) =", real(f(xv, yv))
print "float arguments: real(g(", xv, ",", yv, ")) =", real(g(xv, yv))
print
and its output:
(x, y) |--> 100*abs(x/(245*x^4/(20*x^2 - I*x*(20*y + 1) - 20) - 100*x^2 + I*x + 100))
(x, y) |--> 50*(x*((245*(40*x + 20*I*y + I)*x^4/(20*x^2 + I*x*(20*y + 1) - 20)^2 - 980*x^3/(20*x^2 + I*x*(20*y + 1) - 20) + 200*x + I)*x/(245*x^4/(20*x^2 + I*x*(20*y + 1) - 20) - 100*x^2 - I*x + 100)^2 + 1/(245*x^4/(20*x^2 + I*x*(20*y + 1) - 20) - 100*x^2 - I*x + 100))/(245*x^4/(20*x^2 - I*x*(20*y + 1) - 20) - 100*x^2 + I*x + 100) + x*((245*(40*x - 20*I*y - I)*x^4/(20*x^2 - I*x*(20*y + 1) - 20)^2 - 980*x^3/(20*x^2 - I*x*(20*y + 1) - 20) + 200*x - I)*x/(245*x^4/(20*x^2 - I*x*(20*y + 1) - 20) - 100*x^2 + I*x + 100)^2 + 1/(245*x^4/(20*x^2 - I*x*(20*y + 1) - 20) - 100*x^2 + I*x + 100))/(245*x^4/(20*x^2 + I*x*(20*y + 1) - 20) - 100*x^2 - I*x + 100))/abs(x/(245*x^4/(20*x^2 - I*x*(20*y + 1) - 20) - 100*x^2 + I*x + 100))
rational arguments: f( 137/100 , 7/25 ) = 2.57021481299952
rational arguments: g( 137/100 , 7/25 ) = -8.24860099120441
float arguments: f( 1.37000000000000 , 0.280000000000000 ) = 2.57021481299952
float arguments: g( 1.37000000000000 , 0.280000000000000 ) = -8.24860099120440 + 1.68733316297436e-15*I
float arguments: real(f( 1.37000000000000 , 0.280000000000000 )) = 2.57021481299952
float arguments: real(g( 1.37000000000000 , 0.280000000000000 )) = -8.24860099120440
Clearly, as long as the arguments x and y are rationals, the derivative of f (g in the sample code) yields a real number. However, when the arguments are given as floating point values, in many cases the result has a small imaginary component. By the way, f itself always yields a real number.
Whereas it is possible to take this component away by taking the real part, this work-around apparently does not help in the region_plot function: region_plot(real(g) < 0, (x, 0.7, 1.7), (y, 0, 1)) fails with the same TypeError: unable to coerce to a real number.
Why doesn’t this work-around work? Is it a Sage bug?
Anyway, the only feasible work-around that I have been able to find so far is to replace the abs function with my_abs(z) = sqrt(real(z)^2 + imag(z)^2). Here is the code
my_abs(z) = sqrt(real(z)^2 + imag(z)^2)
f(x, y) = 100*my_abs(x/(245*x^4/(20*x^2 - I*x*(20*y + 1) - 20) - 100*x^2 + I*x + 100))
g(x, y) = diff(f, x)
print f
print
print g
print
xv = 1.37
yv = 0.28
print "rational arguments: f(", QQ(xv), ",", QQ(yv), ") =", N(f(QQ(xv), QQ(yv)))
print "rational arguments: g(", QQ(xv), ",", QQ(yv), ") =", N(g(QQ(xv), QQ(yv)))
print
print "float arguments: f(", xv, ",", yv, ") =", f(xv, yv)
print "float arguments: g(", xv, ",", yv, ") =", g(xv, yv)
print
print "float arguments: real(f(", xv, ",", yv, ")) =", real(f(xv, yv))
print "float arguments: real(g(", xv, ",", yv, ")) =", real(g(xv, yv))
print
print "time required to calculate g:"
timeit('g(xv, yv)')
print
and its output:
(x, y) |--> 100*sqrt(10000*(49*x^7/((24010000*x^12/(400*x^4 + 400*x^2*y^2 + 40*x^2*y - 799*x^2 + 400)^2 + 2401000*x^10*y/(160000*x^8 + 320000*x^6*y^2 + 160000*x^4*y^4 + [... output truncated ...]
(x, y) |--> 200*(5000*(49*x^7/((24010000*x^12/(400*x^4 + 400*x^2*y^2 + 40*x^2*y - 799*x^2 + 400)^2 + 2401000*x^10*y/(160000*x^8 + 320000*x^6*y^2 + 160000*x^4*y^4 + [... output truncated ...]
rational arguments: f( 137/100 , 7/25 ) = 2.57021481299952
rational arguments: g( 137/100 , 7/25 ) = -8.24860099120441
float arguments: f( 1.37000000000000 , 0.280000000000000 ) = 2.57021481299952
float arguments: g( 1.37000000000000 , 0.280000000000000 ) = -8.24860099120455
float arguments: real(f( 1.37000000000000 , 0.280000000000000 )) = 2.57021481299952
float arguments: real(g( 1.37000000000000 , 0.280000000000000 )) = -8.24860099120455
time required to calculate g:
5 loops, best of 3: 43.5 ms per loop
However, this work-around unfortunately yields a substantial performance loss (the calculation of g using the built-in abs function takes about 20 ms). Since I have to calculate this function—and several even more complex ones—many times, this performance decrease by 50% does have a sever impact.
An alternative work-around that has been suggested elsewhere in this forum is to make use of the lambda construct, as in region_plot(lambda x, y: g(x, y) < 0, (x, 0.7, 1.7), (y, 0, 1)), where g is the based on the built-in abs function.
However, this latter work-around is much worse than the above mentioned work-around: the region_plot takes much longer than doing a region_plot without the lambda construct (but using the my_abs function); about 210 seconds versus 5 seconds on my machine (2012 iMac)! Moreover, the resulting plot is very coarse and grainy, rather than nicely smooth.
Here is the code to create the plots:
var('x, y')
assume(x, 'real')
assume(y, 'real')
assume(x > 0)
assume(y > 0)
f(x, y) = 100*abs(x/(245*x^4/(20*x^2 - I*x*(20*y + 1) - 20) - 100*x^2 + I*x + 100))
g(x, y) = diff(f, x)
# g based on the built-in abs function:
P = region_plot(lambda x, y: g(x, y) < 0, (x, 0.7, 1.7), (y, 0, 1), plot_points=100)
P.show()
my_abs(z) = sqrt(real(z)^2 + imag(z)^2)
f(x, y) = 100*my_abs(x/(245*x^4/(20*x^2 - I*x*(20*y + 1) - 20) - 100*x^2 + I*x + 100))
g(x, y) = diff(f, x)
# g based on the my_abs function
P = region_plot(g(x, y) < 0, (x, 0.7, 1.7), (y, 0, 1), plot_points=100)
P.show()
Anyone having suggestions for a better solution?
