conocido como el potencial Saxon-Woods, utilizado en la física nuclear para describir el potencial que actúa entre un núcleo atómico y un neutrón que pasa cerca de dicho núcleo, en donde V0 representa la profundidad del pozo de potencial, a es una longitud que representa el “grosor de la superficie del núcleo”, y R0 = r0A1/3 es el radio del núcleo del átomo siendo r0 = 1.25 fermis y siendo A el número atómico (o número de masa) del átomo. No existe una solución analítica exacta para encontrar los niveles de energía de los nucleones dentro del átomo bajo un potencial Saxon-Woods, es necesario recurrir a una solución numérica.
El punto de partida para poder entender el método de solución numérica que estudiaremos aquí es la definición matemática formal de la derivada de primer orden:
Si en lugar de tomar el límite infinitesimal Δx.→.0 tomamos incrementos finitos, entonces tendremos la siguiente aproximación:
que se acercará más y más al valor exacto de la derivada conforme se disminuye el tamaño de Δx. Podemos despejar esta relación para obtener así:
Esto nos dá el valor de la función en un punto x+Δx cuando se tiene una expresión para la derivada dy/dx y cuando se tiene un punto de partida y(x). Con la finalidad de dar mayor claridad a lo que estamos haciendo, escribiremos lo anterior de la siguiente manera para dos puntos distintos pero cercanos x1 y x2:
El incremento que se está manejando aquí está definido de la siguiente manera:
Una notación alterna expresando lo mismo para el incremento anterior que resulta útil cuando se llevan a cabo evaluaciones numéricas a lo largo de varios puntos es la siguiente:
siendo n un número entero. De este modo, despejando para xn+1 tenemos:
Con esta notación de sub-índices, la expresión que habíamos obtenido arriba se puede escribir de la siguiente manera:
Obsérvese que se ha puesto también un sub-índice en la derivada de primer orden.
Lo anterior se basa en la definición de la derivada de primer orden conocida como la definición del incremento completo. Pero hay otra definición para la derivada de primer orden que para fines de cálculo numérico converge con mayor rapidez y precisión hacia la solución numérica buscada, la cual es completamente equivalente a la definición en la que se usa un incremento completo. Se trata de la definición del medio-incremento expresada simbólicamente de la siguiente manera:
Nuevamente, si en lugar de tomar el límite infinitesimal tomamos incrementos finitos aunque pequeños, lo anterior lo podemos escribir de la siguiente manera:
Con la finalidad de ilustrar mejor este concepto de derivada, a continuación se obtendrá la derivada aproximada de la función y(x) = x3 en el punto x1 = 2. Sabemos del cálculo infinitesimal que la derivada de esta función es dy/dx = 3x2, con lo cual la derivada de la función en el punto x1 = 2 será igual a 12. Supondremos que ignoramos cómo llevar a cabo el proceso de diferenciación. Si utilizamos un incremento Δx = 1, entonces la anterior fórmula nos producirá el siguiente resultado:
La respuesta obtenida no difiere mucho del valor exacto dy/dx = 12. Ahora utilizaremos un incremento Δx = 1/4:
Obsérvese que la respuesta se va acercando al valor exacto conforme el tamaño del incremento Δx se hace más pequeño, lo cual está de acuerdo con la definición formal de una derivada .
Ahora extenderemos lo que hemos visto arriba para poder llevar a cabo la solución numérica de ecuaciones diferenciales de segundo orden. Con este objetivo en mente, empezaremos con una ecuación diferencial de segundo orden extremadamente sencilla, la cual se obtiene aplicando la fórmula de Newton F.=.ma (la cual involucra en la variable aceleración a una derivada de segundo orden con respecto al tiempo) a la definición del peso gravitacional de una masa W = mg siendo g la aceleración ocasionada por la gravedad de la Tierra sobre un objeto en su superficie. Haciendo F = W y usando a la coordenada-y para medir el desplazamiento vertical de un cuerpo en su caída, se obtiene la siguiente ecuación diferencial:
Puesto que la derivada de segundo orden es realmente la derivada de una derivada de primer orden, lo anterior se puede escribir de la siguiente manera:
Ahora bien, usando incrementos finitos en lugar de infinitesimales, esto se puede considerar equivalente en una forma aproximada a una operación como la siguiente:
en donde n es un número entero positivo, con lo cual la ecuación diferencial de segundo orden se nos convierte en lo siguiente:
Despejando lo anterior obtenemos entonces:
Esta es nuestra primera ecuación. La otra ecuación que necesitamos nos viene directamente de lo que es la definición de una derivada de primer orden, excepto que usaremos incrementos finitos en lugar de infinitesimales, razón por la cual será una aproximación y no una igualdad, convirtiéndose en igualdad únicamente cuando el incremento Δt se hace tender a cero :
Despejando de esta última ecuación, obtenemos nuestra segunda ecuación:
Como puede verse, el método numérico que se acaba de describir es un método iterativo propio del análisis numérico, en el cual empezamos con un valor inicial y0 conocido como la semilla para la cual utilizamos n = 0, obteniendo un n+1 valor nuevo en cada paso usando el valor n obtenido previamente. Esto nos permite establecer un procedimiento para ir obteniendo los pares de puntos de una expresión matemática con los cuales podemos obtener una solución gráfica de la misma, el cual consiste en los siguientes pasos:
(1) Tras seleccionar un incremento Δx razonablemente pequeño, y partiendo de los valores iniciales y0 y dy/dx0 (por ejemplo, y0 = 1 y dy/dx0 = 0), utilizamos la relación:
para obtener el siguiente valor y1.
(2) A continuación, dy/dx0 es utilizado para evaluar dy/dx1 usando la definición de la derivada con
(3) Tomando los valores de y1 y dy/dx1 obtenidos en los pasos anteriores, se calcula y2:
(4) A continuación, dy/dx1 es utilizado para evaluar dy/dx2:
(5) Repitiendo los pasos anteriores, se utilizan y2 y dy/dx2 para evaluar y3 y dy/dx3 .
(6) El procedimiento iterativo se repite cuantas veces sea necesario para obtener la cantidad de pares de puntos (x,y) que se consideren necesarios para la construcción de una gráfica adecuada.
(7) Si la cantidad C es una función de la variable independiente, entonces será necesario ir evaluando dicha cantidad en cada punto nuevo de la variable independiente que se vaya utilizando para evaluar los yk y los dy/dxk.
El procedimiento iterativo que se acaba de delinear está basado en la definición de la derivada que recurre a un incremento completo. Modificándolo para basar el procedimiento sobre la definición del medio incremento, los pasos a seguir para llevar a cabo la solución numérica de una ecuación diferencial de segundo orden del tipo d²y/dx² son:
En este caso, la solución numérica (y gráfica) de la ecuación de onda mediante el método del medio-incremento se lleva a cabo en forma iterativa de la siguiente manera:
Los pasos a seguir cuando se emplea el método del medio-incremento son entonces los siguientes:
(1) Tras seleccionar un incremento Δx razonablemente pequeño, y partiendo del valor inicial dy/dx0, se obtiene el valor aproximado de dy/dx1/2, o sea la derivada a una mitad del incremento Δx posterior al valor inicial x0, usando la relación:
(2) Se obtiene el valor de y1 a partir de los valores de dy/dx1/2 y de y0 usando la relación:
(3) Lo anterior completa la primera parte del bucle iterativo. Ahora se obtiene el valor aproximado de dy/dx3/2 usando la relación:
(4) Se obtiene el valor de y2 a partir de los valores de dy/dx3/2 y de y1 usando la relación:
(5) Repitiendo los pasos anteriores, se utilizan y2 y dy/dx3/2 para evaluar y3 y dy/dx5/2.
(6) El procedimiento iterativo se repite cuantas veces sea necesario para obtener la cantidad de pares de puntos (x,y) que se consideren necesarios para la construcción de una gráfica adecuada.
(7) Si la cantidad C es una función de la variable independiente, entonces será necesario ir evaluando dicha cantidad en cada punto nuevo de la variable independiente que se vaya utilizando para evaluar los yk y los dy/dxk.
Para una gran mayoría de casos prácticos que puedan ser considerados, el método del medio-incremento resulta ser mucho más eficiente y mucho más rápido y exacto que el método basado en la definición de la derivada que se apoya en el incremento completo. Naturalmente, posiblemente a estas alturas surge la pregunta en el lector del por qué ésto es así. La respuesta la podemos visualizar mejor comparando la gráfica de la definición clásica de la derivada basada en el incremento completo:
con la definición de la derivada basada en el medio-incremento:
La derivada de una función f(x) en un punto P de la curva trazada por dicha función no es más que la pendiente de la tangente trazada a la curva en dicho punto, y la definición usual de la derivada basada en un incremento Δx se obtiene haciendo Δx.→.0, y conforme esto ocurre la pendiente de la recta trazada entre dos puntos distintos de la curva (la línea de color negro en el primer dibujo) se va acercando más y más a la pendiente de la tangente a la curva en el punto P (la línea de color rojo). Si contrastamos la acción de la derivada basada en el medio-incremento, podemos comprobar en el segundo dibujo que la pendiente entre los dos puntos de la curva que definan al medio-incremento Δx/2 se va acerca mucho más rápidamente que la pendiente entre dos puntos similares cuando se utiliza la definición de derivada basada en el incremento completo Δx (obsérvese que la línea negra entre los dos puntos de la curva que definen al incremento está más “paralela” a la línea roja en el segundo dibujo que en el primero).
Una vez delineado el procedimiento general, lo utilizaremos para el caso que ocupa nuestro interés principal, la ecuación de onda de Schrödinger, cambiando ligeramente la notación a la notación que estaremos utilizando en el procedimiento iterativo para resolver numéricamente la ecuación diferencial de segundo orden expresada en forma compacta como:
teniendo así el siguiente conjunto de ecuaciones:
El esquema iterativo es idéntico al esquema general requiriendo un simple cambio de notación, escribiendo ψ en lugar de la variable que corresponde a la coordenada-x, llevándose a cabo de la siguiente manera:
A continuación usaremos el método numérico para resolver una ecuación de onda en la cual el potencial no es constante sino que está especificado mediante la siguiente relación:
Este potencial corresponde clásicamente al de una masa conectada a un resorte. Con este potencial, la ecuación de Schrödinger:
toma el siguiente aspecto:
Esta es la ecuación diferencial que según la Mecánica Cuántica describe a nivel sub-microscópico el mismo sistema que de acuerdo a la mecánica clásica está descrito por la relación d²x/dt² = -k/m. ¡Qué diferencia! Antes de aplicar el método del medio-incremento para resolver numéricamente esta ecuación de onda, usaremos la siguiente relación que nos conecta a la constante k del resorte con la frecuencia de oscilación:
Con esta substitución, la ecuación de onda toma la forma:
Podemos escribir lo anterior como:
En esta ecuación diferencial, al igual que en muchas otras, se vuelve deseable simplificarlas un poco ocultando las constantes físicas de modo tal que podamos enfocarnos sobre la cuestión puramente matemática. Esto lo podemos lograr aquí haciendo:
De este modo, tenemos la siguiente ecuación diferencial simplificada:
Podemos hacer esto más matemáticamente “puro” escribiéndolo en términos de una variable adimensional que sea un número sin dimensiones físicas:
con lo cual se tiene:
Entonces la ecuación de onda toma la forma:
Podemos simplificar esto aún más definiendo el primer término dentro del paréntesis como una cantidad ε que resulta ser:
Entonces la ecuación diferencial a resolver numéricamente:
Definiendo en el lado derecho de la igualdad la cantidad C de la siguiente manera:
tenemos entonces la siguiente ecuación diferencial a resolver numéricamente:
La prescripción (se vuelve tentador usar la frase “receta de cocina”) para la solución numérica de la ecuación diferencial, haciendo los cambios correspondientes en la notación, es entonces:
Por simplicidad, en el esquema iterativo dado arriba no aparece mostrada la evaluación en cada paso del valor que va tomando la cantidad C, lo cual es necesario en virtud de que es una función de u y de ψ(u). Si empezamos con el siguiente conjunto de valores que corresponden al comienzo en u0.=.0 de una solución par del tipo ψ(u) = Acos(ku):
ψ0 = 1.0___u0 = 0___dψ/du0 = 0
Δu = 0.05___ε = 1.0
Δu = 0.05___ε = 1.0
tenemos entonces la forma en la cual se van generando los pares de puntos para poder llevar a cabo la gráfica:
(u0, ψ0) = (0, 1)
C = (u0² - ε)ψ0 = (0² - 1)(1) = -1
dψ/du1/2 ≈ dψ/du0 + C(Δu/2)
≈ 0 + (-1)(0.05/2)
≈ - 0.025
≈ - 0.025
ψ1 = ψ0 + dψ/du1/2·Δu
= 1 - (- 0.025)(0.05)
= 0.99875
= 1 - (- 0.025)(0.05)
= 0.99875
u1 = u0 + Δu
= 0 + 0.05
= 0.05
= 0 + 0.05
= 0.05
(u1, ψ1) = (0.05, 0.99875)
C = (u1² - ε)ψ1 = (0.05² - 1)(0.99875) = - 0.99625
dψ/du3/2 ≈ dψ/du1/2 + CΔu
≈ - 0.025 + (- 0.99625)(0.05)
≈ - 0.07481
≈ - 0.07481
ψ2 = ψ1 + dψ/du3/2·Δu
= 0.99875 + (- 0.07481)(0.05)
= 0.995
= 0.99875 + (- 0.07481)(0.05)
= 0.995
u2 = u1 + Δu
= 0.05 + 0.05
= 0.10
= 0.05 + 0.05
= 0.10
(u2, ψ2) = (0.10, 0.995)
C = (u2² - ε)ψ2 = (0.10² - 1)(0.995) = - 0.98505
dψ/du5/2 ≈ dψ/du3/2 + CΔu
≈ - 0.07481 + (- 0.98505)(0.05)
≈ - 0.12406
≈ - 0.12406
ψ3 = ψ2 + dψ/du5/2·Δu
= 0.995 + (- 0.12406)(0.05)
= 0.9888
= 0.995 + (- 0.12406)(0.05)
= 0.9888
u3 = u2 + Δu
= 0.10 + 0.05
= 0.15
= 0.10 + 0.05
= 0.15
(u3, ψ3) = (0.15, 0.9888)
C = (u3² - ε)ψ3 = (0.15² - 1)(0.9888) = - 0.96655
dψ/du7/2 ≈ dψ/du5/2 + CΔu
≈ - 0.12406 + (- 0.96655)(0.05)
≈ - 0.172388
≈ - 0.172388
ψ4 = ψ3 + dψ/du7/2·Δu
= 0.9888 + (- 0.172388)(0.05)
= 0.9802
= 0.9888 + (- 0.172388)(0.05)
= 0.9802
u4 = u3 + Δu
= 0.15 + 0.05
= 0.20
= 0.15 + 0.05
= 0.20
(u4, ψ4) = (0.20, 0.9802)
Los demás pares sucesivos de puntos que se quieran obtener se pueden ir calculando de la misma manera. Continuando con la evaluación numérica de pares de coordenadas adicionales, tendremos entonces que la gráfica que corresponde a los distintos valores numéricos de ψ(u) en función de u obtenidos con las condiciones iniciales dadas arriba viene siendo la siguiente (esta gráfica y las que le siguen se pueden apreciar un poco mejor ampliando cada una de las imágenes):
Quienes tengan acceso a una calculadora científica programable de bolsillo se darán cuenta de que el procedimiento iterativo que se acaba de delinear se puede programar sin mayores problemas. Y quienes tienen acceso a un programa computacional de graficados pueden obtener varias gráficas para distintos valores del parámetro ε que es a fin de cuentas la verdadera variable. Como puede apreciarse en la gráfica, para el valor ε.=.1.00 tanto la función ψ(u) como la derivada de la función dψ(u)/du se aproximan a cero conforme u.→.+∞. Estos son precisamente los dos criterios esenciales que cuando son cumplidos ambos sabemos que hemos encontrado una eigensolución admisible para la ecuación diferencial de onda. En pocas palabras, el valor ε.=.1.00 nos dá la primera eigensolución (numérica) a la ecuación diferencial:
¿Y cómo podemos estar seguros de que la solución que hemos encontrado es efectivamente una solución válida? Pues obteniendo otros conjuntos de pares de puntos para valores de ε tanto mayores que 1.00 como menores que 1.00. Si probamos con una selección ε = 1.2, o sea si resolvemos el problema con el siguiente conjunto de valores iniciales:
ψ0 = 1.0___u0 = 0___dψ/du0 = 0
Δu = 0.05___ε = 1.2
Δu = 0.05___ε = 1.2
obtendremos entonces una nueva colección de pares de puntos que graficados nos producirán el siguiente resultado:
Con la condición ε.=.1.2 hemos encontrado que la función ψ(u) no sólo nunca se hace cero, sino que en cierto momento la curva cambia de concavidad volviéndose cóncava hacia arriba antes del punto en el cual se debería de haber llevado a cabo el cambio de concavidad para hacer que la función ψ(u) fuera “aterrizando” en forma suave y gradual justo en la cantidad requerida para tocar el eje horizontal en la forma en la que lo hace una solución admisible, cumpliendo con los dos criterios básicos. Para que la solución sea admisible, el cambio en la concavidad se debe llevar a cabo justo en el punto en donde se debe llevar a cabo para que se puedan cumplir los dos criterios de solubilidad, ni antes ni después. En este caso, hemos encontrado que para un valor ligeramente mayor al valor ε.=.1.0 la función ψ(u) eventualmente se volverá infinitamente grande en el sentido positivo para u.→.+∞.
Hagamos ahora otro experimento numérico. En este caso, obtendremos otro conjunto de pares de puntos para un valor de ε menor que 1.00. Si probamos con una selección ε = 0.98, o sea si resolvemos el problema con el siguiente conjunto de valores iniciales:
ψ0 = 1.0___u0 = 0___dψ/du0 = 0
Δu = 0.05___ε = 0.98
Δu = 0.05___ε = 0.98
obtendremos entonces una nueva colección de pares de puntos que graficados nos producirán el siguiente resultado:
Con la condición ε.=.0.98 encontramos ahora que la función ψ(u) efectivamente sí se hace igual a cero, aproximadamente en u = 2.35 pero la pendiente de la curva dψ(u)/du en ese punto no es igual a cero lo cual significa que la curva atraviesa el eje horizontal iniciando su caída hacia el infinito negativo. De este modo, para un valor menor al valor de ε.=.1.0 la función ψ(u) eventualmente se volverá infinitamente grande en el sentido negativo para u.→.+∞.
Si exploramos otros valores inferiores a ε.=.1.00, no encontraremos algo que nos produzca un comportamiento como el que hemos obtenido arriba en el cual tanto la función ψ(u) como la derivada de la función dψ(u)/du se vayan aproximando a cero conforme u.→.+∞; esto es, no hay otros eigenvalores de la ecuación diferencial por debajo de ε.=.1.00. Esto nos deja como única alternativa explorar valores mayores que ε.=.1.00. Al explorar valores cada vez mayores de ε por la vía numérica por encima del valor ε.=.1.00, eventualmente encontraremos que la solución ε.=.1.00 no es la única, y aunque en el valor ε.=.2.00 no hay solución alguna eventualmente encontraremos que hay otra solución “crítica” de ε para ε.=.3.00. Esta solución la podemos obtener utilizando un conjunto de valores iniciales como el siguiente (obsérvese que para poder obtener el valor inicial de dψ/du este ya no puede ser tomado como igual a cero sino como dψ/du0.≈.1.65):
ψ0 = 0.0___u0 = 0___dψ/du0 = 1.65
Δu = 0.02___ε = 3.00
Δu = 0.02___ε = 3.00
con lo cual obtendremos entonces una nueva colección de pares de puntos que graficados nos producirán el siguiente resultado:
Nuevamente, tanto la función ψ(u) como la derivada de la función dψ(u)/du se aproximan a cero suavemente al mismo tiempo conforme u.→.+∞, con lo cual sabemos que hemos encontrado otra solución a la ecuación diferencial en ε = 3.00. A diferencia de la solución numérica encontrada en el caso ε.=.1.00 en donde el valor inicial de la pendiente dψ(u)/du era igual a cero en el origen, confirmado con la gráfica que obtuvimos, en este caso para la segunda solución la pendiente dψ(u)/du no es igual a cero, sino que tiene cierto valor que inicialmente tenemos que determinar al tanteo. ¿Y cómo podemos saber si el valor inicial de la pendiente dψ(u)/du se puede tomar como cero o se tiene que tomar como diferente de cero en el origen? Si tuviéramos a la mano la expresión analítica para la función de onda ψ(u) de la ecuación diferencial que estamos resolviendo numéricamente, la respuesta a la pregunta sería sencilla, ya que todo lo que tendríamos que hacer sería tomar la derivada de la función de onda y evaluarla en el punto u.=.0. Por ejemplo, para la función de onda:
ψ(u) = Asen(ku)
la derivada es:
dψ(u)/du = kAcos(ku)
y en el punto de origen u = 0 la pendiente tendrá el valor:
dψ(0)/du = kAcos(0) = kA
El problema es que cuando decidimos resolver una ecuación diferencial por la vía numérica no tenemos una expresión para la función de onda, siendo esta precisamente la razón por la cual recurrimos a la vía de la solución numérica de la ecuación diferencial, cuando hemos llegado a la conclusión de que no existe una solución analítica a la ecuación diferencial que estamos tratando de resolver. Esta situación nos advierte de antemano que la solución numérica de las ecuaciones diferenciales es tanto un arte como una ciencia, y además de las “pistas” que el investigador pueda aprender en los libros del análisis numérico eventualmente irá desarrollando algunas mañas que le permitirán tomar atajos.
Encima del valor ε.=.3.00 eventualmente encontraremos que hay otro valor “crítico” de ε para el cual tanto la función de onda ψ(u) como la pendiente de la función de onda dψ(u)/du se irán acercando simultáneamente a cero justo en la combinación requerida conforme u.→.+∞. Este valor es ε.=.5.00. Y encima de este valor “crítico” hay otro valor “crítico” que es ε.=.7.00. Juntando estos valores “críticos” nos damos cuenta de que hay un patrón matemático bien definido para las eigensoluciones de la energía:
La expresión general que nos genera este conjunto de valores de energías discretizadas es:
Esta solución que hemos obtenido, por la vía numérica, es de hecho la expresión general para los eigenvalores de la energía del oscilador armónico simple.
Otro de los valores “críticos” de ε para los cuales tanto la función de onda ψ(u) como la pendiente de la función de onda dψ(u)/du se irán acercando simultáneamente a cero justo en la combinación requerida conforme u.→.+∞. es ε.=.17.0, el cual equivale a resolver el problema numéricamente con el siguiente conjunto de valores iniciales:
ψ0 = 1.0___u0 = 0___dψ/du0 = 0
Δu = 0.05___ε = 17.0
Δu = 0.05___ε = 17.0
La gráfica obtenida bajo estas condiciones es la siguiente:
Se sabe de antemano por consideraciones teóricas que el valor de ε.=.17 para la ecuación diferencial sobre la cual hemos estado trabajando dá una solución admisible exacta a la ecuación de onda, la cual corresponde de hecho al noveno eigenestado. Y aunque ciertamente vemos aquí el comportamiento oscilatorio esperado que corresponde a este eigenestado, la función de onda ψ(u) en vez de irse acercando asintóticamente al eje horizontal lo cruza y empieza a desplomarse hacia el infinito en el sentido negativo conforme u.→.+∞. ¿Qué es lo que está sucediendo aquí? Desafortunadamente, el efecto discrepante se debe al error de redondeo numérico que eventualmente nos alcanza hasta interferir con nuestro análisis. Para compensar por el error de redondeo, una cosa que podemos hacer es disminuír el incremento Δu de 0.05 a 0.01 ó mejor aún, a 0.001. Pero esto nos puede terminar elevando el tiempo de procesamiento de modo tal que algo que antes nos llevaba un minuto en procesar con una computadora de escritorio nos puede terminar costando 100 minutos, varias horas, o inclusive varios días. Eventualmente tenemos que llegar a un compromiso sobre la precisión que deseamos obtener. Una buena estrategia consiste en empezar con unos valores “toscos” empleando incrementos razonablemente grandes de Δu, y cuando creamos estar cerca de una solución podemos ir disminuyendo el incremento Δu según se vaya requiriendo para ir aumentando gradualmente nuestra precisión sin incurrir en pérdidas exageradas de tiempo de procesamiento. Aquí la experiencia que vayamos acumulando será nuestra mejor guía. Desafortunadamente, otro obstáculo insalvable con el que tenemos que lidiar es la precisión numérica limitada de las calculadoras programables de bolsillo, y eventualmente, con la precisión numérica limitada hasta de las computadoras más poderosas si es que tenemos acceso a ellas. Si tuviésemos a nuestra disposición una computadora capaz de llevar a cabo cálculos numéricos con una precisión infinitamente grande, sin errores de redondeo, siempre podríamos obtener soluciones numéricas exactas a todo tipo de ecuaciones diferenciales que no muestren jamás divergencia alguna en las gráficas. Pero tenemos que reconocer nuestras limitaciones tecnológicas, y hacer lo mejor que podamos para poder salir adelante de alguna manera. Como ya se mencionó, la solución numérica de ecuaciones diferenciales es tanto un arte como una ciencia, y se requiere de cierta experiencia y de cierto “colmillo” para saber cómo movernos en esas aguas. Si hay alguna consolación en todo esto es que el procedimiento numérico de solución es directo y no presenta mayores problemas una vez que se ha entendido su mecánica. Más aún, la disponibilidad comercial de paquetes computacionales como Mathematica y MathCAD permite la generación y el graficado de lo anterior de manera automática, permitiendo hacer en cuestión de minutos un trabajo que antes del advenimiento de las computadoras de escritorio se podía llevar varias semanas o inclusive meses de duro trabajo.
En los procedimientos de solución numérica empleados en el ejemplo de arriba, vimos cómo las condiciones ψ0 = 1.0 y dψ/dx0 = 0 nos generan una gráfica que será simétrica con respecto al eje vertical que pasa por el origen (por el punto x = 0), una simetría como la de un espejo en virtud de la cual si llevamos a cabo la solución numérica para valores negativos de x comprobaremos que para la función de onda que corresponde a este tipo de solución se debe cumplir la condición ψ(-x) = ψ(x). Este tipo de función de onda es lo que se conoce como una función par (llamada también gerade que significa lo mismo traducido del alemán). Por otro lado, las condiciones ψ0 = 0 y dψ/dx0 ≠ 0 nos generan una gráfica que también será en cierto modo simétrica con respecto al eje vertical que pasa por el origen, pero esta no es una simetría como la de un espejo ya que todo el lado de la gráfica que vá del lado izquierdo (correspondiendo a valores negativos de x) se invierte verticalmente, por lo cual para la función de onda que corresponde a este tipo de solución se debe cumplir la condición ψ(-x) = -ψ(x). A este tipo de función se le conoce como una función impar (llamada también ungerade que significa lo mismo traducido del alemán). Como regla general, los estados con paridad par ψ(-x) = ψ(x) no se desvanecen en el centro del potencial V(x) pero en cambio tienen una derivada de primer orden dψ/dx que es igual a cero, mientras que los estados con paridad impar ψ(-x) = -ψ(x) se vuelven cero en el centro del potencial pero no tienen una derivada de primer orden dψ/dx que se desvanezca en el centro del potencial V(x).
Habiendo visto cómo funcionan los métodos de solución numérica para la ecuación de Schrödinger en el caso del oscilador armónico simple, tenemos a continuación otro caso unidimensional en el cual el potencial V(x) representa un potencial rectangular de pozo finito que podemos dibujar de la siguiente manera:
Como puede verse, el pozo de potencial tiene una profundidad V0, y en el fondo del pozo al potencial se le asigna un valor igual a cero. El pozo en sí tiene una anchura a, y si situamos el origen de las coordenadas justo a la mitad de la distancia que hay entre las paredes del pozo entonces cada pared estará situada en el punto x = -a/2 y en el punto x = +a/2. Siguiendo el modo de pensar de la mecánica clásica, se considera que una partícula que está dentro del pozo se encuentra en cierto modo “atrapada” o ligada, razón por la cual forma parte de un estado ligado, mientras que una partícula que se encuentre fuera del pozo con una energía potencial mayor que V0 será esencialmente una partícula libre, razón por la cual forma parte de un estado no-ligado. Puesto que dentro del pozo la energía potencial es igual a cero, la ecuación de Schrödinger:
toma el siguiente aspecto dentro del pozo:
Nuevamente, resulta conveniente convertir esta ecuación física en una ecuación adimensional como lo hicimos anteriormente mediante un cambio de variable:
Con este cambio de variable, la ecuación diferencial que tenemos que resolver numéricamente será;
en donde:
De nueva cuenta, en la solución numérica que será llevada a cabo, nos limitaremos a la obtención de pares de puntos para los cuales u tendrá un valor positivo, en el entendido que por la simetría de la curva obtenida será relativamente fácil invertir la gráfica para tener la solución tratándose de valores negativos de u. Debe destacarse el hecho de que aquí estamos definiendo a ε como el cociente entre la energía de la partícula y la altura V0 del pozo de potencial, de modo tal que para una partícula atrapada dentro del pozo estaremos probando varios valores de ε mayores que cero pero menores a la unidad. Para las dimensiones atómicas usuales encontradas en el mundo sub-microscópico, un valor típico de β es el de β = 64. La ecuación diferencial que tenemos entre manos será válida dentro del pozo, esto es, para valores positivos de u menores que u = +1/2. En el caso de valores positivos de u mayores que u = +1/2, la ecuación diferencial tomará otro aspecto y tiene que ser resuelta por separado. Con esto en mente, las condiciones iniciales para la búsqueda de los eigenestados del pozo de potencial serán:
ψ0 = 1.0___u0 = 0___dψ/du0 = 0
Δu = 0.02___β = 64
Δu = 0.02___β = 64
En este caso, el parámetro que estaremos variando para la búsqueda de valores “críticos”será ε. Tras una búsqueda que parte de cero hacia arriba (en el sentido positivo), eventualmente confirmamos que un eigenestado se encontrará aproximadamente en:
ε = 0.0987
La gráfica de ψ(u) que corresponde a los puntos generados a través del análsis numérico sobre este problema para las condiciones iniciales dadas es la siguiente (la línea roja representa la localización de la pared del pozo en el sentido positivo de la coordenada-u):
¿Y cómo podemos estar tan seguros de que el valor ε = 0.0987 es un eigenvalor que representa un eigenestado del sistema para una partícula dentro del pozo? Si probamos con otro valor inferior, digamos ε = 0.094, la curva no sólo no cruzará el eje horizontal en ψ(u) = 0, ni siquiera lo tocará ya que antes de que ello ocurra la curva empezará a girar hacia arriba manifestando una divergencia en el sentido positivo. Obsérvese que en este eigenestado algo de la función de onda de la partícula atraviesa la pared derecha del pozo (en el lado contrario, para valores negativos de u, ocurre exactamente lo mismo). Clásicamente, esto no puede ocurrir, ya que la partícula sólida no puede “atravesar” las paredes del pozo finito, ni siquiera parcialmente. Pero en la Mecánica Cuántica la partícula sólida ya no es tan sólida, se puede comportar también como una onda de materia, y debemos estar preparados para este tipo de sorpresas.
Como es de esperarse, un pozo de potencial de profundidad finita no admite un número infinitamente grande de estados discretizados. Para que el pozo de potencial pueda acomodar una cantidad infinitamente grande de estados eigen, se requiere que el pozo sea un pozo de profundidad infinita, o lo que es lo mismo, que posea paredes infinitamente altas que actúen como barrera impenetrable. De cualquier modo, para un pozo de potencial finito como el que aquí tenemos es importante continuar con la búsqueda de otros estados eigen que pueda haber por allí esperando ser descubiertos. Variando las condiciones iniciales, eventualmente encontramos un segundo eigenestado, el cual está ubicado aproximadamente en:
ε = 0.383
bajo las siguientes condiciones iniciales:
ψ0 = 0.0___u0 = 0___dψ/du0 = 5.0
Δu = 0.02___β = 64
Δu = 0.02___β = 64
La gráfica de ψ(u) que corresponde a los puntos generados a través del análsis numérico sobre este problema para las condiciones iniciales dadas es la siguiente (nuevamente, la línea roja representa la localización de la pared del pozo en el sentido positivo de la coordenada-u):
El procedimiento señalado nos permite obtener numéricamente los tres eigenestados de una partícula situada dentro de un pozo de potencial en una región en la cual el potencial V(x) es igual a cero. Pero no nos dice nada acerca de la región situada “fuera” de las paredes del pozo en la cual el potencial es diferente de cero, en donde la partícula se encuentra en un estado no-ligado (unbounded). Arriba de esta región con respecto al nivel en el cual el potencial es V0, tenemos que echar mano nuevamente de la ecuación de Schrödinger que para este caso será, después de una simplificación para convertirla en una ecuación diferencial adimensional:
La solución numérica que intentaremos llevar a cabo sobre esta ecuación diferencial para una partícula situada fuera del pozo de potencial será efectuada con las siguientes condiciones iniciales:
ψ0 = 1.0___u0 = 0___dψ/du0 = 0
Δu = 0.02___β = 64
Δu = 0.02___β = 64
empleando en un comienzo un valor de ε igual a 1.4 (obsérvese que a diferencia de los casos anteriores para una partícula situada dentro del pozo de potencial, aquí ε tomará valores mayores que la unidad).
Obsérvese con atención cómo la amplitud de la función de onda de la partícula es mayor fuera de la región que corresponde al exterior de pozo que la amplitud de la misma dentro de la región que corresponde al interior del pozo, pese a que la partícula está fuera del pozo. Es como si el pozo de alguna manera influyera en lo que clásicamente vendría siendo una partícula libre. Esto lo podemos apreciar mejor trazando la gráfica completa para incluír tanto los valores negativos como los valores positivos de u, poniendo debajo de la misma el pozo:
Es importante darse cuenta de que, a diferencia de lo que ocurre con los estados ligados (en el interior del pozo), no hay dificultad alguna en obtener una solución a la ecuación de Schrödinger para un parámetro de energía ε mayor que la unidad (correspondiendo a una energía mayor que E/V0), ya que la solución permanece oscilatoria sin importar cuán grande sea el valor de u, y sin importar tampoco cuál sea el valor de E; no hay tendencia alguna en la función de onda ψ a mostrar una divergencia hacia el infinito cuando la energía E de la partícula está por encima del borde del pozo. Puesto que clásicamente la partícula no está ligada al interior del pozo, se puede concluír que cuando la relación entre la energía potencial V de una partícula y su energía total E es tal que la partícula no está ligada a una región limitada del espacio, entonces en la Mecánica Cuántica todos los valores de E son permitidos. Ya no hay eigenestados, y los valores permitidos de E forman lo que llamamos un continuo o continuum. Todos los valores de E ó ε en el continuum arriba de V0 o de la unidad son iguales en el sentido de que todos ellos están permitidos. Sin embargo, algunos de los valores de E en el continuum definitivamente son más iguales que otros, y esto lo descubriremos probando el siguiente conjunto de condiciones iniciales:
ψ0 = 1.0___u0 = 0___dψ/du0 = 0
Δu = 0.02___β = 64
Δu = 0.02___β = 64
empleando en el arranque el análisis numérico un valor de ε igual a 2.5. La gráfica de ψ(u) que corresponde a los puntos generados a través del análsis numérico sobre este problema para las condiciones iniciales dadas es la siguiente:
Obsérvese que en este caso sucede algo muy curioso. La amplitud de la función de onda de la partícula es la misma en la región que corresponde al exterior del pozo que la amplitud en la región que corresponde al interior del pozo. Nuevamente, esto lo podemos apreciar mejor trazando la gráfica completa para incluír tanto los valores negativos como los valores positivos de u, poniendo debajo de la misma el pozo:
Esto que acabamos de descubrir ocurre no sólo para el valor de ε igual a 2.5. Podemos encontrar una infinidad de valores (discretos) de ε para los cuales la amplitud de la función de onda será la misma tanto en la región que corresponde al exterior del pozo como en la región que corresponde al interior del pozo, como si la presencia del pozo hubiera dejado de tener efecto alguno para estos valores específicos de ε. Estos estados, ubicados para una partícula que está fuera del pozo, son conocidos como estados virtuales. Son una especie de “eigenestados” sin serlo realmente. Inclusive para el valor de ε igual a 1.4 que ya habíamos probado antes, tenemos también un estado virtual, si utilizamos las condiciones iniciales correctas, que en este caso implica utilizar una función de onda impar para la cual ψ(-u) = -ψ(u) en lugar de una función de onda par para la cual ψ(-u) = ψ(u), o sea las condiciones iniciales:
ψ0 = 0___u0 = 0___dψ/du0 = 9.45
Δu = 0.02___β = 64
Δu = 0.02___β = 64
empleando en el arranque el análisis numérico un valor de ε igual a 1.4. La gráfica de ψ(u) que corresponde a los puntos generados a través del análsis numérico sobre este problema para las condiciones iniciales dadas es la siguiente:
Nuevamente, podemos ver con mayor claridad cómo la amplitud de la función de onda se mantiene constante tanto fuera como dentro de las regiones que corresponden al pozo de potencial incluyendo tanto los valores negativos como los valores positivos de u y poniendo debajo de la misma gráfica el pozo:
Para estas energías especiales, la curvatura de la parte interior de la función de onda ψ tiene un valor tal que justo en la orilla del pozo se cumple la condición dψ/du = 0. Y el resultado es tal que la amplitud de la oscilación en el exterior tendrá la misma amplitud que la oscilación en el interior.
Probando otras combinaciones de condiciones iniciales para el caso en el cual la energía E de la partícula es mayor que la energía potencial V0, o sea para una partícula que no está ligada permanentemente al pozo de potencial, no tardaremos en descubrir que cuando la amplitud de la onda que describe a la partícula al estar sobre el pozo es diferente de la amplitud de la onda al estar en una región fuera del pozo la amplitud de la onda que corresponde a la región interior al pozo es invariablemente menor que la amplitud de la onda que corresponde a la región exterior al pozo, como podemos apreciarlo en la solución numérica dada arriba para el caso en el cual tenemos un valor de ε igual a 1.4 con una función par. Esto sugiere que, para fines comparativos, podría resultar ilustrativo definir un cociente de amplitudes como el siguiente:
Al referirnos tanto a la “amplitud interior” como la “amplitud exterior”, estaríamos hablando de los valores absolutos (positivos) de dichas amplitudes, esto con la finalidad de darle una vuelta al asunto de que una onda senoidal tiene tantas excursiones positivas como negativas correspondiendo las excursiones negativas a una amplitud negativa. Para el caso en el cual la amplitud de la onda que corresponde al interior del pozo es igual a la amplitud de la onda que corresponde al exterior del pozo, el cociente de amplitudes será igual a la unidad, y como se acaba de mencionar, será menor que la unidad en la gran mayoría de los casos exceptuando los “estados virtuales” para los cuales el cociente es igual a la unidad. Aunque esta es una definición válida, hay otra que puede ser más útil por estar basada en los cuadrados de las amplitudes de las ondas:
Este ajuste tiene que ver directamente con un hecho que estudiaremos más a fondo cuando tratemos el asunto de la interpretación probabilista de la función de onda. Pero el punto relevante a notar aquí es que aquellos valores de ε tales que produzcan una curvatura en la función de onda con la cual se satisfaga la condición dψ/du = 0 justo en u0 = 1/2 serán aquellos que producirán amplitudes iguales en la función de onda tanto en la oscilación exterior como en la oscilación interior al pozo. De cualquier modo, asentaremos aquí que para una partícula libre con una energía E mayor que la energía V0 del pozo de potencial, cuando la amplitud de la oscilación de la función de onda en la región que corresponde al interior del pozo es igual a la amplitud de la oscilación de la función de onda en la región que corresponde al exterior del pozo la probabilidad de encontrar a la partícula sobre la región que corresponde al interior del pozo debe ser igual a la probabilidad de encontrar a la partícula sobre la región que corresponde al exterior del pozo. En cambio, cuando la amplitud de la función de onda cuando la partícula libre está sobre la región que corresponde al interior del pozo es mucho menor que la amplitud de la función de onda cuando la partícula libre está sobre la región que corresponde al exterior del pozo entonces las probabilidades de encontrar a la partícula libre sobre la región que corresponde al interior del pozo deben ser muy bajas. Esto significa que, para los estados virtuales que corresponden a amplitudes iguales en el interior y el exterior para una partícula libre, es como si el pozo de potencial “allí abajo” no existiera. Esto sugiere que el continuum no está carente de cierto grado de estructura. Por otro lado, puesto que la única condición que se requiere para que haya un estado virtual es que se satisfaga la condición dψ/du = 0 justo en u0 = 1/2, habrá un número infinito de estados virtuales. En contraste, sólo hay tres estados ligados.
Puesto que ε se definió simplemente como una relación entre la energía de la partícula y la energía potencial del pozo, esto es:
podemos hacer el siguiente diagrama de niveles energía que resume los resultados obtenidos arriba por la vía numérica para el pozo de potencial con el que hemos estado trabajando:
Al principio de esta entrada, se señaló cómo para ciertos problemas no existe otra solución más que la solución numérica, citándose como ejemplo de ello el potencial Saxon-Woods. Un ejemplo de este potencial así como su gráfica se muestran a continuación:
Como puede apreciarse, este potencial se asemeja mucho al potencial de pozo finito rectangular que acabamos de ver; ya que es “casi” rectangular. Sin embargo, no lo es, y la diferencia es lo suficiente como para que no haya una solución analítica exacta a diferencia del caso del potencial de pozo finito rectangular en donde sí la hay. Para el ejemplo dado, se puede comprobar por la vía de la solución numérica que, utilizando un valor de β = 2ma²V0/ħ² = 64, este potencial Saxon-Woods admite en el interior del pozo únicamente tres estados ligados. El potencial Saxon-Woods es un caso especial de un potencial un poco más general conocido como el potencial Hulthén:
Para un valor específico de q el potencial Hulthén generalizado se reduce a tres tipos muy populares de potenciales. Para q = 0 se reduce al potencial exponencial, para q = 1 se reduce al potencial Hulthén estándard, y para q = -1 se reduce al potencial Saxon-Woods mencionado arriba. A su vez, el potencial Hulthén puede ser considerado emparentado con una familia de potenciales ampliamente utilizados dentro de los cuales cae el potencial Eckart:
para el cual, para variar, se conocen soluciones analíticas exactas en la ecuación de Schrödinger (los parámetros A y B son parámetros arbitrarios, cada uno de los cuales puede ser positivo, negativo, o cero). En el caso específico del potencial Hulthén estándard, una gráfica del mismo nos muestra el gran parecido que tiene con el potencial de Coulomb:
Si hacemos una gráfica poniendo al potencial Coulomb y al potencial Hulthén estándard lado-a-lado, encontraremos que la diferencia entre ambos potenciales se manifiesta de inmediato en las diferencias que hay en los niveles de sus respectivos eigenestados:
Por lo que hemos visto aquí, la gran variedad que hay en la Mecánica Ondulatoria subyace en la amplia variedad de potenciales concebibles que podemos considerar para su estudio (entre los cuales podemos citar al potencial de Yukawa, al potencial de Morse y al potencial Majorana, entre muchos otros). Y en muchos de esos potenciales excepto los más sencillos, no existe solución analítica alguna para la ecuación de Schrödinger, sólo podemos recurrir a la solución numérica. La ventaja de darle una solución numérica a la ecuación de Schrödinger es que, exceptuando algunos casos patológicos, este método casi siempre funciona. La desventaja es que al no llevarse a cabo una solución analítica podemos perder de vista conclusiones importantes, como en el caso de un potencial esféricamente simétrico que resulta esencial para poder entender lo que sucede en el átomo de hidrógeno, y que a su vez resulta esencial para poder entender las propiedades físicas y químicas de los elementos. La abstracción simbólica, a fin de cuentas, puede tener sus ventajas sobre la simple cuestión numérica.
