Después de de una larga recopilación de datos dimos con la solución esta cuestión que tienen algunos los lectores. Te dejamos la respuesta y nuestro deseo es que sea de mucha ayuda.
Solución:
Solicitó enfoques alternativos a lo que hizo, así que aquí hay uno:
Un enfoque completamente diferente a la ecuación de Schrödinger unidimensional independiente del tiempo sería utilizar técnicas de matriz. La idea es eliminar la necesidad de NDSolve enteramente. Para problemas de estado ligado, puede hacer esto eligiendo una base que satisfaga la condición de la función de onda que desaparece en el infinito y expandiendo todas las cantidades en esa base.
Se ha señalado una forma particularmente agradable de hacer esto en:
HJ Korsch y M Glück, “Calcular valores propios cuánticos de forma fácil” Eur. J. Phys. 23 (2002) 413-426.
Utiliza la base del oscilador armónico, pero también introduce la relación entre los operadores de posición y momento por un lado, y los operadores del oscilador armónico por el otro. Esto lo hace muy útil como método de enseñanza.
Implementé ese enfoque aquí, con algunas mejoras para satisfacer las tolerancias dadas y hacer gráficos:
Se supone que la energía potencial $ U (x) $ es un polinomio de grado finito en una sola variable. Dado que estamos usando una base de oscilador armónico para encontrar los estados propios del hamiltoniano, un polinomio de grado finito se traduce en un conjunto de potencias finitas de la matriz de posición. xOp. Al aumentar la dimensión de estas matrices mientras se mantienen las potencias fijas en sus valores finitos, podemos lograr la convergencia de los autovalores y autovectores de energía cerca del mínimo potencial.
Cuanto mayor sea la dimensión de la matriz, más valores propios inferiores convergerán al true valores propios.
Primero defino los operadores de posición y momento en base a las funciones propias del oscilador armónico, y una matriz que representa la energía cinética. La dimensión de las matrices es dim:
Definir matrices básicas
xOp[dim_] :=
SparseArray[Band[1, 2] -> #, Band[2, 1] -> #, #, # &@dim] &[
Table[Sqrt[i], i, dim - 1]/Sqrt[2]]
eKin[dim_] :=
1/2 SparseArray[Band[1, 1] -> Table[i - 1/2, i, dim],
Band[1, 3] -> #, Band[3, 1] -> #, #, # &@dim] &[
Table[-Sqrt[i (i + 1)/2], i, dim - 2]/Sqrt[2]]
Espectro y funciones propias
Este es el corazón del cálculo.
La función llamada spectrum determina los autovalores y autovectores de un hamiltoniano con la energía cinética anterior y un término de energía potencial de tipo polinomial. El número de valores propios que se devolverán se da como primer argumento, n, y su error máximo está dado por tolerance (este último puede omitirse si el valor predeterminado es suficiente):
spectrum[n_, tolerance_: .00001][potential_, var_] /;
PolynomialQ[potential, var] := Module[
x, h, e, v,
c = CoefficientList[potential, var],
min = [email protected][
Minimize[potential, var],
Print[
"No absolute potential minimum found. Make sure the leading
power has a positive coefficient!"];
Abort[],
Minimize::natt
]
,
FixedPoint[
#[[1]] + 1,
Eigenvalues[
h = eKin[#[[1]]] +
N[Sum[c[[i]] MatrixPower[xOp[#[[1]]], i - 1], i, 2,
Length[c]] + (c[[1]] - min) IdentityMatrix[#[[
1]]]], -n] &, n + 1, 2 [email protected][n] - 1,
SameTest -> ([email protected][#1[[2]] - #2[[2]]] < tolerance &)];
e, v = Reverse /@ Eigensystem[h, -n];
e + min, ψe /@ v
]
El segundo conjunto de argumentos anterior (potential, var) son la energía potencial U(x) y el nombre de la variable de posición, normalmente sería x.
Tuve que hacer un cambio por compatibilidad con Mathematicaversión 10 porque arroja un error cuando intentas aplicar MatrixPower a una matriz numéricamente singular con potencia 0. Las matrices de posición pueden volverse numéricamente singulares para grandes dimensiones, por lo que tengo que tratar el poder cero de xOp por separado. Esto no era necesario en la versión 8. Dado que ese término solo determina un desplazamiento de energía constante, se agrupa junto con otro término constante que introduzco para cambiar el mínimo potencial a cero (de modo que el Eigenvalues se ordenan en orden ascendente, comenzando con el estado fundamental - este cambio por min se deshace al final del cálculo).
Resultados de salida en base a la posición
Para obtener la función de onda, implementamos la suma sobre los estados base para obtener las funciones de onda ψe (aquí, e no tiene nada que ver con la paridad, el método funciona independientemente de la simetría de paridad). El formato de salida se define a continuación para producir una representación abreviada que no escupe la lista completa de coeficientes, solo enumera su número.
Format[ψe[l_List]] := ψe[
Row["[LeftSkeleton]", Length[l], "[RightSkeleton]"]]
ψe[l_List][x_] := 1/(Pi^(1/4))Exp[-x^2/ 2] l.(1/Sqrt[2^# #!] HermiteH[#, x] &[Range[Length[l]] - 1])
plotSpectrum[data_, potential_, var_, min_, max_] :=
Show[
Apply[
Plot[
Evaluate[#1 + (Last[#1] - First[#1])/(2 + Length[#1])/
2 Through[#2[var]]], var, min, max,
Prolog -> Gray, Map[Line[min, #, max, #] &, #1],
Frame -> True
] &,
data
],
Plot[potential, var, min, max]
]
Ejemplos de
Compruebe que obtenemos los valores propios numéricos correctos para el potencial del oscilador armónico $ U (x) = frac 1 2 x ^ 2 $ en unidades adimensionales:
spectrum[10][1/2 x^2, x]
0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, ψe[<<15>>], el[<<15>>], el[<<15>>], el[<<15>>], el[<<15>>], el[<<15>>], el[<<15>>], el[<<15>>], el[<<15>>], el[<<15>>]
La primera lista contiene los valores propios, la segunda lista tiene las funciones propias en forma abreviada, lista para ser graficada. Observe de las formas cortas que el número de coeficientes usados en ψe es 15 aquí, aunque especifiqué que solo quiero encontrar 10 valores propios en el espectro. Esto se debe a que estos 10 Los valores propios no serán precisos para los valores dados. tolerance a menos que la matriz hamiltoniana esté diagonalizada en un espacio de Hilbert más grande de 15 funciones de base. Utilizando FixedPoint, el procedimiento incrementó iterativamente el número de funciones base hasta que ninguna de las primeras 10 los valores propios cambiaron en más de la tolerancia dada en la siguiente iteración.
Ejemplo 2: oscilador cuártico
Ahora veremos los estados propios de un potencial cuártico, que se muestra aquí:
Plot[-3 x^2 + x^4, x, -2, 2]
Este es un potencial de doble pozo, y los estados propios de tierras bajas estarán localizados excepto por un pequeño acoplamiento de túnel. La energía potencial sigue siendo una polinomio en $ x $, por lo que el método también debería funcionar aquí.
ψ1 = spectrum[5][-10 x^2 + x^4, x] // Last // Last
el[<<72>>]
Plot[Evaluate[ψ1[x]], x, -5, 5]
Si graficamos todos los estados propios juntos, las amplitudes se comprimen, pero puede discernir la distinción entre estados pares e impares:
With[v = -10 x^2 + x^4, n = 10,
plotSpectrum[
spectrum[n][v, x], v, x, -4, 4]
]
Como puede ver en el ejemplo del oscilador cuártico, el número de estados base (72) es más grande para potenciales más complicados. Pero para el caso de prueba del oscilador armónico, no es un problema obtener una precisión muy alta para niveles muy excitados:
First[spectrum[100][1/2 x^2, x]]
0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5 , 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5 , 50.5, 51.5, 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5, 59.5, 60.5, 61.5, 62.5, 63.5, 64.5, 65.5, 66.5, 67.5, 68.5, 69.5, 70.5, 71.5, 72.5, 73.5, 74.5 , 75.5, 76.5, 77.5, 78.5, 79.5, 80.5, 81.5, 82.5, 83.5, 84.5, 85.5, 86.5, 87.5, 88.5, 89.5, 90.5, 91.5, 92.5, 93.5, 94.5, 95.5, 96.5, 97.5, 98.5, 99.5
Como mencionó Jens, la discretización espacial de la ecuación es otra alternativa para el problema del estado ligado. Aquí está mi implementación muy simple de este enfoque. La idea básica es expresar la ecuación en una cuadrícula. Los diferenciales se pueden expresar como diferencias finitas. Por ejemplo, la derivada de segundo orden se puede expresar como
$$ frac d ^ 2 psi dx ^ 2 approx frac psi (x + Delta x) + psi (x- Delta x) -2 psi (x) Delta x ^ 2 $$
Entonces, el operador $ frac d ^ 2 dx ^ 2 $ se puede expresar como una matriz tridiagonal (condición de límite cero)
$$ left[ beginmatrix
-2 & 1 & & 0\
1 & ddots & ddots & \
& ddots & ddots & 1 \
0 & & 1 & -2 endmatrix right]
$$
y el operador potencial es una matriz diagonal
$$ left[ beginmatrix
V(x_1) & & & 0\
& ddots & & \
& & ddots & \
0 & & & V(x_N) endmatrix right]
$$
entonces la ecuación en forma de matriz es
$$ left[ beginmatrix
2a_0+V(x_1) & -a_0 & & 0\
-a_0 & ddots & ddots & \
& ddots & ddots & -a_0 \
0 & & -a_0 & 2a_0+V(x_N) endmatrix right]
cdot left[beginmatrix
Psi_1\
vdots\
vdots\
Psi_N
endmatrixright]
= E izquierda[beginmatrix
Psi_1\
vdots\
vdots\
Psi_N
endmatrixright]
$$
donde $ a_0 = frac hbar ^ 2 2m ( Delta x) ^ 2 $.
Diagonalizar el hamiltoniano dará los estados propios.
Código
TISE1D[U_Function, xmin_, xmax_, N0Grid_: 101,
BoundaryCondition_String: "zero"] := Module[
Δx = (xmax - xmin)/(N0Grid - 1), Hmtx, Tmtx, Vmtx,
Tmtx = -(1/(2 (Δx)^2))
SparseArray[i_, i_ -> -2, i_, j_ /; Abs[i - j] == 1 ->
1, N0Grid, N0Grid];
Vmtx = DiagonalMatrix[U /@ Range[xmin, xmax, Δx]];
Hmtx = Tmtx + Vmtx;
If[BoundaryCondition == "periodic",
Hmtx[[1, -1]] = Hmtx[[-1, 1]] = -(1/(2 (Δx)^2));
];
Sort[[email protected][Hmtx], (#1[[1]] < #2[[1]]) &]
]
Ejemplos de
1. oscilador armónico
V1[x_] = 1/2. x^2;
Transpose[TISE1D[Function[x, V1[x]], -10, 10, 400]][[1, 1 ;; 4]]
(*0.499921,1.49961,2.49898,3.49804*)
ListPlot[TISE1D[Function[x, V1[x]], -10, 10, 400][[1 ;; 4, 2]],
Joined -> True, PlotRange -> -5, 5, All, DataRange -> -10, 10,
Axes -> False, Frame -> True]
2.Potencial irreflexivo
V2[x_] = With[[HBar] = 1, m = 1, a = 1.4, -(([HBar]^2 a^2)/m) Sech[a x]^2];
Para este potencial, la función de onda de estado límite analítico es $ psi (x) = A * sech (ax) $ y la energía del estado fundamental es $ E_0 = - frac a ^ 2 hbar ^ 2 2 m $ .
grd = TISE1D[Function[x, V2[x]], -10, 10][[1, 2]];
With[a = 1.4, TISE1D[Function[x, V2[x]], -10, 10, 200][[1, 1]], -1.4^2/2]
(*-0.980757,-0.98*)
Show[ListPlot[Abs[grd]^2, PlotRange -> All, Joined -> True,
DataRange -> -10, 10],
Plot[Max[Abs[grd]^2]*Sech[1.4 x]^2, x, -10, 10, PlotRange -> All,
PlotStyle -> Red]]
3.potencial $ V (x) = - frac b hbar ^ 2 m text sech ^ 2 (ax) $.
V3[x_] = With[[HBar] = 1, m = 1, a0 = 10., a = 0.81, b = 1.94, -(([HBar]^2 b)/m) Sech[a x]^2];
bdE = Select[Transpose[TISE1D[Function[x, V3[x]], -10, 10, 200]][[1]]*27.211, # < 0 &]
(*-35.1018, -8.64135*)
ListPlot[TISE1D[Function[x, V2[x]], -10, 10, 200][[1 ;; 2, 2]],
PlotRange -> All, Joined -> True, DataRange -> -10, 10,
Axes -> False, Frame -> True]
Si el potencial es $ ( tanh (x) +1) ( tanh (x) -1) $, puede obtener la solución analítica usando Mathematica de la siguiente manera:
[I have omitted some of the detail - e.g. the asymptotic expansions - because the details are analogous to the simple harmonic oscillator case in my previous answer (see above).]
Defina el potencial.
u[x_] = (1 + Tanh[x]) (-1 + Tanh[x])
Defina la ecuación diferencial.
diffeq = -(1/2) [Psi]''[x] + u[x] [Psi][x] == e [Psi][x]
Resuelve la ecuación diferencial.
soln = DSolve[diffeq, [Psi], x][[1, 1]] // Simplify
(* [Psi] -> Function[x,
C[1] LegendreP[1, I Sqrt[2] Sqrt[e], Tanh[x]] +
C[2] LegendreQ[1, I Sqrt[2] Sqrt[e], Tanh[x]]] *)
Verifique esta solución.
diffeq /. soln // FullSimplify
(* True *)
[This is where I have omitted the asymptotic expansion details.] los LegendreQ pieza diverge como $ x longrightarrow infty $, así que use C[2] -> 0 para suprimirlo. los LegendreP pieza diverge como $ x longrightarrow - infty $ con una excepción cuando e -> -(1/2).
Entonces la solución de estado ligado es
[Psi][x] /. soln /. C[2] -> 0, e -> -(1/2)
(* 1/2 C[1] Sqrt[1 - Tanh[x]^2] *)
Aquí hay un Manipulate que ilustra cómo la solución de la ecuación diferencial depende de la energía.
Manipulate[
Plot[LegendreP[1, I Sqrt[2] Sqrt[e], Tanh[x]] // Chop, x, -5, 5,
AxesLabel -> "x", "[Psi]"], e, -1/2, "e", -1, 0]
Comentarios y calificaciones del post
Si crees que te ha sido útil nuestro artículo, sería de mucha ayuda si lo compartieras con el resto programadores y nos ayudes a dar difusión a este contenido.