Supposedly real function returning a complex result

Tony-64
91 ●3 ●7 ●13

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?

add a comment

2 Answers

Sort by » oldest newest most voted

answered 7 years ago

mforets
1413 ●9 ●27 ●44 http://marcelo-forets.fr/

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!

link

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 ( 7 years ago )

add a comment

answered 7 years ago

dan_fulea
5650 ●5 ●44 ●97

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

link

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 ( 7 years ago )

add a comment

Your Answer

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

Add Answer

Supposedly real function returning a complex result

2 Answers

Comments

Comments

Your Answer

Question Tools

Stats

Supposedly real function returning a complex result savecancel

2 Answers

Comments

Comments

Your Answer

Question Tools

Stats

Supposedly real function returning a complex result