Ask Your Question
3

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

asked 2018-04-04 11:46:06 +0100

Günter Rote gravatar image

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 flag offensive close merge delete

Comments

1

Welcome to Ask Sage! Thank you for your question!

slelievre gravatar imageslelievre ( 2018-04-04 16:19:32 +0100 )edit

1 Answer

Sort by » oldest newest most voted
2

answered 2018-04-04 16:18:56 +0100

slelievre gravatar image

updated 2018-04-04 16:27:57 +0100

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

Comments

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))
Sébastien gravatar imageSébastien ( 2018-04-04 23:04:24 +0100 )edit

@Sébastien -- I agree!

slelievre gravatar imageslelievre ( 2018-04-05 19:36:39 +0100 )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: 2018-04-04 11:46:06 +0100

Seen: 572 times

Last updated: Apr 04 '18