0

tengo un problema con éste código, el cual calcula los eigenestados para el oscilador armónico cuántico y lo evalúa en un tiempo dado. A su vez, intento graficar todos los subplots de estados para cada valor de n. Pero me aparece este error:

for i in xx: rePsi.eval(i)
AttributeError: 'Mul' object has no attribute 'eval'*

Este es mi código:

import numpy as np
import matplotlib.pylab as plt
import sympy as sym
from math import *

w=1
m=1
hbar=1
x = sym.symbols ('x')
a=-1
b=1
phi_n=0
psi=0
plt.ylim(-1,1 )
plt.xlim(-3,3)

print('give me number of states?')
nn = int(input())
if (nn<0):
    print 'must be >=0\n'
if (nn>=0):
    print('time to evaluate')
    t = float(input())
    if (t<0):
        print 'Time must be >=0'

print '   '

mr = float((m*w/(hbar*pi))**0.25)
gaussq = sym.exp(-m*w*x**2/(2*hbar))

fig = [0 for i in range(nn)]
fig1 = plt.figure(1)

xx = np.arange(-3., 3.,0.1)
for n in range (0,nn,2):
    H_root = sym.sqrt(1/(2**n*float(sym.factorial(n))))
    Hermite = (-1)**n*sym.exp(m*w*x**2/hbar)*sym.diff(sym.exp(-m*w*x**2/hbar),x,n)
    Energy = hbar*w*(n+0.5)
    phi_n = sym.simplify(H_root*mr*gaussq*Hermite)
    arg_int = sym.simplify(phi_n*gaussq*mr)
    c_n = sym.simplify(sym.integrate(arg_int, (x,a,b)))
    psi_n = c_n*phi_n*sym.exp(-1j*w*t*(n+0.5))
    psi_n_re = sym.re(psi_n)
    psi_n_im = sym.im(psi_n)
    psi += sym.simplify(psi_n)
    rePsi = sym.re(psi) 
    imPsi = sym.im(psi)
    print type (rePsi)
    print type (imPsi)
    for i in xx: rePsi.eval(i)
    for i in xx: imPsi.eval(i)
    list_rePsi= rePsi
    list_imPsi= imPsi
    plt.subplot(nn, 1, (n+2)/2)
    plt.plot(xx,list_repsi)
    plt.plot(xx,list_impsi)

plt.show()
print '\n'+str(psi)
1

1 respuesta 1

1

Efectivamente el método eval no existe para las expresiones SymPy, lo que debes hacer es substituir la variable (símbolo x) en todo caso. Para ello lo más simple es usar el métododo subs y si quieres el valor numérico evalf junto a subs. Esta substitución debes hacerla para cada valor de xx.

El código con algunas modificaciones, principalmente en la validación de la entrada y en la creación de las gráficas podría quedar así:

import numpy as np
import matplotlib.pylab as plt
import sympy as sym
from math import ceil


w = 1
m = 1

hbar = 1
x = sym.symbols ('x')
a = -1
b = 1
phi_n = 0
psi = 0
plt.ylim(-1, 1)
plt.xlim(-3, 3)



try:
    nn = int(raw_input('Give me number of states: '))
except ValueError:
    raise ValueError('Number of states must be integuer')
else:
    if nn < 0:
        raise ValueError('Number of states must be >= 0')

try:
    t = float(input('Time to evaluate: '))
except ValueError:
    raise ValueError('Number of states must be float')
else:
    if nn < 0:
        raise ValueError('Time must be >= 0')


mr = float((m * w / (hbar * sym.pi)) ** 0.25)
gaussq = sym.exp(-m * w * x ** 2 / (2 * hbar))

fig = plt.figure()
fig.subplots_adjust(hspace=0.7)
xx = np.arange(-3.0, 3.0, 0.1)

for n in xrange(0, nn, 2):
    h_root = np.sqrt(1 / (2 ** n * float(sym.factorial(n))))
    hermite = (-1 ** n * sym.exp(m * w * x ** 2 / hbar) *
              sym.diff(sym.exp(-m * w * x ** 2 / hbar), x, n))

    energy = hbar * w * (n + 0.5)
    phi_n = sym.simplify(h_root * mr * gaussq * hermite)
    arg_int = sym.simplify(phi_n * gaussq * mr)
    c_n = sym.simplify(sym.integrate(arg_int, (x, a, b)))
    psi_n = c_n * phi_n * sym.exp(-1j * w * t * (n + 0.5))
    psi_n_re = sym.re(psi_n)
    psi_n_im = sym.im(psi_n)
    psi += sym.simplify(psi_n)
    re_psi = sym.re(psi) 
    im_psi = sym.im(psi)

    list_re_psi= [re_psi.subs(x, i) for i in xx]
    list_im_psi= [im_psi.subs(x, i) for i in xx]

    ax = fig.add_subplot(ceil(nn / 2.0), 1, (n + 2) / 2)
    ax.plot(xx, list_re_psi)
    ax.plot(xx, list_im_psi)
    ax.set_title("n = {}".format(n))


plt.show()
print('\n' + str(psi))

introducir la descripción de la imagen aquí

Matplotlib no necesita que la expresión sea explícitamente evaluada, pero si queremos que las listas contengan los valores numéricos y no expresiones podemos usar evalf como se comentó antes:

list_im_psi= [im_psi.evalf(subs={x: i}) for i in xx]

Hay que tener en cuenta subs + evalf es lento, si bien en este caso el array xx es pequeño y el impacto es escaso, teniendo en cuenta que es la integración la parte más costosa con diferencia. Lo óptimo es usar sympy.lambdify y aplicar la función de forma vectorizada sobre el array:


re_psi_func = sym.lambdify(x, re_psi, "numpy")
list_re_psi = re_psi_func(xx)

el problema es que a menos en mi caso con SymPy 1.3 y Numpy 1.15.0 falla a la hora de convertir sympy.im a numpy.imag en el proceso, en lo que aparentemente es un bug.

Tu Respuesta

By clicking “Publica tu respuesta”, you agree to our terms of service and acknowledge you have read our privacy policy.

¿No es la respuesta que buscas? Examina otras preguntas con la etiqueta o formula tu propia pregunta.