ASKSAGE: Sage Q&A Forum - Individual question feedhttp://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Mon, 19 Jun 2017 16:09:17 -0500Supposedly real function returning a complex resulthttp://ask.sagemath.org/question/37979/supposedly-real-function-returning-a-complex-result/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?Fri, 16 Jun 2017 17:20:31 -0500http://ask.sagemath.org/question/37979/supposedly-real-function-returning-a-complex-result/Answer by dan_fulea for <p>Hi all,</p>
<p>I have this function:</p>
<pre><code>f(x, y) = 100*abs(x/(245*x^4/(20*x^2 - I*x*(20*y + 1) - 20) - 100*x^2 + I*x + 100))
</code></pre>
<p>and need to create a plot showing where its derivative with respect to x is negative:</p>
<pre><code>region_plot(diff(f, x) < 0, (x, 0.7, 1.7), (y, 0, 1))
</code></pre>
<p>Unfortunately, region_plot fails with a TypeError: unable to coerce to a real number.</p>
<p>Apparently, rounding errors cause the derivative to have a small complex part. Here is some sample code to demonstrate the effect:</p>
<pre><code>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
</code></pre>
<p>and its output:</p>
<pre><code>(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
</code></pre>
<p>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.</p>
<p>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.</p>
<p>Why doesn’t this work-around work? Is it a Sage bug?</p>
<p>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</p>
<pre><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
</code></pre>
<p>and its output:</p>
<pre><code>(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
</code></pre>
<p>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.</p>
<p>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.</p>
<p>However, this latter work-around is much worse than the above mentioned work-around: the region_plot takes <strong>much</strong> 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.</p>
<p>Here is the code to create the plots:</p>
<pre><code>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()
</code></pre>
<p>Anyone having suggestions for a better solution?</p>
http://ask.sagemath.org/question/37979/supposedly-real-function-returning-a-complex-result/?answer=37999#post-id-37999The following is a solution based roughly on the algebraic idea, that when dealing with rational functions and the signs of their derivatives we can always do computations so far, that the sign of a polynomial has to be computed. Polynomials are quick enough to compute, and this is the last wall of a quick improvement.
So let us start with the given function $f$, compute $|f|^2$ as in the question and the other answer and write it in the form
$$ |f|^2 = \frac AB\ . $$
var('x, y')
assume( x, 'real' )
assume( y, 'real' )
f = ( 245*x^4 / (20*x^2 - I*x*(20*y + 1) - 20) - 100*x^2 + I*x + 100 )
ref = real(f)
imf = imag(f)
print ( "...ref and imf have the same denominator? %s"
% bool( ref.denominator() == imf.denominator() ) )
absf2 = ( ref^2 + imf^2 ).factor()
A = absf2. numerator()
B = absf2.denominator()
print ( "...just a check: Is abs(f)^2 == A/B ? %s"
% bool( ref^2 + imf^2 == A/B ) )
(The check i for my nerves, sometimes a factor like $-2017$ gets lost. A $-1$ is also enough to invert the picture.)
The computed values of $A$ (manually rearranged) and $B$ are:
sage: A
3080025*x^8 + 4000000*x^6*y^2 + 409800*x^6*y - 14029110*x^6
- 7999600*x^4*y^2 - 799960*x^4*y + 22999201*x^4 + 4000000*x^2*y^2
+ 400000*x^2*y - 15989600*x^2 + 4000000
sage: B
400*x^4 + 400*x^2*y^2 + 40*x^2*y - 799*x^2 + 400
Note that the denominator $B$ is $\ge0$, since it is
sage: ( (20*x^2 - I*x*(20*y + 1) - 20) * (20*x^2 + I*x*(20*y + 1) - 20) ).expand()
400*x^4 + 400*x^2*y^2 + 40*x^2*y - 799*x^2 + 400
but we will not need this. (This was just an other check.)
Now we need the differential of $\sqrt{A/B}$. This is
$$ \left(\ \sqrt{\frac AB}\ \right)' =\frac 1{2\sqrt{\dots}}\frac 1{B^2}(A'B-AB')\ . $$
So we need the sign of the polynomial
$$ A'B-AB'\ . $$
The code for its plot on the square $[0,2]^2$ is then given by the further lines:
G = ( diff( A, x ) * B - A * diff( B, x ) ).factor()
P = region_plot( G < 0
, (x, 0, 2)
, (y, 0, 2)
, plot_points = 100
, aspect_ratio = 1 )
P . show()
N.B.
Despite of my effort to control the sign, i've got the colors inverted...
One more sign check:
sage: diff( real(f)^2 + imag(f)^2 , x ).subs(x=0.1).subs(y=0.1)
-3969.59625416928
sage: G.subs(x=0.1).subs(y=0.1)
-6.10388687137070e8
Sun, 18 Jun 2017 07:08:24 -0500http://ask.sagemath.org/question/37979/supposedly-real-function-returning-a-complex-result/?answer=37999#post-id-37999Comment by Tony-64 for <p>The following is a solution based roughly on the algebraic idea, that when dealing with rational functions and the signs of their derivatives we can always do computations so far, that the sign of a polynomial has to be computed. Polynomials are quick enough to compute, and this is the last wall of a quick improvement. </p>
<p>So let us start with the given function $f$, compute $|f|^2$ as in the question and the other answer and write it in the form
$$ |f|^2 = \frac AB\ . $$</p>
<pre><code>var('x, y')
assume( x, 'real' )
assume( y, 'real' )
f = ( 245*x^4 / (20*x^2 - I*x*(20*y + 1) - 20) - 100*x^2 + I*x + 100 )
ref = real(f)
imf = imag(f)
print ( "...ref and imf have the same denominator? %s"
% bool( ref.denominator() == imf.denominator() ) )
absf2 = ( ref^2 + imf^2 ).factor()
A = absf2. numerator()
B = absf2.denominator()
print ( "...just a check: Is abs(f)^2 == A/B ? %s"
% bool( ref^2 + imf^2 == A/B ) )
</code></pre>
<p>(The check i for my nerves, sometimes a factor like $-2017$ gets lost. A $-1$ is also enough to invert the picture.)</p>
<p>The computed values of $A$ (manually rearranged) and $B$ are:</p>
<pre><code>sage: A
3080025*x^8 + 4000000*x^6*y^2 + 409800*x^6*y - 14029110*x^6
- 7999600*x^4*y^2 - 799960*x^4*y + 22999201*x^4 + 4000000*x^2*y^2
+ 400000*x^2*y - 15989600*x^2 + 4000000
sage: B
400*x^4 + 400*x^2*y^2 + 40*x^2*y - 799*x^2 + 400
</code></pre>
<p>Note that the denominator $B$ is $\ge0$, since it is</p>
<pre><code>sage: ( (20*x^2 - I*x*(20*y + 1) - 20) * (20*x^2 + I*x*(20*y + 1) - 20) ).expand()
400*x^4 + 400*x^2*y^2 + 40*x^2*y - 799*x^2 + 400
</code></pre>
<p>but we will not need this. (This was just an other check.)</p>
<p>Now we need the differential of $\sqrt{A/B}$. This is
$$ \left(\ \sqrt{\frac AB}\ \right)' =\frac 1{2\sqrt{\dots}}\frac 1{B^2}(A'B-AB')\ . $$
So we need the sign of the polynomial
$$ A'B-AB'\ . $$
The code for its plot on the square $[0,2]^2$ is then given by the further lines:</p>
<pre><code>G = ( diff( A, x ) * B - A * diff( B, x ) ).factor()
P = region_plot( G < 0
, (x, 0, 2)
, (y, 0, 2)
, plot_points = 100
, aspect_ratio = 1 )
P . show()
</code></pre>
<p>N.B.</p>
<p>Despite of my effort to control the sign, i've got the colors inverted...</p>
<p>One more sign check:</p>
<pre><code>sage: diff( real(f)^2 + imag(f)^2 , x ).subs(x=0.1).subs(y=0.1)
-3969.59625416928
sage: G.subs(x=0.1).subs(y=0.1)
-6.10388687137070e8
</code></pre>
http://ask.sagemath.org/question/37979/supposedly-real-function-returning-a-complex-result/?comment=38029#post-id-38029Nice trick! It certainly gives me some food for further thought about how to optimize the rest of the application. Obviously, the function f represents only a small part of the stuff I have to visualize ;-)
Anyway, with respect to the question about region_plot(real(g) < 0) failing with TypeError: unable to coerce to a real number, is this a bug? If yes, how and where should I report it?Mon, 19 Jun 2017 16:09:17 -0500http://ask.sagemath.org/question/37979/supposedly-real-function-returning-a-complex-result/?comment=38029#post-id-38029Answer by mforets for <p>Hi all,</p>
<p>I have this function:</p>
<pre><code>f(x, y) = 100*abs(x/(245*x^4/(20*x^2 - I*x*(20*y + 1) - 20) - 100*x^2 + I*x + 100))
</code></pre>
<p>and need to create a plot showing where its derivative with respect to x is negative:</p>
<pre><code>region_plot(diff(f, x) < 0, (x, 0.7, 1.7), (y, 0, 1))
</code></pre>
<p>Unfortunately, region_plot fails with a TypeError: unable to coerce to a real number.</p>
<p>Apparently, rounding errors cause the derivative to have a small complex part. Here is some sample code to demonstrate the effect:</p>
<pre><code>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
</code></pre>
<p>and its output:</p>
<pre><code>(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
</code></pre>
<p>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.</p>
<p>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.</p>
<p>Why doesn’t this work-around work? Is it a Sage bug?</p>
<p>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</p>
<pre><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
</code></pre>
<p>and its output:</p>
<pre><code>(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
</code></pre>
<p>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.</p>
<p>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.</p>
<p>However, this latter work-around is much worse than the above mentioned work-around: the region_plot takes <strong>much</strong> 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.</p>
<p>Here is the code to create the plots:</p>
<pre><code>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()
</code></pre>
<p>Anyone having suggestions for a better solution?</p>
http://ask.sagemath.org/question/37979/supposedly-real-function-returning-a-complex-result/?answer=37989#post-id-37989here is a solution using [fast_callable](http://doc.sagemath.org/html/en/reference/misc/sage/ext/fast_callable.html):
sage: (x, y, z) = var('x, y, z')
sage: assume(x, 'real', y, 'real')
sage: from sage.ext.fast_callable import fast_callable
sage: cabs = fast_callable(sqrt(real(z)^2 + imag(z)^2), vars=[z])
sage: cg = fast_callable(diff(100*cabs(x/(245*x^4/(20*x^2 - I*x*(20*y + 1) - 20) - 100*x^2 + I*x + 100)), x), vars=[x, y])
sage: region_plot(cg(x, y) < 0, (x, 0.7, 1.7), (y, 0, 1), plot_points=100)
it seems to work and in a reasonable time, but someone more knowledgeable may have a better approach!
Sat, 17 Jun 2017 14:03:44 -0500http://ask.sagemath.org/question/37979/supposedly-real-function-returning-a-complex-result/?answer=37989#post-id-37989Comment by Tony-64 for <p>here is a solution using <a href="http://doc.sagemath.org/html/en/reference/misc/sage/ext/fast_callable.html">fast_callable</a>:</p>
<pre><code>sage: (x, y, z) = var('x, y, z')
sage: assume(x, 'real', y, 'real')
sage: from sage.ext.fast_callable import fast_callable
sage: cabs = fast_callable(sqrt(real(z)^2 + imag(z)^2), vars=[z])
sage: cg = fast_callable(diff(100*cabs(x/(245*x^4/(20*x^2 - I*x*(20*y + 1) - 20) - 100*x^2 + I*x + 100)), x), vars=[x, y])
sage: region_plot(cg(x, y) < 0, (x, 0.7, 1.7), (y, 0, 1), plot_points=100)
</code></pre>
<p>it seems to work and in a reasonable time, but someone more knowledgeable may have a better approach!</p>
http://ask.sagemath.org/question/37979/supposedly-real-function-returning-a-complex-result/?comment=38028#post-id-38028Thanks for the example. Whereas cg calculates about a factor of 5 faster on my machine, region_plot(cg) takes about the same time toe execute as region_plot(g). Apparently, region_plot() doesn't spend the bulk of its time doing the raw calculations.Mon, 19 Jun 2017 15:53:22 -0500http://ask.sagemath.org/question/37979/supposedly-real-function-returning-a-complex-result/?comment=38028#post-id-38028