Análisis numérico de soluciones autosemejantes de un flujo dispersivo de curvas planas

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) An´alisis num´erico de soluciones autosemejantes de un flujo dispersivo de curvas planas Francisco de la Hoz1 1Dpto. de Matem´atica Aplicada, Estad´ıstica e Investigaci´on Operativa, Facultad de Ciencia y Tecnolog´ıa, Universidad del Pa´ıs Vasco - Euskal Herriko Unibertsitatea, Aptdo. 644, 48080 Bilbao (Spain). E-mail: [email protected]. Palabras clave: Flujo geom´etrico de la KdV, estudio num´erico Resumen En [2], Perelman y Vega estudian el siguiente flujo geom´etrico de curvas planas reversible en el tiempo, que puede desarrollar singularidades en tiempo finito    zt=−zsss +3 2¯zsz2 ss, |zs|2= 1, t 6= 0. (1) con sel par´ametro de arco. Denotando por kla curvatura de z, esta satisface la mKdV kt+ksss +3 2k2ks= 0.(2) Perelman y Vega consideran soluciones autosemejantes de (2) de la forma k(s, t) = 2 (3t)1/3uµs (3t)1/3¶, t > 0; (3) lo cual conduce a estudiar la EDO uxx −xu + 2u3=µ, x ∈R, µ ∈R.(4) En esta comunicaci´on, consideraremos µ= 0. Aunque necesitamos conocer u(0) y ux(0), para resolver (4), imponiendo limx→∞ u(x) = 0, los datos iniciales para (4) forman una familia uniparam´etrica, que obtendremos num´ericamente. Adem´as, daremos evidencia num´erica de que las soluciones de (4) correspondientes satisfacen −π 2≤Z∞ −∞ u(x)dx ≤π 2.(5) Por (3), a cada ule corresponde un dato inicial zpara (1) en t= 1. Considerando datos iniciales sin intersecciones, mostraremos num´ericamente su evoluci´on, as´ı como la formaci´on de una singularidad en t= 0. 1 Francisco de la Hoz 1. Introducci´on En esta comunicaci´on, consideraremos el siguiente flujo geom´etrico reversible en el tiempo de curvas planas que puede desarrollar singularidades en tiempo finito    zt=−zsss +3 2¯zsz2 ss, |zs|2= 1, t 6= 0, (1) con sel par´ametro de arco. Este flujo es reversible en el tiempo porque si z(s, t) es soluci´on, tambi´en lo es z(−s, −t). Dicho flujo, fue obtenido por Goldstein y Petrich [1] al considerar la evoluci´on de un parche de vorticidad en el plano sujeto a las ecuaciones de Euler. Tambi´en es conocido como el flujo geom´etrico de la KdV. La evoluci´on del flujo viene determinada por la evoluci´on de su curvatura, k(s, t), salvo un movimiento r´ıgido que varia con el tiempo y que ser´a fijado por las condiciones iniciales. Dicha curvatura satisface la KdV modificada kt+ksss +3 2k2ks= 0; (2) obs´ervese que el ´ultimo sumando es una derivada perfecta kt=−ksss −1 2(k3)s; (3) por lo cual Z+∞ −∞ k(s, t)ds (4) es una cantidad conservada. Asimismo, definiendo θ(s, t) = θ(−∞, s) + Zs −∞ k(s0, t)ds0,(5) la ecuaci´on (3) se convierte en θt(s, t) = −θsss(s, t)−1 2(θs)3(s, t).(6) Siguiendo a Perelman y Vega [2] a lo largo de esta comunicaci´on, consideraremos soluciones autosemejantes de la KdV modificada de la forma k(s, t) = 2 (3t)1/3uµs (3t)1/3¶, t > 0.(7) lo que nos lleva a estudiar la EDO [2] uxx −xu + 2u3=µ, x ∈R, µ ∈R,(8) con µuna constante de integraci´on. En esta comunicaci´on consideraremos el caso µ= 0. A partir de u(x) obtendremos k(s, 1) y, por tanto, θ(s, 1), que utilizaremos como dato inicial para (6), calculando la evoluci´on num´erica de θ(s, t), yendo hacia atr´as en el tiempo. 2 An´alisis num´erico de soluciones de un flujo dispersivo Es inmediato obtener z(s, t) a partir de θ(s, t), haciendo zs(s, t) = exp(iθ(s, t)), salvo un movimiento r´ıgido que viene determinado por z(0, t) = −2(3t)1/3£iu0(0) + u2(0)¤zs(0, t),(9) con zs(0, t) una constante de m´odulo unidad que no depende del tiempo. Tanto θ(s, t) como k(s, t) y z(s, t) desarrollar´an una singularidad en t= 0. Debido a que las u(x) consideradas, que satisfacen limx→∞ u(x)→0, forman una familia uniparam´etrica, tenemos que las z(s, t) correspondientes son tambi´en familias uniparam´etricas de soluciones de (1). 2. Integraci´on de u(x) Volvamos a uxx =xu −2u3, x ∈R.(10) Para integrar esta ordinaria de segundo orden, la reescribimos como (ux=v, vx=xu −2u3,(11) y utilizamos un Runge-Kutta de cuarto orden, necesitando dos datos iniciales u(0) y ux(0) = v(0). Si ues soluci´on de (10), entonces, para xnegativos suficientementes grandes, u(x) tiene el mismo comportamiento oscilatorio que las soluciones de la ecuaci´on de Airy uxx −xux= 0. Asimismo, cuando x→ ∞, las soluciones de (10) se caracterizan por ser muy sensibles respecto a peque˜nas variaciones de datos iniciales, que pueden hacer que pasemos de soluciones u(x)→ ∞, cuando x→ ∞ a soluciones u(x)→ −∞, cuando x→ ∞ y a la inversa. Consideremos, por ejemplo, u(0) = 0.024, as´ı como cuatro elecciones posibles de ux(0) e integremos num´ericamente (10). Las soluciones para ux(0) = −0.018 y ux(0) = −0.017, pintadas con trama discontinua (figura 1), tienden respectivamente a −∞ y +∞; asimismo, las soluciones ux(0) = −0.0175 yux(0) = −0.0174, pintadas con trama continua, tienden respectivamente a −∞ y∞. No obstante, vemos que para estos dos valores, la explosi´on se produce algo m´as tarde. −5 0 5 10 15 20 −3 −2 −1 0 1 2 3 u(x) Figura 1: Dependencia sensitiva de u(x) respecto de las condiciones iniciales. 3 Francisco de la Hoz Si continu´aramos con el proceso, entre ux(0) = −0.0175 y ux(0) = −0.0174 existe un ux(0) tal que limx→∞ u(x) = 0. Dicho valor, obtenido mediante Newton-Rapson, es aproximadamente ux(0) = −0.0174881944 ... En la figura 1, por cuestiones de escala, no se aprecian bien las oscilaciones en el eje real negativo; la figura 2, para dicho ux(0) limite, muestra dichas oscilaciones. −20 −15 −10 −5 0 5 10 15 20 −0.02 −0.01 0 0.01 0.02 0.03 u(x), u(0) = .024, ux(0) = −.0174881944... Figura 2: u(x), con limx→∞ u(x) = 0. Debido a la dependencia sensitiva respecto de las condiciones iniciales, tan s´olo podemos calcular num´ericamente tantos decimales del ux(0) deseado como el εde la m´aquina, lo que implica que ´unicamente estamos retrasando el tiempo de explosi´on de u(x). Dicho de otra forma, las condiciones que hacen lim x→∞ u(x) = 0,(12) son sumamente inestables. Sin embargo, para nuestros prop´ositos, basta con considerar u(x)≡0 a partir de un xsuficientemente grande, ya que u(x)→0 de forma exponencial. Puesto que u(x) es integrable entre (−∞,0], debido a que las oscilaciones en x < 0 se cancelan, la condici´on (12) es equivalente a ¯¯¯R+∞ −∞ u(x)dx¯¯¯<∞. Figura 3: Pares (u(0), ux(0)), con lim x→∞u(x) = +∞(negro) y lim x→∞u(x) = −∞ (blanco). 4 An´alisis num´erico de soluciones de un flujo dispersivo En la figura 3, los puntos de la regi´on en negro, tomados como condiciones iniciales (u(0), ux(0)) de (10), hacen que las soluciones de (10) correspondientes u(x)→+∞, cuando x→ ∞. Para los puntos de la regi´on en blanco, las soluciones de (10) correspondientes u(x)→ −∞, cuando x→ ∞. La frontera entre ambas regiones son los pares (u(0), ux(0)) cuyas soluciones satisfacen (12), que en adelante denominaremos pares admisibles. Dicha frontera es una curva conexa no acotada de dimensi´on uno que divide al plano en dos mitades. Cada uno de los puntos de dicha curva determinan una ´unica k(s, t), una ´unica θ(s, t) y una ´unica z(s, t) autosemejantes que son soluciones de (2), (6) y (1), respectivamente. Por tanto, tenemos sendas familias uniparam´etricas de soluciones de (2), (6) y (1) que desarrollan una singularidad en tiempo finito. 3. Un m´etodo num´erico para la evoluci´on de (1) Las k(s, t), θ(s, t) y z(s, t) consideradas se caracterizan por k(−∞, t) = 0, k(+∞, t) = 0, θ(−∞, t) = θ−, θ(+∞, t) = θ+, lim s→−∞ |z(s, t)−seiθ− |= 0,lim s→+∞|z(s, t)−seiθ− |= 0. El utilizar θ(s, t) en la evoluci´on, en lugar de k(s, t), tiene la ventaja de preservar (4) y |zs(s, t)|= 1. La EDP para θ, con condiciones de frontera, es        θt(s, t) = −θsss(s, t)−1 2(θs)3(s, t), s ∈R, θ(−∞, t) = θ−, θ(+∞, t) = θ+. (13) Al no poder considerar todo s∈R, nos limitaremos a s∈[sa, sb], con sa<< −1 y sb>> 1        θt(s, t) = −θsss(s, t)−1 2(θs)3(s, t), s ∈[sa, sb], θ(sa, t) = θ−, θ(sb, t) = θ+. (14) Esta ´ultima ecuaci´on se puede integrar mediante un Runge-Kutta de cuarto orden con factor integrante. 3.1. C´alculo de k(s, 1) yθ(s, 1) Ya que kno es peri´odica, consideraremos s∈[sa, sb], sa<< −1, sb>> 1, de tal forma que k(sa,1) sea suficientemente peque˜no y k(s, 1) = 0, ∀s > sb−ε. Consideramos N1+ 1 puntos equiespaciados en [sa, sb], sj=sa+j∆s, ∆s=sb−sa N1 , j = 0, . . . , N1. Las elecciones de sa,sbyN1se hacen de tal forma que 0 ∈ {sj}N1 j=0. 5 Francisco de la Hoz Para integrar θ(s, 1) con exactitud espectral a partir de k(s, 1), obtenida a su vez de (7), interesa forzar que k(s, 1) sea peri´odica. Para ello, localizamos el sjcon menor ´ındice tal que k(sj,1) >0 y ks(sj,1) >0, es decir, escogemos el primer sjen el que k(sj,1) sea positiva y creciente. Acoplaremos a k(s, 1) en s=sjuna funci´on que tenga un contacto de orden uno y decaimiento exponencial; posteriormente, mediante un filtro, aumentaremos la regularidad del contacto. La funci´on que acoplamos a k(s, 1) en s=sjes kleft(s, 1) = k(sj,1) exp ·ks(sj,1) k(sj,1) (s−sj)¸, s ≤sj.(15) −117 −116 −115 −114 −113 −112 −0.5 0 0.5 zoom de k(s) −250 −200 −150 −100 −50 0 50 100 −2 −1 0 1 2k(s) Figura 4: Periodizaci´on de k(s). En la figura 4 vemos c´omo hemos acoplado k(s, 1) con la exponencial adecuada. Gr´aficamente vemos que en la intersecci´on ambas curvas tienen la misma pendiente, pero la exponencial decae r´apidamente a cero. De ahora en adelante, nos referiremos a k(s, 1) como la uni´on de ambas curvas. Asimismo, designaremos su punto de uni´on como sjoint. Como estamos considerando k(s, 1) peri´odica, cuando realmente no lo es, a˜nadimos ceros en ambos extremos (gr´afica inferior), de tal forma que la informaci´on tarde en llegar a los bordes. Al final del proceso, tendremos N+ 1 puntos, siendo Nuna potencia de dos, para poder aplicar eficientemente la FFT. Asimismo, redefinimos sa≡s0ysb≡sN. Por ´ultimo, regularizamos kaplicando un filtro espectral suave ˆ k(ξ) = ˆ k(ξ) exp "−10 µ2.5|ξ| N¶25#,(16) que regulariza k(s, 1) en un entorno del punto de empalme sjoint, siendo inocuo lejos de dicho entorno. El haber obtenido una k(s, 1) regular y peri´odica permite integrar θ(s, 1) con exactitud espectral, mediante la FFT. Como tenemos un grado de libertad, fijaremos θ+= lims→∞ θ(s, 1) = θ(sb,1) ≡0. En la figura 5, observamos que θ(s, 1), asint´oticamente, no tiende al valor correcto cuando s→ −∞. Sin embargo, podemos determinar con bastante exactitud lims→−∞ θ(s, 1) como la media del primer m´aximo y el primer m´ınimo de θ(s, 1), con s > sjoint, cantidad que en la gr´afica superior pintamos con trazado m´as fino. Esto da una buena aproximaci´on de (4). Por ´ultimo, obrando exactamente como en el caso de la curvatura, forzamos que θ(s, 1) tienda exponencialmente al θ−predicho, volviendo a filtrar θ(s, t) al final del proceso. 6 An´alisis num´erico de soluciones de un flujo dispersivo −250 −200 −150 −100 −50 0 50 100 −6 −4 −2 0 2θ(s, 1) −122 −120 −118 −116 −114 −112 −110 −108 −106 −3.1 −3 −2.9 −2.8 −2.7 zoom de θ(s, 1) Figura 5: θ(s, t). 3.1.1. Experimentos num´ericos El objetivo del m´etodo que hemos desarrollado es intentar describir la formaci´on de una singularidad en tiempo finito de un flujo regular de curvas. En efecto, en t= 0, (1) desarrolla una singularidad en forma de esquina [2] z(s, 0) = (z0+seiθ+, s ≥0, z0+seiθ− , s ≤0.(17) Aqu´ı nos estamos limitando a s∈[sa, sb] y hemos fijado el ´angulo de la curva en saysb, por lo que impedimos que entre energ´ıa, que se define como Zsb sa k2(s0, t)ds0. En la soluci´on exacta, dicha energ´ıa crece conforme nos acercamos a t= 0, concentr´andose en un entorno de s= 0. Por el contrario, las simulaciones muestran que con nuestro m´etodo la energ´ıa se conserva con varios decimales de exactitud. Sin embargo, el m´etodo es v´alido cualitativamente, ya que la energ´ıa finita tiende tambi´en a concentrarse en un entorno de s= 0, aproxim´andose k(s, t) a una delta de Dirac, cuando t→0+. Asimismo, hemos observado que esta aproximaci´on mejora al aumentar el tama˜no del soporte de k(s, 1). −6 −5 −4 −3 −2 −1 0 −6 −5 −4 −3 −2 −1 0z(s, t), con t = 1, 0.5, 0.1, 0.01, 0.001 Figura 6: z(s, t). Una vez obtenido θ(s, t), es inmediato recuperar z(s, t), mediante zs(s, t) = exp(iθ(s, t)) y (9). En la figura 6, hemos superpuesto las gr´aficas de z(s, t) en los instantes de tiempo 7 Francisco de la Hoz t= 1, t= 0.5, t= 0.1, t= 0.01 y t= 0.001, en un entorno de s= 0. Se observa claramente el car´acter autosemejante de las soluciones. 4. Rk(s0, t)dt en funci´on de u(0) yux(0) Hemos integrado (10) con un conjunto bastante grande de pares admisibles (u(0), ux(0)) tomados como condiciones iniciales; a partir de cada soluci´on u(x) hemos definido su k(s, 1) correspondiente y hemos estimado R∞ −∞ k(s0,1)ds0. En la figura 7, pintamos cada par admisible (u(0), ux(0)) junto con el valor de la integral de k(s, 1) correspondiente; la curva as´ı pintada es antisim´etrica respecto del origen. Puesto que max µZ∞ −∞ k(s0,1, u(0), u0(0))ds0−π¶<3·10−3,(18) tenemos que −π−ε < Z∞ −∞ k(s0,1, u(0), u0(0))ds0< π +ε, 0< ε << 1.(19) Por tanto, hay a nuestro juicio, evidencia num´erica seria para aceptar que Z∞ −∞ u(x)dx ∈³−π 2,π 2´⇐⇒ Z∞ −∞ k(s0, t)ds0∈(−π, π),∀t. (20) −5 05 −10 0 10 −4 −3 −2 −1 0 1 2 3 4 −4 −2 0 2 4 −4 −3 −2 −1 0 1 2 3 4 Figura 7: Integral de la curvatura. Para valores de la integral entre [−3−ε, 3 + ε], es evidente por la gr´afica derecha, que representa la primera y tercera componentes, que s´olo hay una condici´on inicial de (10) adecuada. Los valores de la integral parecen converger de forma exponencial a π. Nos inclinamos a pensar que hay una ´unica condici´on inicial para cada valor (−π, π), si bien, en caso de tener convergencia exponencial, ser´a mucho m´as dif´ıcil dar evidencia num´erica. Referencias [1] R. E. Goldstein, D. M. Petrich, Soliton’s, Euler’s Equations, and Vortex Patch Dynamics, Phys. Rev. Letters. 69, 4 (1992) 555-558. [2] G. Perelman, L. Vega, Self-similar planar curves related to modified Korteweg-de Vries equation, To appear in J. of Diff. Eqns. 8