How do I plot a real function whose computation involves complex intermediate results?

I have an algebraic expressing ev for an eigenvalue of a 4x4 matrix and I want to plot it. The expression has square roots and cube roots, and sometimes there are complex intermediate results although the final result is real. With float, the imaginary part will be non-zero due to numeric errors, but I can just ignore it. Still, when I try to plot this function, I get error messages.

var("u")
ev = (1/2*u + 1/2*sqrt(2/3*u^2 - 1/9*(u^4 - 18*u)/(1/27*u^6 - u^3 +
1/6*sqrt(-32/3*u^6 + 252*u^3 + 9) + 1/2)^(1/3) - (1/27*u^6 - u^3 +
1/6*sqrt(-32/3*u^6 + 252*u^3 + 9) + 1/2)^(1/3) + 6/sqrt((u^4 +
3*(1/27*u^6 - u^3 + 1/6*sqrt(-32/3*u^6 + 252*u^3 + 9) + 1/2)^(1/3)*u^2 -
18*u + 9*(1/27*u^6 - u^3 + 1/6*sqrt(-32/3*u^6 + 252*u^3 + 9) +
1/2)^(2/3))/(1/27*u^6 - u^3 + 1/6*sqrt(-32/3*u^6 + 252*u^3 + 9) +
1/2)^(1/3))) + 1/6*sqrt((u^4 + 3*(1/27*u^6 - u^3 + 1/6*sqrt(-32/3*u^6 +
252*u^3 + 9) + 1/2)^(1/3)*u^2 - 18*u + 9*(1/27*u^6 - u^3 +
1/6*sqrt(-32/3*u^6 + 252*u^3 + 9) + 1/2)^(2/3))/(1/27*u^6 - u^3 +
1/6*sqrt(-32/3*u^6 + 252*u^3 + 9) + 1/2)^(1/3)))

print ev.subs(u=2.6)
print ev.subs(u=3.0)
print ev.subs(u=3.0).real()
plot(ev.subs(u=x).real(),(x,2.6,3))


Plotting will draw only from 2.6 to 2.62 and then abort with the following error message

verbose 0 (3749: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 189 points.
verbose 0 (3749: plot.py, generate_plot_points) Last error message: 'math domain error'


I tried this in the sage cell server. (On my outdated office computer installation, sage 7.4, I get a different error "AssertionError".)

edit retag close merge delete

1

( 2018-04-04 09:19:32 -0500 )edit

Sort by » oldest newest most voted

Since the imaginary parts you get are just tiny error terms, the best is probably to ignore them by selecting the real part, as you understood.

The problem you are getting is that Sage's floating point numbers don't behave exactly like floats, in particular when raising negative numbers to fractional powers.

The "plot" command uses Python floats, not Sage "RealNumber"s.

So the trick with your definition of ev, would be to plot as follows:

sage: plot(lambda x: ev.subs(u=RR(x)).real(), (2.6, 3))


To pinpoint where the error occurs, observe the difference between:

sage: (-1.0)^(1/3)
0.500000000000000 + 0.866025403784439*I

sage: RR(-1)^(1/3)
0.500000000000000 + 0.866025403784439*I

sage: RealNumber(-1)^(1/3)
0.500000000000000 + 0.866025403784439*I


and

sage: float(-1)^(1/3)
Traceback (most recent call last)
...
ValueError: negative number cannot be raised to a fractional power

sage: RDF(-1)^(1/3)
Traceback (most recent call last)
...
ValueError: negative number cannot be raised to a fractional power

more

2

In this case, I would define a plain Python function instead of a symbolic expression :

sage: def ev(u):
....:     u = RR(u)
....:     ev = (1/2*u + 1/2*sqrt(2/3*u^2 - 1/9*(u^4 -  ....  + 1/2)^(1/3)))
....:     return ev.real()


So that the last plot command looks simple :

sage: plot(ev, (2.6, 3))

( 2018-04-04 16:04:24 -0500 )edit

@Sébastien -- I agree!

( 2018-04-05 12:36:39 -0500 )edit