Método variacional

Jorge Garza. Agosto de 2014

Estamos interesados en resolver la ecuación de Schrodinger asociada a

(-((ℏ^2)/(2m))(((d^2) )/(dx^2))+k (x^4))ψ=ϵψ

Al hacer el cambio de variable z = αx la ecuación es transformada en

-(((d^2)ψ )/(dz^2))+(z^4)ψ=λ ψ

o también

(((d^2)ψ )/(dz^2))+(λ-(z^4))ψ=0

Método variacional usando un parámetro

Comencemos con la función de prueba

psi[z_]=(1/Pi^(1/4))Exp[-z^2/2];

Comprobemos que la función está normalizada

Integrate[psi[z]^2,{z,-Infinity,Infinity}]

Construyamos el valor esperado

Integrate[-psi[z]D[psi[z],{z,2}]+psi[z] z^4 psi[z],{z,-Infinity,Infinity}]

%//N

Ahora, usemos una función de onda más flexible

norma=1/Sqrt[Integrate[Exp[-2α z^2],{z,-Infinity,Infinity},Assumptions->Re[α]>0]]

psi[z_]=norma*Exp[-α z^2]

Integrate[psi[z]^2,{z,-Infinity,Infinity},Assumptions->Re[α]>0]

energia=Integrate[-psi[z]D[psi[z],{z,2}]+psi[z] z^4 psi[z],{z,-Infinity,Infinity},Assumptions->Re[α]>0]

Obtengamos el valor de alfa que genere la energía más baja

solucion=Solve[D[energia,α]==0,α]

alfa=α/.solucion[[2]]//N

energia/.α->alfa

psi[y]/.α->alfa

Combinación lineal de funciones de base

Definamos al conjunto de funciones de base

f[i_,z_]:=(1/Sqrt[2^i Factorial[i] Sqrt[Pi]])HermiteH[i,z]Exp[-z^2/2]

Cambiemos el valor de n para ver cada una de las funciones

n=10;
f[n,z]

Verifiquemos si son funciones ortonormales

traslape[i_,j_]:=Integrate[f[i,z]f[j,z] ,{z,-Infinity,Infinity}]

n=5;
m=5;
traslape[n,m]

Ahora construyamos los elementos de matriz del operador hamiltoniano

h[i_,j_]:=Integrate[-f[i,z]D[f[j,z],{z,2}]+f[i,z] z^4 f[j,z],{z,-Infinity,Infinity}]

h[0,2]

Generemos la matriz H dependiendo del número de funciones de base

Clear[matrizH];
total=5;
matrizH=Table[N[h[i,j]],{i,0,total-1},{j,0,total-1}];

matrizH//MatrixForm

Con la matriz H construida se obtienen sus valores y vectores propios

valores=Eigenvalues[matrizH]

Es importante mencionar que el valor más bajo se encuentra en el último lugar de la variable valores.

vectores=Eigenvectors[matrizH];

funcion[k_,y_]:=Sum[vectores[[k,i]]*f[i-1,y],{i,1,total}];

estado=5;
Plot[{funcion[estado,y]+valores[[estado]],y^4,valores[[estado]]},{y,-5,5},PlotRange->{valores[[estado]]-1,valores[[estado]]+1}]

Partícula sometida a un potencial lineal

Ahora tratemos un potencial de la forma
V(x)=(Piecewise[{{∞, x≤0}, {x, 0<x<L}, {∞, x≥L}}])
donde m es una constante.
En este caso tratemos a las funciones solución de la partícula en una caja
(Subscript[ψ, n])(x)=(((2/L))^(1/2))sen(nπx/L) con n=1,2,..,∞
 para representar la función variacional. Usemos unidades tales que ℏ=1,m=1 y además L=1

fpc[n_,x_]:=Sqrt[2]Sin[n*Pi*x];

¿Son funciones ortonormales?

v[x_]:=x;

hpc[i_,j_]:=Integrate[-fpc[i,x]D[fpc[j,x],{x,2}],{x,0,1}]+Integrate[fpc[i,x] v[x] fpc[j,x],{x,0,1}]

Clear[matrizHpc];
total=4;
matrizHpc=Table[N[hpc[i,j]],{i,1,total},{j,1,total}];

solucionpc[k_,y_]:=Sum[vectorespc[[k,i]]*fpc[i,y],{i,1,total}];

valorespc=Eigenvalues[matrizHpc]

vectorespc=Eigenvectors[matrizHpc];

estado=4;
solucionpc[estado,x]
Plot[{solucionpc[estado,x]^2+valorespc[[estado]],v[x],valorespc[[estado]]},{x,0,1},PlotRange->{valorespc[[estado]]-2,valorespc[[estado]]+2}]