ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Thu, 05 Apr 2018 19:36:39 +0200How do I plot a real function whose computation involves complex intermediate results?https://ask.sagemath.org/question/41875/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".)Wed, 04 Apr 2018 11:46:06 +0200https://ask.sagemath.org/question/41875/how-do-i-plot-a-real-function-whose-computation-involves-complex-intermediate-results/Comment by slelievre for <p>I have an algebraic expressing <strong>ev</strong> 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.</p>
<pre><code>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))
</code></pre>
<p>Plotting will draw only from 2.6 to 2.62 and then abort with the following error message</p>
<pre><code>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'
</code></pre>
<p>I tried this in the sage cell server. (On my outdated office computer installation, sage 7.4, I get a different error "AssertionError".)</p>
https://ask.sagemath.org/question/41875/how-do-i-plot-a-real-function-whose-computation-involves-complex-intermediate-results/?comment=41877#post-id-41877Welcome to Ask Sage! Thank you for your question!Wed, 04 Apr 2018 16:19:32 +0200https://ask.sagemath.org/question/41875/how-do-i-plot-a-real-function-whose-computation-involves-complex-intermediate-results/?comment=41877#post-id-41877Answer by slelievre for <p>I have an algebraic expressing <strong>ev</strong> 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.</p>
<pre><code>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))
</code></pre>
<p>Plotting will draw only from 2.6 to 2.62 and then abort with the following error message</p>
<pre><code>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'
</code></pre>
<p>I tried this in the sage cell server. (On my outdated office computer installation, sage 7.4, I get a different error "AssertionError".)</p>
https://ask.sagemath.org/question/41875/how-do-i-plot-a-real-function-whose-computation-involves-complex-intermediate-results/?answer=41876#post-id-41876Since 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 powerWed, 04 Apr 2018 16:18:56 +0200https://ask.sagemath.org/question/41875/how-do-i-plot-a-real-function-whose-computation-involves-complex-intermediate-results/?answer=41876#post-id-41876Comment by slelievre for <p>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.</p>
<p>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.</p>
<p>The "plot" command uses Python floats, not Sage "RealNumber"s.</p>
<p>So the trick with your definition of <code>ev</code>, would be to plot as follows:</p>
<pre><code>sage: plot(lambda x: ev.subs(u=RR(x)).real(), (2.6, 3))
</code></pre>
<p>To pinpoint where the error occurs, observe the difference between:</p>
<pre><code>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
</code></pre>
<p>and</p>
<pre><code>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
</code></pre>
https://ask.sagemath.org/question/41875/how-do-i-plot-a-real-function-whose-computation-involves-complex-intermediate-results/?comment=41883#post-id-41883@Sébastien -- I agree!Thu, 05 Apr 2018 19:36:39 +0200https://ask.sagemath.org/question/41875/how-do-i-plot-a-real-function-whose-computation-involves-complex-intermediate-results/?comment=41883#post-id-41883Comment by Sébastien for <p>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.</p>
<p>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.</p>
<p>The "plot" command uses Python floats, not Sage "RealNumber"s.</p>
<p>So the trick with your definition of <code>ev</code>, would be to plot as follows:</p>
<pre><code>sage: plot(lambda x: ev.subs(u=RR(x)).real(), (2.6, 3))
</code></pre>
<p>To pinpoint where the error occurs, observe the difference between:</p>
<pre><code>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
</code></pre>
<p>and</p>
<pre><code>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
</code></pre>
https://ask.sagemath.org/question/41875/how-do-i-plot-a-real-function-whose-computation-involves-complex-intermediate-results/?comment=41880#post-id-41880In 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))Wed, 04 Apr 2018 23:04:24 +0200https://ask.sagemath.org/question/41875/how-do-i-plot-a-real-function-whose-computation-involves-complex-intermediate-results/?comment=41880#post-id-41880