7.2: Discretización de ecuaciones diferenciales parciales

Con derivadas discretizadas, las ecuaciones diferenciales pueden formularse como sistemas discretos de ecuaciones. Lo discutiremos usando un ejemplo específico: la discretización de la ecuación de onda de Schrödinger independiente del tiempo en 1D.

7.2.1 Derivar una ecuación de diferencia finita

La ecuación de onda de Schrödinger independiente del tiempo 1D es la ecuación diferencial ordinaria de segundo orden

\[-\frac{\hbar^2}{2m} \, \frac{d^2\psi}{dx^2} + V(x) \psi(x) = E \psi(x),\]

donde\(\hbar\) está la constante de Planck dividida por\(2\pi\),\(m\) es la masa de la partícula,\(V(x)\) es el potencial,\(\psi(x)\) es la función de onda cuántica de un estado propio de energía de la partícula, y\(E\) es la energía correspondiente. La ecuación diferencial suele tratarse como un problema propio, en el sentido de que se nos da\(V(x)\) y buscamos encontrar los posibles valores de la función propia\(\psi(x)\) y el valor propio de la energía\(E\). Para mayor comodidad, adoptaremos unidades donde\(\hbar=m = 1\):

\[-\frac{1}{2}\, \frac{d^2\psi}{dx^2} + V(x) \psi(x) = E \psi(x).\]

Para discretizar esta ecuación diferencial, simplemente la evaluamos en\(x = x_n\):

\[-\frac{1}{2}\, \psi''(x_n) + V_n \psi_n = E \psi_n,\]

donde, por concisión, denotamos

\[V_n \equiv V(x_n).\]

Luego reemplazamos la segunda derivada\(\psi''(x_n)\) por una aproximación discreta, específicamente la regla de tres puntos:

\[-\frac{1}{2h^2}\, \Big[\psi_{n+1} - 2\psi_n + \psi_{n-1} \Big] + V_n \psi_n = E \psi_n.\]

Este resultado se denomina ecuación de diferencia finita, y sería válido para todos\(n\) si el número de puntos de discretización es infinito. Sin embargo, si hay un número finito de puntos de discretización\(\{x_0, x_1, \dots, x_{N-1}\}\), entonces la fórmula de diferencia finita falla en los puntos límite\(n=N-1\),\(n=0\) y, donde implica el valor de la función en los puntos “inexistentes”\(x_{-1}\) y\(x_{N}\). Veremos cómo manejar este problema en la siguiente sección.

Dejando de lado los límites, la ecuación de diferencia finita describe una ecuación matricial:

\[\left\{-\frac{1}{2h^2}\begin{bmatrix} \ddots & \ddots \\ \ddots & -2 & 1\\ & 1 & -2 & 1\\ & & 1 & -2 & \ddots \\ & & & \ddots & \ddots & \end{bmatrix} + \begin{bmatrix} & \ddots & \\ & & V_{n-1} & \\ & & & V_n & \\ & & & & V_{n+1} & \\ & & & & & \ddots \end{bmatrix} \right\}\begin{bmatrix}\vdots \\ \psi_{n-1} \\ \psi_n \\ \psi_{n+1} \\ \vdots\end{bmatrix} = E \begin{bmatrix}\vdots \\ \psi_{n-1} \\ \psi_n \\ \psi_{n+1} \\ \vdots\end{bmatrix}.\]

El operador de segunda derivada está representado por una matriz tridiagonal con\(-2\) en cada elemento diagonal, y\(1\) en los elementos directamente por encima y por debajo de la diagonal. El operador potencial está representado por una matriz diagonal, donde los elementos a lo largo de la diagonal son los valores del potencial en cada punto de discretización. De esta manera, la ecuación de onda de Schrödinger se reduce a un problema de autovalor discreto.

7.2.2 Condiciones de contorno

Ahora tenemos que averiguar cómo manejar los límites. Supongamos que\(\psi(x)\) se define sobre un intervalo finito,\(a \le x \le b\). Como recordamos de la teoría de ecuaciones diferenciales, la solución a una ecuación diferencial no está totalmente determinada por la propia ecuación diferencial, sino también por las condiciones límite que se imponen. Así, tenemos que especificar cómo\(\psi(x)\) se comporta en los puntos finales del intervalo. Mostraremos cómo se hace esto para un par de las condiciones de límite más comunes; otras elecciones de condiciones de límite se pueden manejar usando el mismo tipo de razonamiento.

Condiciones de contorno de Dirichlet

Bajo las condiciones de límite de Dirichlet, la función de onda desaparece en los límites:

\[\psi(a) = \psi(b) = 0.\]

Físicamente, estas condiciones de límite se aplican si dejamos que el potencial explote en las regiones externas,\(x>b\) y\(x<a\), por lo tanto, obligando a que la función de onda se limite estrictamente al intervalo\(a \le x \le b\).

Aún no hemos expuesto cómo se\(\{x_0, \dots, x_{N-1}\}\) distribuyen los puntos de discretización dentro del intervalo; tomaremos esta decisión en tándem con la implementación de las condiciones límite. Considera el primer punto de discretización,\(x_{0}\), donde sea que esté. La ecuación de diferencia finita en este punto es

\[-\frac{1}{2h^2}\, \Big[\psi_{-1} - 2\psi_0 + \psi_{1} \Big] + V_0 \psi_0 = E \psi_0.\]

Esto implica la función de onda en\(x_{-1}\), que se encuentra justo fuera de nuestro conjunto de puntos de discretización. Pero si elegimos los puntos de discretización para que\(x_{-1}=a\), luego\(\psi_{-1}=0\) bajo condiciones de límite de Dirichlet, así la fórmula de diferencia finita anterior se reduzca a

\[-\frac{1}{2h^2}\, \Big[- 2\psi_0 + \psi_{1} \Big] + V_0 \psi_0 = E \psi_0.\]

En cuanto al otro límite, la ecuación de diferencia finita en\(x_{N-1}\) implica\(\psi_{N}\). Si elegimos los puntos de discretización para que\(x_{N}=b\), entonces la fórmula de diferencia finita se convierta en

\[-\frac{1}{2h^2}\, \Big[ \psi_{N-2} - 2\psi_{N-1} \Big] + V_{N-1} \psi_{N-1} = E \psi_{N-1}.\]

A partir de esto, se concluye que los puntos de discretización deben estar igualmente espaciados, con\(x_{0}\) una distancia\(h\) a la derecha del límite izquierdo\(a\) y\(x_{N-1}\) una distancia\(h\) a la izquierda del límite derecho\(b\). Esto se muestra en la siguiente figura:

Figura\(\PageIndex{1}\): Posición de los puntos de discretización para las condiciones de contorno de Dirichlet en\(x=a\) y\(x=b\).

Dado que hay puntos de\(N\) discretización, el intervalo debe contener\((N+1)\) múltiplos de\(h\). Por lo tanto,

\[h = \frac{b - a}{N + 1} \;\; \Rightarrow \;\; x_n \,=\, a + h (n+1) \,=\, \frac{a(N-n)+b(n+1)}{N+1}.\]

Habiendo hecho las elecciones anteriores, la ecuación matricial se convierte en

\[\left\{-\frac{1}{2h^2}\begin{bmatrix} -2 & 1 \\ 1 & -2 & \ddots \\ & \ddots & \ddots & 1 \\ & & 1 & -2\end{bmatrix} + \begin{bmatrix} V_0 \\ & V_1 \\ & & \ddots\\ & & & V_{N-1} \end{bmatrix} \right\}\begin{bmatrix}\psi_0 \\ \psi_1 \\ \vdots \\ \psi_{N-1}\end{bmatrix} = E \begin{bmatrix}\psi_0 \\ \psi_1 \\ \vdots \\ \psi_{N-1}\end{bmatrix}.\]

Puedes comprobar por ti mismo que la primera y última fila de esta ecuación son las ecuaciones correctas de diferencia finita en los puntos límite, correspondientes a las condiciones de contorno de Dirichlet.

Condiciones de límite de Neuman

Las condiciones de límite de Neumann son otra opción común de condiciones de límite. Anotan que los primeros derivados desaparecen en los límites:

\[\psi'(a) = \psi'(b) = 0.\]

Un ejemplo de tal condición límite se encuentra en la electrostática, donde la primera derivada del potencial eléctrico va a cero en la superficie de una superficie metálica cargada.

Seguimos la misma estrategia que antes, averiguando los puntos de discretización en tándem con las condiciones límite. Consideremos nuevamente la ecuación de diferencia finita en el primer punto de discretización:

\[-\frac{1}{2h^2}\, \Big[\psi_{-1} - 2\psi_0 + \psi_{1} \Big] + V_0 \psi_0 = E \psi_0.\]

Para implementar la condición de que la primera derivada se desvanezca en el límite, invocamos la regla del punto medio. Supongamos que el punto límite\(x=a\) cae entre los puntos\(x_{-1}\) y\(x_{0}\). Entonces, de acuerdo con la regla del punto medio,

\[\frac{\psi_0 - \psi_{-1}}{h} \approx \psi'(a) = 0.\]

Con esta elección, por lo tanto, podemos hacer el reemplazo\(\psi_{-1} = \psi_0\) en la ecuación de diferencia finita, que luego se convierte

\[-\frac{1}{2h^2}\, \Big[- \psi_0 + \psi_{1} \Big] + V_0 \psi_0 = E \psi_0.\]

Del mismo modo, para aplicar la condición de límite de Neumann en\(x=b\), dejamos que el límite caiga entre\(x_{N-1}\) y\(x_{N}\), de modo que la ecuación de diferencia finita se convierta

\[-\frac{1}{2h^2}\, \Big[ \psi_{N-2} - \psi_{N-1} \Big] + V_{N-1} \psi_{N-1} = E \psi_{N-1}.\]

La distribución resultante de los puntos de discretización se muestra en la siguiente figura:

Figura\(\PageIndex{2}\): Posición de los puntos de discretización para las condiciones límite de Neumann en\(x=a\) y\(x=b\).

A diferencia del caso Dirichlet, el intervalo contiene\(N\) múltiplos de\(h\). De ahí que obtengamos una fórmula diferente para las posiciones de los puntos de discretización

\[h = \frac{b - a}{N} \;\; \Rightarrow \;\; x_n \,=\, a + h \left(n+\frac{1}{2}\right) \,=\, \frac{a(N-n-\tfrac{1}{2})+b(n+\tfrac{1}{2})}{N}.\]

La ecuación matricial es:

\[\left\{-\frac{1}{2h^2}\begin{bmatrix} -1 & 1 \\ 1 & -2 & \ddots \\ & \ddots & \ddots & \ddots \\ && \ddots & -2 & 1 \\ & & & 1 & -1\end{bmatrix} + \begin{bmatrix} V_0 \\ & V_1 \\ & & \ddots\\ & & & V_{N-1} \end{bmatrix} \right\}\begin{bmatrix}\psi_0 \\ \psi_1 \\ \vdots \\ \psi_{N-1}\end{bmatrix} = E \begin{bmatrix}\psi_0 \\ \psi_1 \\ \vdots \\ \psi_{N-1}\end{bmatrix}. \]

Debido a las condiciones de límite de Neumann y la regla del punto medio, la matriz tridiagonal tiene en\(-1\) lugar de\(-2\) en sus entradas de esquina. Nuevamente, se puede verificar que la primera y última fila de esta ecuación matricial corresponden a las ecuaciones correctas de diferencia finita para los puntos límite.