Ask Your Question
0

Contour Plot not working with pseudo-piecewise function

asked 13 years ago

Jonathan.B gravatar image

Hi, I'm fairly new to Sage and don't know if this is a bug or just my own limited understanding. I want to plot a function 'f(g1, g2)' but only in some relevant area of the plot. All other points should display/evaluate to zero. This should work fine in Octave or numpy but Sage for some reason gives me an error-message which I fail to understand.

Function assignment like so:

def f(g1, g2):
if g1*g2 > 0.0:
    if abs(g2) < 1/abs(g1):
        return sqrt(0.5*((1+g1*g2)*abs(g1+g2))/(sqrt((1.0-g1*g2)^3)*g1*g2))
else:
    return 0

Then the Plot:

var('g1 g2')
contour_plot(f(g1,g2), (g1,-4,4), (g2,-4,4),plot_points=500, cmap='afmhot_r', colorbar=True, contours=srange(1,6,0.5))

and at last the error-message:

ValueError: zero-size array to ufunc.reduce without identity

Strangely, if I assign a single number to 'contours' (or leave it out to get the default value) the error-message does not show but instead I get an empty/white plot.

If I do not use a conditional/if-statement in my function-assignment the contour-plot shows just fine (except that it didn't do what I wanted). Does anyone know how to fix this or why it doesn't work or at least why I'm being silly for even trying it that way?

Thanks in advance

Preview: (hide)

1 Answer

Sort by » oldest newest most voted
3

answered 13 years ago

DSM gravatar image

updated 13 years ago

I can't reproduce the error (4.7.1), but I can see at least one problem.

The way the scoping works, f(g1, g2) is evaluated before the plot is made, so your plot will be zero (or maybe None) everywhere. Stick a "print g1,g2" inside the function to convince yourself this is true. You probably want "f" instead of "f(g1,g2)".

The following seems to work for me. Note that I've also changed it so that it returns 0 instead of None sometimes like your original code did-- if that's what you want, switch back.

var("g1 g2")
def f(g1, g2):
    if g1*g2 > 0.0:
        if abs(g2) < 1/abs(g1):
            return sqrt(0.5*((1+g1*g2)*abs(g1+g2))/(sqrt((1.0-g1*g2)^3)*g1*g2))
    return 0

z=contour_plot(f, (g1,-4,4), (g2,-4,4),plot_points=500, cmap='afmhot_r', colorbar=True, contours=srange(1,6,0.5))
show(z)

image description

Preview: (hide)
link

Comments

Awesome, thanks for the quick reply. I had never considered the scope being an issue here.

Jonathan.B gravatar imageJonathan.B ( 13 years ago )

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: 13 years ago

Seen: 1,248 times

Last updated: Oct 09 '11