Plotting question

plot
complex

asked 2011-01-05 17:20:55 +0200

mhfrey
157 ●5 ●12 ●22

updated 2011-05-12 16:16:46 +0200

Kelvin Li
503 ●5 ●12 ●17

I am trying to plot Sage calculation listed below. I get the following error:

verbose 0 (3989: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 200 points. verbose 0 (3989: plot.py, generate_plot_points) Last error message: 'unable to simplify to float approximation'

What am I doing wrong?

Thank you for your help.

Gather all for Full H off resonance

Ix=1/2*matrix(SR,2,2,[[0,1],[1,0]])
Iy=1/(2*I)*matrix(SR,2,2,[[0,1],[-1,0]])
Iz=1/2*matrix(SR,2,2,[[1,0],[0,-1]])
Al=vector(SR,[1,0])
Be=vector(SR,[0,1])

g = 26.7522128e7 # 1H Gamma rad/s/T
rf_field = 100 # B1 rf field in kHz
Brf = rf_field*1e3*2*pi/g*2 # B1 in T for Wnut/2pi= 100 kHz
Thetarf = pi/2 # angle between BO and Brf

Wnut = (1/2*g*Brf*sin(Thetarf))
t360 = 1/((1/2*g*Brf*sin(Thetarf))/(2*pi))# 360 pulse
tip = 180 # Tip angle
tp = tip/360*t360
Bo = 500e6*2*pi/g# T Spectrometer Field Strength
wref = g*Bo # spectrometer reference frequency for Protons, rad/S
Bo.n(),(wref/(2*pi)).simplify()/1e6;(Wnut/(2*pi)).simplify()/1e3,tp.simplify()*1e6

Ratio = var ('Ratio')# Wnut/O0, Ratio creates offset, 0=on res
Boo=(Wnut*Ratio+wref)/g
woo = g*Boo
O0=woo-wref

Phip = var('Phip') # Phip = phase of rf pulse
Hrfrf = (O0)*Iz + Wnut*(Ix*cos(Phip)+Iy*sin(Phip))
Pabrf=((abs(Be*exp(-i*Hrfrf(Phip=0)*tp)*Al))^2)



for Ratio in srange (-1,1,0.25,include_endpoint=True):
    Pabrf(Ratio=Ratio).n(digits=5) #(pi)x on Alpha off resonance Full H

plot (Pabrf,(Ratio, -1,1))

edit retag flag offensive close merge delete

Comments

looks like nmr to me

Evgeny ( 2011-01-05 17:41:52 +0200 )edit

two quick tips that don't address your actual problem - you can use "var('Ratio')" instead of "Ratio = var('Ratio')", and the loop in the second and third last lines is unnecessary if you want to use plot()

Felix Lawrence ( 2011-01-05 18:51:30 +0200 )edit

add a comment

1

answered 2011-01-06 10:58:03 +0200

mhfrey
157 ●5 ●12 ●22

updated 2011-01-06 12:04:40 +0200

I have been able to get most things working with:

list_plot([(r,Pabrf(Ratio=r).n()) for r in srange(-1, 1, 0.1, include_endpoint=True)], plotjoined=True, frame=true).show(ymin=0,axes_labels=('$\Omega _{0}/\omega _{nut}$',"$ P _{a\Rightarrow b}$"))

I have an additional questions about controlling plot features:

How do I turn off or turn on the lines at the origin(o)?
How do I turn off the minor ticks?
How do I get the \rightarrow to display instead of the \Rightarrow? Only the \Rightarrow displays properly.

Thank you again for your help.

edit flag offensive delete link

Comments

Not really sure about any of these. Have a look through the plotting documentation http://sagemath.org/doc/reference/sage/plot/plot.html . If you want lots of control over your graphs then it might be worth checking out matplotlib/pylab, which is included with Sage but has separate documentation. Personally I've found sage's plotting functions convenient for quick and dirty figures, but use matplotlib when I want finesse.

Felix Lawrence ( 2011-01-07 01:53:47 +0200 )edit

add a comment

1

answered 2011-01-05 18:58:06 +0200

jk
43 ●1 ●1 ●7

'unable to simplify to float approximation'

Came here because I have a similar problem

Apparently the plotting command provided by Sage evaluate expressions by coercing to float, i.e. plot calls

float(Pabrf)

This fails on your expression because it's too complicated.

Can someone explain why Sage doesn't evaluate expressions numerically when plotting? What would be a straightforward why to do this?

edit flag offensive delete link

add a comment

1

answered 2011-01-05 18:59:35 +0200

Felix Lawrence
624 ●8 ●11 ●18

One solution is to use list_plot instead of plot. plot takes a function, whereas list_plot accepts a list of (x,y) coordinates to plot. list_plot therefore gives you a bit more control and flexibility. If you replace your last three lines by:

Pabrfs= []
for r in srange (-1,1,0.25,include_endpoint=True):
    #(pi)x on Alpha off resonance Full H
    Pabrfs.append((r,Pabrf(Ratio=r).n(digits=5)))

list_plot (Pabrfs)

then it should work. Alternatively you can make this a bit shorter using a list comprehension:

#(pi)x on Alpha off resonance Full H
Pabrfs = [(r,Pabrf(Ratio=r).n()) for r in srange(-1,1,0.25,include_endpoint=True)]

list_plot(Pabrfs)

edit flag offensive delete link

add a comment

1

answered 2011-05-24 19:03:33 +0200

nate_iverson
13 ●1 ●2 ●6

Try

def Pabrf:
    return N((abs(Be*exp(-i*Hrfrf(Phip=0)*tp)*Al))^2)
plot(Pabrf)

edit flag offensive delete link

Comments

you could also use a lambda function

niles ( 2011-05-25 02:46:42 +0200 )edit

add a comment

0

answered 2020-07-05 10:47:53 +0200

ortollj
1319 ●17 ●55 ●91

updated 2020-07-05 10:52:15 +0200

the first code in the question does not work for me with Ubuntu 18.04 notebook SageMath 9.1, but I noticed that if we ask for the graph of a function outside its domain of definition it causes a verbose and moreover the graph refuses to extend at the requested values

r=1
C_d(x)= -sqrt(r^2-x^2)
gf=plot(C_d,(x,-2,2),color='black') #produce verbose and do not span
#gf=plot(C_d,(x,-1,1),color='black') # no verbose 

show(gf,aspect_ratio=1)