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 f':\R⟶\R es una función suave. El uso de la derivada no es la notación usual, sin embargo, esta simplifica la definición del Hamiltoniano, la cual requiere de la antiderivada, f, del potencial. La ecuación NLS es utilizada en la modelación de algunos fenómenos físicos que surgen por ejemplo en el estudio de fibras ópticas no lineales, en mecánica cuántica, en el estudio de la física de plasmas, etc. Para una introducción a la ecuación no lineal de Schrödinger y sus aplicaciones referimos al lector al trabajo realizado por Rojas en ¹. Algunos aspectos básicos de los solitones ópticos espaciales han sido resumidos por López y Col. en ²; mientras que, Agüero y Col. presentan un recuento histórico del desarrollo del estudio de las ondas solitarias en ³.

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 ¹¹; 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 ¹² 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 ¹³, 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 ¹⁴, 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 ¹⁵.

Shabat y Zakharov, ¹⁶, 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 ¹⁷; 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, ¹⁰. 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 Th={xk}k=0N de [a, b], con puntos no necesariamente uniformemente distribuidos; en cada celda I _k = (x _k , x _k +1), y para todo t∈[0,  T], se considera el espacio, Ppk(Ik), de polinomios de valor complejo, de grado menor o igual a p _k , con pk≥1; restringidos a la celda I _k . El método LDG, busca, para todo t∈(0,  T], una aproximación (qh(x,t),ψh(x,t))∈Pp(Th)×Pp(Th) de (q(x,t),u(x,t)); tal que las siguientes ecuaciones se verifiquen para toda función de prueba r,v∈Pp(Ik)

donde u¯ representa el conjugado complejo de u. Para cada nodo x _k de la partición Th, los flujos numéricos ψh^ y qh^^ se definen como combinaciones convexas, sea θk∈[0, 1]:

En la práctica, el valor de θk es el mismo en todos los nodos y se consideran únicamente los siguientes valores 0,1/2 y 1. Nos referiremos a estos valores particulares como flujo izquierdo, central y derecho, respectivamente.

En cada nodo, x _k , el término de estabilización, sψh, se define como sigue

donde ηk≥0. Cuando ηk=0, el método es también conocido como LDG de disipación mínima (md-LDG) por sus siglas en inglés); el cual fue analizado para problemas transitorios unidimensionales de convección-difusión lineal con coeficientes constantes, por Castillo ¹⁸ para la versión h y por Castillo y Col.¹⁹ para la version hp. En ²⁰, Cockburn y Dong extendieron el análisis de convergencia para un problema estacionario de convección difusión en mallas triangulares. El primer análisis de convergencia para la versión estabilizada en el caso puramente elíptico, es decir, cuando ηk>0, fue presentada en Castillo y Col.²¹; Perugia y Schötzau ²² realizaron un análisis de convergencia para la version hp; mientras que Cockburn y Col.²³ presentaron un resultado de superconvergencia para mallas Cartesianas. El régimen transitorio multidimensional de convección difusión fue considerado por Cockburn y Dawson ²⁴.

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 qh,  ψh satisface para todo t∈(0, T] y para todo r, v ∈Pp Th×Pp Th las siguientes ecuaciones

En lo que sigue, se denotará en mayúsculas, el vector de coeficientes de la expansión de una función de Pp Th, en una base específica; en nuestra implementación se utilizaron los polinomios de Legendre. El problema discreto se reduce a un sistema global de m ecuaciones diferenciales ordinarias:

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 ψ↦-ψxx; la cual se probó en ²¹ ser una matriz de coeficientes reales simétrica no negativa; y fh(⋅) es la representación discreta del operador complejo de potencial no lineal, tal que, para todo v∈Pp(Th) se tiene

donde U;V representa el producto interior usual en Cm.

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 ψ'h, no es casual, ya que es propia de la formulación mixta del método LDG, Ecs. (4a) y (4b). Además, se ha demostrado que en mallas Cartesianas, ^18,23, esta variable puede alcanzar orden de convergencia óptimo, O(hp+1); y en ¹³ se observó, numéricamente, una mejor precisión que su análogo discreto ψ'h para mallas no estructuradas.

Proposición 1 (Conservación de energía) La soluciónqhx, t,ψh x, tde la formulación semi-discreta del método LDG, Ecs. (8a) y (8b), conserva energía, para la ecuación no lineal de Schrödinger (1) con condiciones de borde periódicas.

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 AΨh;Ψh  son, ambas, cantidades reales, se deduce que

Proposición 2 (Conservación del Hamiltoniano) La solución(qh(x,t),ψh(x,t))de la formulación semi-discreta del método LDG, Ecs. (8a) y (8b), conserva el Hamiltoniano discreto (11).

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 ²⁵ para una ecuación no lineal de Klein-Gordon; Delfour, Fortin y Payré ⁴ 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 ⁹, 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: τ denotará el incremento en tiempo; el cual se considera constante;

y

Siguiendo Sanz-Serna ⁹, 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 Fh(⋅) es la discretización del operador no lineal tal que para todo v∈PpTh se tenga

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ónqh, ψhobtenida por el método de Crank-Nicolson modificado y con discretización espacial LDG conserva energía.

Prueba. Sea v=ψhn+(1/2), entonces, de acuerdo a la definición del método MCN se tiene

Nótese que por la definición de la función auxiliar F(⋅,⋅) se tiene

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ónqh, ψhobtenida por el método de Crank-Nicolson modificado y con discretización espacial LDG conserva el Hamiltoniano discreto.

Prueba. Sea v=δτψhn+(1/2), entonces, considerando, esta vez, la parte real se tiene

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 τ=(1/2)h(p+1)/2 para equiparar los órdenes de aproximación en la variable espacial y temporal. En cada paso de tiempo se utilizó el método de Newton para resolver el problema no lineal. Puesto que el sistema lineal no es de gran tamaño, este se resuelve por el método de eliminación de Gauss. Para cada ejemplo se consideró una cuadratura Gaussiana con suficientes puntos de cuadratura de tal manera que el cálculo del operador no lineal fuese lo más preciso posible.

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

iψt=-ψxx-2|ψ|2ψ, (21)

y condición inicial

ψox=sechx+10 exp(2i(x+10)).

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.

Figura 1. Evolución del solitón en el intervalo de tiempo [0, 10], para p = 2. 

En la Fig. 2 se compara la evolución de la la energía y del Hamiltoniano para los flujos izquierdo, central y derecho.

Figura 2. Comparación con respecto a los flujos izquierdo, central y derecho en el problema del solitón. 

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.

Figura 3. Comparación con respecto al grado de aproximación para el problema del solitón. 

Figura 4. Evolución de la colisión en el intervalo de tiempo [0, 5], para p = 2. 

Figura 5. Comparación con respecto a los flujos izquierdo, central y derecho para el problema de la doble colisión. 

Figura 6. Comparación con respecto al grado de aproximación para el problema de doble colisión. 

Figura 7. Comportamiento del error en el Hamiltoniano con el método original de Crank-Nicolson (CN) y su version modificada (MCN). 

Figura 8. Comparación con respecto a los flujos para el problema del potencial de alto orden. 

Figura 9. Evolución de la solución del problema del potencial de alto orden en el intervalo de tiempo [0, 5], para p = 2. 

Figura 10. Efecto de la tolerancia en el método de Newton en la variación del error en la energía y el Hamiltoniano. 

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.

Cuadro I. Tiempo de ejecución para aproximaciones de grado p = 3, 9, 15 e igual número de grados de libertad (1600). 

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.