PACS: 02.60.Cb; 02.70.Dh; 05.45.Yv
1. Introducción
Se considera la siguiente ecuación de Schrödinger no lineal (NLS por sus siglas en inglés)
en el dominio de espacio y tiempo (a, b) × (0, T] y condiciones periódicas; donde el potencial
Aunque existen soluciones explícitas para cierto tipo de potencial; esto no es cierto en general; por lo que el estudio de métodos numéricos para esta ecuación es imprescindible. Si bien es cierto, diversas discretizaciones espaciales han sido utilizadas; como por ejemplo diferencias finitas 4,5,6,7,8, elemento finito 9,10 y métodos espectrales 11; este trabajo se concentra en una técnica de discretización espacial particular, llamada “Local Discontinuous Galerkin” (LDG por sus siglas en inglés), la cual fue propuesta, originalmente, por Cockburn y Shu, en 12 para problemas transitorios de difusión y convección no lineal. Este es un método particular dentro de la clase de métodos de elemento finito conocida como métodos Galerkin Discontinuo, cuya particularidad es la de no imponer continuidad entre celdas adyacentes, lo cual resulta apropiado para métodos adaptativos en espacio y grado de aproximación; es decir, para la versión hp del método de elemento finito. En 13, Castillo realizó un estudio comparativo para un problema elíptico en 2D, el cual reveló cierta superioridad del método LDG en comparación a otros métodos discontinuos similares.
En 14, Xu and Shu analizaron por primera vez el método LDG aplicado a la ecuación NLS (1) y a un sistema acoplado de ecuaciones de Schrödinger no lineales. Aunque en dicho artículo no se obtuvieron estimados de error óptimos, estos fueron obtenidos posteriormente en 15.
Shabat y Zakharov, 16, utilizaron por primera vez la técnica de dispersión inversa para mostrar que existe un conjunto numerable de leyes de conservación para la ecuación no lineal de Schrödinger. Dicho conjunto de leyes de conservación es una condición necesaria de integrabilidad del problema como se mostró en 17; siendo la conservación de la energía (2a) y del Hamiltoniano (2b), de las más importantes.
La conservación de estos invariantes no debe ser desestimada, pues está íntimamente relacionada con la estabilidad de la solución como fue mostrado por Griffiths, Mitchell y Morris, 10. La conservación de la energía garantiza que la aproximación se mantenga acotada en cada paso de tiempo; descartando el posible problema de “blow-up” (crecimiento al infinito en tiempo finito). Por otro lado, estos invariantes representan cantidades de interés físico, por lo que es deseable conservarlas a nivel discreto. Por ejemplo, en mecánica cuántica estos representan masa y energía, respectivamente; mientras que, en óptica, potencia y energía electromagnética.
2. Método LDG para la ecuación NLS
En la formulación del método LDG, la ecuación no lineal de Schrödinger Ec. (1); se reescribe como un sistema de ecuaciones diferenciales parciales de primer grado
Dada una partición
donde
En la práctica, el valor de
En cada nodo, x
k
, el término de estabilización,
donde
Para efectos de implementación, es más adecuado considerar la representación matricial D, B y S, respectivamente, de las siguientes formas sesquilineales
las cuales se deducen, de la formulación del método Ecs. (4a) y (4b), sumando sobre todas las celdas. La solución
En lo que sigue, se denotará en mayúsculas, el vector de coeficientes de la expansión de una función de
donde M es la matriz de masa, en el caso unidimensional M = D, y A = B
T
D
-1
B+S no es más que la representación de la discretización por el método LDG del operador diferencial de segundo orden
donde
3. Problema semi-discreto
A continuación se analiza la conservación de la energía y del Hamiltoniano para el método semi-discreto; es decir, asumiendo la variable temporal continua y utilizando el método LDG como discretización espacial. Para ello se deriva el Hamiltoniano discreto a partir de la formulación variacional de la siguiente forma:
Esta cantidad debe interpretarse como una aproximación del Hamiltoniano (2b). La aparición de la variable auxiliar q
h
, en lugar de
Proposición 1 (Conservación de energía)
La solución
Prueba. Nótese que
Además, por Ec. (9) se tiene
Por lo tanto considerando la parte imaginaria en la ecuación anterior y observando que
y
Proposición 2 (Conservación del Hamiltoniano)
La solución
Prueba. Por (9) se tiene
Por lo que, considerando la parte real de la expresión anterior y notando que M es una matriz simétrica definida positiva, se deduce la siguiente ecuación
Como A es una matriz real simétrica entonces
Por otro lado, nótese que también se tiene
Substituyendo las expresiones obtenidas en las Ecs. (14) y en la Ec. (15), se obtiene la relación deseada
4. Discretización temporal
A continuación se analiza las propiedades de conservación de la energía y del Hamiltoniano en tiempo discreto en combinación con la discretización espacial discontinua LDG. Para ello, se considera como método de discretización en tiempo el método de Crank-Nicolson modificado que ha sido ampliamente utilizado en el pasado. Motivados por el esquema propuesto por Strauss and Vazquez en 25 para una ecuación no lineal de Klein-Gordon; Delfour, Fortin y Payré 4 presentan una modificación para un caso particular de la Ec. (1), utilizando como discretización espacial el método de diferencias finitas; mientras que, Sanz-Serna en 9, analiza la convergencia de este método para el problema (1) utilizando el método de elemento finito clásico; es decir, aproximaciones continuas.
Para la descripción de estos métodos se adoptará la siguiente notación:
y
Siguiendo Sanz-Serna 9, se considera la siguiente función auxiliar
El método modificado de Crank-Nicolson (MCN) en combinación con la discretización espacial LDG, (9) se reduce a
donde
Esta formulación difiere de la forma original del método de Crank-Nicolson (CN) en la representación del término no lineal, el cual se calcularía de la siguiente forma
Como se mostrará a continuación, esta modificación posee las propiedades de conservación deseadas para los análogos discretos de la energía y del Hamiltoniano.
Proposición 3 (Conservación de energía)
La solución
Prueba. Sea
Nótese que por la definición de la función auxiliar
Considerando, únicamente, la parte imaginaria en la Ec. (17), y puesto que la matriz A es simétrica, se obtiene
Proposición 4 (Conservación del Hamiltoniano)
La solución
Prueba. Sea
Siendo A una matriz simétrica se tiene
Por otro lado también se tiene
El resultado se obtiene por substitución de las expresiones de las Ecs. (19) y (20) en la Ec. (18). □
5. Experimentos numéricos
Con el propósito de validar los resultados teóricos demostrados en las secciones anteriores, se presenta, a continuación, una serie de experimentos numéricos. La implementación se realizó en el ambiente de Matlab, utilizando aritmética compleja; es decir, que la Ec. (1) no se descompuso como un sistema no lineal en función de la parte real e imaginaria. Se utilizó un incremento en tiempo constante, donde
Debido al movimiento de la solución de onda a lo largo del tiempo, la selección de una malla unifome es natural. Se utilizaron mallas con una distribución uniforme de tamaño h = 0.25.
5.1 Evolución de un solitón
Como primer ejemplo se considera la evolución de un solitón el cual se describe mediante la siguiente ecuación
y condición inicial
La evolución del solitón, en el dominio [a, b] × [0, T] = [-100, 100] × [0, 10]. se muestra en la Fig. 1, para aproximaciones de grado 2 y flujo central.
En la Fig. 2 se compara la evolución de la la energía y del Hamiltoniano para los flujos izquierdo, central y derecho.
Aunque solamente se muestran resultados para polinomios de grado p = 2, se obtuvieron resultados similares para distintos grados de aproximación. En la Fig. 3 se compara la evolución del error en la energía y el Hamiltoniano, para polinomios de grado p = 1 - 4, con respecto al valor inicial de estas cantidades. Si bien es cierto la teoría predice conservación, desde el punto de vista numérico, es importante monitorear el error. Nótese que, aunque se observan pequeñas fluctuaciones en el error, las cuales se deben a la falta de precisión en la solución de los problemas no lineales en cada iteración; estas ocurren en escala de precisión de máquina, 10-14 y 10-13, respectivamente; lo cual ratifica, nuevamente, los resultados teóricos. Este comportamiento también se muestra en el siguiente experimento y se detalla aún más en las explicaciones de la Fig. 10.
5.2. Doble colisión de solitones
En este ejemplo se examina la conservación de los invariantes discretos en el transcurso de una colisión de dos solitones.
El modelo está dado por la Ec. (21) del ejemplo anterior, con la siguiente condición inicial
La evolución de la doble colisión, en el intervalo de tiempo [0, 5], se muestra en la Fig. 4, para aproximaciones de grado 2 y flujo central.
La energía y el Hamiltoniano discreto, para este problema, se muestran en la Fig. 5, para todos los flujos y aproximaciones de grado 4. La comparación con respecto al grado de aproximación se muestra en la Fig. 6, la cual muestra conservación con un error de orden 10-13 con respecto a las aproximaciones iniciales de estas cantidades.
En la Fig. 7 se muestra el comportamiento del error en el Hamiltoniano discreto para el método de Crank-Nicolson (CN) y su versión modificada (MCN), para el experimento de la doble colisión. Se aprecia la falta de conservación del Hamiltoniano al utilizar método original de Crank Nicolson. La discretización espacial se realizó con polinomios cuadráticos y utilizando el flujo central.
5.3 Potencial de alto orden
En este experimento se muestra la conservación de ambos invariantes para potenciales de alto orden con la ecuación
en [-30, 30] × (0, 5], y condición inicial
5.4. Aproximaciones de alto orden
El Cuadro I muestra la ventaja de utilizar polinomios de alto orden. Como ejemplo particular se consideró el problema de doble colisión de solitones. Para cada grado se utiliza una malla de tal manera que el número total de grados de libertad sea de 1600. Nótese que entre mayor grado de aproximación, menor la cantidad total de celdas; y mucho menor tiempo de ejecución. En este ejemplo se alcanzó un factor de 10, en la reducción del tiempo de ejecución, con respecto al de aproximaciones cúbicas. Los cálculos se realizaron en una Laptop DELL con procesador Intel Core i5 y 4Gb de memoria RAM, bajo el sistema operativo Linux. Se ha incluido el número promedio de iteraciones en cada paso de tiempo realizado por el método de Newton. Nótese el aumento del mismo con el grado de aproximación.
Es importante resaltar el efecto de la tolerancia del criterio de convergencia en el método de Newton sobre el cálculo de los invariantes. En la Fig. 10 se muestra la evolución del error en la energia y el Hamiltoniano para tres valores distintos en la tolerancia. La gráfica muestra una pérdida de precisión sustancial en el cálculo de ambos invariantes, a medida que dicha tolerancia aumenta.
6. Conclusiones
La conservación de la energía y del Hamiltoniano para una ecuación de Schrödinger no lineal general, fue demostrada para discretizaciones espaciales basadas en el método LDG. Además, se demostró la conservación del análogo discreto de estas cantidades para el problema completamente discreto utilizando el método modificado de Crank-Nicolson como técnica de integración en tiempo. Los experimentos numéricos validaron el análisis teórico presentado.
Debido al costo computacional excesivo que requiere la iteración de Newton, estamos estudiando otros tipos de técnicas que permitan acelarar esta iteración; las cuales pueden ser de gran beneficio, particularmente en problemas multidimensionales. Estas serán presentadas en un artículo futuro.