Ask Your Question
1

Numerical integration and plot failing

asked 2019-06-06 02:30:11 -0600

psniady gravatar image

updated 2019-06-30 12:45:42 -0600

slelievre gravatar image

I have a function which can be efficiently calculated on some interval on the real line and if I feed Sage with concrete arguments, Sage gives concrete numerical values. Nevertheless, if I try to integrate it numerically (or plot it), something goes terribly wrong and Sage is not able to calculate the values.

In my case I get the error message

negative number cannot be raised to a fractional power

The function itself is a solution to the cubic equation so it involves square root and a cubic root.

is it some bug in sage and if there is some workaround?

MINIMAL NON_WORKING EXAMPLE

y = var('y')
z = var('z');

cauchy = y.substitute(solve(z*y^3 +y^2 - 2*z*y+2, y)[0])

def F(u):
    return arg(cauchy).substitute(z=u).n()

numerical_integral(F, 1, 2)
edit retag flag offensive close merge delete

Comments

Indeed :

sage: var("y,z")
(y, z)
sage: cauchy(z)=solve(z*y^3 +y^2 - 2*z*y+2 ,y)[0].rhs()
sage: def foo(u):return(arg(cauchy(u)).n())
sage: [foo(t) for t in (1,1.1..2)]
[0.828222717321238,
## Snip...
 0.377665318352514]

Therefore, this function can be evaluated. But :

sage: plot(foo,(1,2))
verbose 0 (3635: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 200 points.
verbose 0 (3635: plot.py, generate_plot_points) Last error message: 'negative number cannot be raised to a fractional power'
Launched png viewer for Graphics object consisting of 0 graphics primitives

... not in plot. Trying to plot symbolically yelds :

TypeError: Cannot evaluate symbolic expression to a numeric value.

A bug, IMHO

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2019-06-06 03:44:05 -0600 )edit

4 answers

Sort by » oldest newest most voted
2

answered 2019-06-30 20:06:30 -0600

Here’s a workaround for plotting:

list_plot( [[t, F(t)] for t in (1,1.1..2)], plotjoined=True )

with a live example. For integration you could likewise evaluate the function at discrete points, add all interior evaluations plus half of the values at the endpoints and multiply by the step size (i.e., a discrete integral).

edit flag offensive delete link more
1

answered 2019-07-01 18:36:27 -0600

Juanjo gravatar image

A different approach that works even in a Python 2 based Sage:

sage: var("y,z");
sage: cauchy(z) = y.subs(solve(z*y^3 + y^2 - 2*z*y + 2, y)[0])
sage: F(z) = arctan2(imag(cauchy(z)),real(cauchy(z)))
sage: numerical_integral(F, 1, 2)
(0.6116545934235456, 6.790730127365591e-15)
sage: plot(F(z),(z,1,2))

image description

The plot is really slow. It may take more than 20 seconds. See this SageMath cell.

edit flag offensive delete link more
0

answered 2019-07-06 10:27:55 -0600

Emmanuel Charpentier gravatar image

updated 2019-07-06 10:31:10 -0600

NOTE : This has replaced, (not complemented) my original answer concluding that this was a (Python) bug. Apologies for the confusion...

As pointed out by Serge Lelièvre, this seems a limitation of Python 2 floats, which does not exist in Python 3, where one can do :

sage: def F(u): return arg((u-3)^(1/3))
sage: F(1)
arg((-2)^(1/3))
sage: F(1).n()
1.04719755119660
sage: F(2).n()
1.04719755119660
sage: plot(F,(1,2))
Launched png viewer for Graphics object consisting of 1 graphics primitive
sage: integrate(F(x),x,1,2)
1/3*pi

See this sage-devel thread for further illumination...

In the original example:

sage: cauchy(z)=y.substitute(solve(z*y^3 +y^2 - 2*z*y+2, y)[0])
sage: def G(u): return arg(cauchy(u))
sage: numerical_integral(lambda u:G(u), 1, 2)
(0.6116545934235456, 6.790730127365591e-15)
edit flag offensive delete link more
0

answered 2019-06-06 04:18:50 -0600

slelievre gravatar image

updated 2019-06-06 14:04:11 -0600

This works for me in (Python 3 based) Sage 8.8.beta7:

sage: y, z = var('y z')
sage: cauchy = y.substitute(solve(z*y^3 + y^2 - 2*z*y + 2, y)[0])
sage: def F(u):
....:     return arg(cauchy).substitute(z=u).n()
....:
sage: numerical_integral(F, 1, 2)
(0.6116545934235456, 6.790730127365591e-15)

Here are links to Jupyter notebooks on CoCalc showing

This also fails with Py2-based SageMath 8.8.beta7.

So the "easy" solution for you is: switch to Python3-based Sage!

For this, follow the instructions to install from source, paying attention to the discussion about Python 3.

Or if you want to turn your existing SageMath into a Python3-based one, you could open a terminal, cd to your Sage directory, and type

make distclean
make configure
./configure --with-python=3
make

which will work if you have all the prerequisites.

edit flag offensive delete link more

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: 2019-06-06 02:30:11 -0600

Seen: 135 times

Last updated: Jul 06