Oscilador Armónico

Jorge Garza. Enero de 2014

La ecuación de Schroedinger asociada al oscilador armónico tiene la forma
-((ℏ^2)/2m)(((d^2)ψ)/(d(x^2)))+(k/2) (x^2)ψ=Eψ
La cual se puede escribir como
(((d^2)ψ)/(d(x^2)))-( m k/(ℏ^2)) (x^2)ψ=-(2m E/(ℏ^2))ψ

Usando el cambio de variable x=α y (y=x/α) se tiene
(1/(α^2))(((d^2)ψ)/(d(y^2)))-((α^2)m k/(ℏ^2)) (y^2)ψ=-(2m E/(ℏ^2))ψ
(((d^2)ψ)/(d(y^2)))-((α^4)m k/(ℏ^2)) (y^2)ψ=-(2m E (α^2)/(ℏ^2))ψ
Definiendo a ϵ=2m E (α^2)/(ℏ^2)
(((d^2)ψ)/(d(y^2)))-((α^4)m k/(ℏ^2)) (y^2)ψ=-ϵ ψ

y eligiendo a α tal que
(α^4)m k/(ℏ^2)=1
(α^4)=(ℏ^2)/(m k) con lo que
(α^2)=ℏ/(((m k))^(1/2)) 
como
ϵ=2m E (α^2)/(ℏ^2) entonces
ϵ=(2m E /(ℏ^2))(ℏ/(((m k))^(1/2)) )=2(m^(1/2))E/(ℏ (k^(1/2)))=2(((m/k))^(1/2))E/ℏ

Entonces la ecuación de Schroedinger toma la forma
(((d^2)ψ)/(d(y^2)))- (y^2)ψ=-ϵ ψ
(((d^2)ψ)/(d(y^2)))+(ϵ- (y^2))ψ= 0

Vamos a proponer como solución de la ecuación de Schroedinger el siguiente producto de funciones

phi[y_]:=f[y]*Exp[-y^2/2];

Con esta propuesta el lado izquierdo de la ecuación de Schroedinger tendrá la forma

phi''[y]+(epsilon-y^2)phi[y]

Factorizando a la exponencial

ecuacion1=Collect[%,{Exp[-y^2/2],f[y]}]

Como la ecuación de Schroedinger debe de ser igual a cero, retenemos la parte que multiplica a la exponencial

ecuacion2=ecuacion1/Exp[-y^2/2]

Vamos a proponer a la función f (y) como serie de potencias, ∑(Subscript[c, n])(y^n) con n=0,1,..N. En este caso usaremos 30 términos (N=30),

terminos=30;
f[y_]:=Sum[c[n]*y^n,{n,0,terminos}];

f[y]

La ecuacion2 tendrá la forma

ecuacion2

El coeficiente que multiplica la respectiva potencia de y se obtiene de

serie=CoefficientList[ecuacion2,y]//MatrixForm

De aquí obtenemos una relación de recurrencia

relacion[n_]:=Solve[serie[[1,n]]==0,c[n+1]]

con n >= 1

Table[relacion[i],{i,1,10}]//TableForm

Verifique en su cuaderno que las siguientes definiciones definiciones concuerdan con la relación de recurrencia anterior

c[n_]:=(2*n-3-epsilon)/((n-1)n)c[n-2]/;n>1
c[n_]:=c[1]/;n==1
c[n_]:=c[0]/;n==0

Es claro que las potencias pares de la variable y quedan en términos de c[0] y que las potencias impares quedan en términos del coeficiente c[1]. Por este motivo hagamos las siguientes definiciones

phipar[y_]:=Sum[c[n]*y^n,{n,0,terminos,2}]Exp[-y^2/2];
phiimpar[y_]:=Sum[c[n]*y^n,{n,1,terminos-1,2}]Exp[-y^2/2];

phipar[y]

De este resultado es claro que para algunos valores de epsilon la serie se trunca. Tomemos un valor particular de epsilon tal que algunos de los términos de la serie se cancele

valor=1;
Print["phi(y) = ",phipar[y]/.epsilon->valor];

El coeficiente de normalización tendrá el siguinete valor

normalizacionpar=Sqrt[Integrate[(phipar[y]*phipar[y]/c[0]^2/.epsilon->valor),{y,-Infinity,Infinity},Assumptions->Re[y^2]>0]]

funcionnormapar[y]=Simplify[(phipar[y]/.epsilon->valor)/(c[0]*normalizacionpar)]

potencialpar=Integrate[funcionnormapar[y]*y^2*funcionnormapar[y],{y,-Infinity,Infinity},Assumptions->Re[z^2]>0]

cineticapar=-Integrate[funcionnormapar[y]*D[funcionnormapar[y],{y,2}],{y,-Infinity,Infinity},Assumptions->Re[y^2]>0]

energiapar=potencialpar+cineticapar

La función tendrá la forma

Plot[{(funcionnormapar[y]/.y->variable)+energiapar,variable^2},{variable,-5,5},PlotRange->{0,3}]

Plot[{(funcionnormapar[y]*funcionnormapar[y]/.y->variable)+energiapar,variable^2,energiapar},{variable,-5,5},PlotRange->{0,3}]

Ejercicio. Encuentre los valores de la energía cinética y potencial para algunos valores de ϵ que trunquen la serie, y analice el comportamiento de las respectivas funciones de onda.

Los nodos en la función de onda se pueden encontrar a partir de ecuacion2 = 0

NSolve[(phipar[y]/.epsilon->5)==0,y]

El mismo procedimiento se aplica a las funciones impares

phiimpar[y]

valor=7;
Print["phi(y) = ",phiimpar[y]/.epsilon->valor];

normalizacionimpar=Sqrt[Integrate[(phiimpar[y]*phiimpar[y]/c[1]^2/.epsilon->valor),{y,-Infinity,Infinity},Assumptions->Re[y^2]>0]];

funcionnormaimpar[y]=Simplify[(phiimpar[y]/.epsilon->valor)/(c[1]*normalizacionimpar)]

potencialimpar=Integrate[funcionnormaimpar[y]*y^2*funcionnormaimpar[y],{y,-Infinity,Infinity},Assumptions->Re[y^2]>0]

cineticaimpar=-Integrate[funcionnormaimpar[y]*D[funcionnormaimpar[y],{y,2}],{y,-Infinity,Infinity},Assumptions->Re[y^2]>0]

energiaimpar=potencialimpar+cineticaimpar

Plot[{(funcionnormaimpar[y]/.y->variable)+energiaimpar,variable^2},{variable,-5,5},PlotRange->{0,8}]

Plot[{(funcionnormaimpar[y]*funcionnormaimpar[y]/.y->variable)+energiaimpar,variable^2,energiaimpar},{variable,-5,5},PlotRange->{0,8}]

NSolve[(phipar[y]/.epsilon->valor)==0,y]