Ask Your Question
0

Supposedly real function returning a complex result

asked 2017-06-16 17:20:31 -0500

Tony-64 gravatar image

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?

edit retag flag offensive close merge delete

2 answers

Sort by » oldest newest most voted
2

answered 2017-06-17 14:03:44 -0500

mforets gravatar image

here is a solution using fast_callable:

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!

edit flag offensive delete link more

Comments

Thanks 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.

Tony-64 gravatar imageTony-64 ( 2017-06-19 15:53:22 -0500 )edit
1

answered 2017-06-18 07:08:24 -0500

dan_fulea gravatar image

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.

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
edit flag offensive delete link more

Comments

Nice 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?

Tony-64 gravatar imageTony-64 ( 2017-06-19 16:09:17 -0500 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2017-06-16 17:20:31 -0500

Seen: 64 times

Last updated: Jun 18