Aproximación numérica de quinto orden de las ecuaciones de Hamilton-Jacobi

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) Aproximaci´on num´erica de quinto orden de las ecuaciones de Hamilton-Jacobi S. Serna1 1Department of Mathematics, University of California Los Angeles, Los Angeles, CA 90095-1555, USA. E-mail: [email protected]la.edu. Palabras clave: Hamilton-Jacobi, Weighted Power-ENO, m´etodo de los “level sets” Resumen En este trabajo aproximamos la soluci´on de viscosidad de las ecuaciones de Hamilton- Jacobi asociadas al problema de la reinicializaci´on de curvas de nivel. Utilizamos para ello un m´etodo de quinto orden de precisi´on espacial ´optimo para la aproximaci´on de las ecuaciones de Hamilton-Jacobi. Como aplicaci´on calculamos la aproximaci´on a alto orden de funciones distancia euclidea signadas de curvas en R2. 1. Introducci´on Consideramos la aproximaci´on num´erica de las soluciones de viscosidad del problema de valores iniciales para las ecuaciones de Hamilton-Jacobi de la forma φt+H(x, φ, ∇φ) = 0, φ(x, 0) = φ0(x)x∈R t > 0 (1) donde Hes una funci´on de φno decreciente, ([2, 3]). Las ecuaciones de Hamilton-Jacobi aparecen en numerosas aplicaciones que abarcan desde el control ´optimo hasta la din´amica de fluidos. El estudio de m´etodos num´ericos eficientes y precisos para resolverlas ha sido un campo de investigaci´on muy activo en los ´ultimos veinte a˜nos, ([6, 7, 10, 13, 14]). La soluci´on de viscosidad de la ecuaci´on de Hamilton-Jacobi desarrolla discontinuidades en derivada en tiempo finito. Los esquemas num´ericos en general y de primer orden en particular tienden a suavizar estas singularidades. Por esto es necesario dise˜nar m´etodos num´ericos de alto orden en espacio para aproximar las singularidades de la soluci´on de forma precisa ([6, 10, 13]). La aproximaci´on num´erica del movimiento de interfases gobernado por un sistema de ecuaciones en derivadas parciales es un problema importante en f´ısica computacional 1 S. Serna debido a la acumulaci´on de ruido y difusi´on introducidos por el m´etodo num´erico y a la dificultad de reproducir los cambios topol´ogicos que sufre la interfase a lo largo de la evoluci´on. Este fen´omeno aparece en la din´amica de fluidos polif´asicos compresibles e incompresibles [1, 5, 8]. La descripci´on de interfases en movimiento mediante la curva de nivel cero de una funci´on impl´ıcita que evoluciona en tiempo con la din´amica es uno de los fundamentos del m´etodo de los conjuntos de nivel, ([4, 8]). El m´etodo de los conjuntos de nivel establece una clase importante de ecuaciones de Hamilton-Jacobi donde el campo de velocidades que mueve las curvas de nivel es funci´on de la curvatura media [9]. La funci´on distancia eucl´ıdea signada a la interfase es la funci´on impl´ıcita que mejor describe la curva de nivel cero. Una aproximaci´on de alta precisi´on de esta funci´on distancia permite realizar c´alculo de integrales precisas sobre la interfase. El concepto de reinicializaci´on fue introducido en [15] para mantener la precisi´on de la interfase en movimiento descrita como la curva de nivel cero de una funci´on impl´ıcita y reducir su degradaci´on cuando evoluciona en tiempo. La ecuaci´on de reinicializaci´on de los conjuntos de nivel propuesta en [15] es una ecuaci´on de Hamilton-Jacobi que permite aproximar la funci´on distancia eucl´ıdea signada. En este trabajo proponemos calcular con alta precisi´on espacial mediante la ecuaci´on de reinicializaci´on, funciones distancia de curvas mediante el esquema num´erico de alto orden Weighted PowerENO propuesto en [13] para ecuaciones de Hamilton-Jacobi de evoluci´on. En la secci´on 2 se describe el esquema num´erico y en la secci´on 3 se presenta un conjunto de experimentos donde se calcula de manera aproximada la funci´on distancia signada de algunas curvas en R2. 2. Esquema num´erico En esta secci´on presentamos un esquema en diferencias finitas de alto orden para las ecuaciones de Hamilton-Jacobi. Describimos el Hamiltoniano num´erico mon´otono debido a Osher y Sethian ([9]) y un m´etodo de reconstrucci´on de quinto orden ´optimo para la aproximaci´on de la soluci´on de las ecuaciones de Hamilton-Jacobi propuesto en [13]. La estrategia consiste en primero aproximar con diferencias finitas las derivadas parciales de la soluci´on a alto orden y despues insertarlas en el Hamiltoniano num´erico. La integraci´on en tiempo se realiza mediante un m´etodo de alto orden Runge-Kuta ([16]). 2.1. Un Hamiltoniano num´erico mon´otono Describimos el esquema num´erico b´asico de primer orden para el problema de valores iniciales de la ecuaci´on de Hamilton-Jacobi en dos dimensiones espaciales φt+H(φx, φy) = 0, φ(x, y, 0) = φ0(x, y), t > 0 (2) mediante una malla computacional definida por los incrementos ∆x, ∆yy ∆t, siendo φn j,k la aproximaci´on num´erica a la soluci´on de viscosidad de (2) φ(xj, xk, tn) = φ(j∆x, k∆y, n∆t) (3) 2 Aproximaci´on num´erica de las ecuaciones de Hamilton-Jacobi Eliminamos la dependencia de x, y yφpara simplicidad en la exposici´on. Consideramos el esquema de primer orden, φn+1 j,k =φn j,k −∆t g(Dx −φj,k, Dx +φj,k, Dy −φj,k, Dy +φj,k) (4) donde ges un flujo Lipschitz continuo y mon´otono en cada una de las variables y consistente con el Hamiltoniano de la ecuaci´on H,g(u, u, v, v) = H(u, v) y Dx −φj,k =φj,k −φj−1,k ∆x;Dx +φj,k =φj+1,k −φj,k ∆x; (5) Dy −φj,k =φj,k −φj,k−1 ∆y;Dy +φj,k =φj,k+1 −φj,k ∆y.(6) Una funci´on gcon estas propiedades se denomina Hamiltoniano num´erico mon´otono asociado a la ecuaci´on original. En este trabajo aproximamos la soluci´on de ecuaciones de Hamilton-Jacobi con un Hamiltoniano H(u, v) que cumpla H(u, v) = h(u2, v2) tal que h1·h2>0, donde hjes la derivada parcial de hrespecto del argumento jde la funci´on h. Para este tipo de Hamiltonianos se puede utilizar el Hamiltoniano num´erico de Osher-Sethian [9] definido como: g(u−, u+, v−, v+) =    h([max((u−)+,(u+)−)]2,[max((v−)+,(v+)−)]2)h1≥0 h([max((u−)−,(u+)+)]2,[max((v−)−,(v+)+)]2) otherwise donde (a)+= max(a, 0) y (a)−=−min(a, 0). Este Hamiltoniano num´erico es mon´otono y define un esquema num´erico de primer orden que converge a la soluci´on de viscosidad, ([9]). 2.2. Reconstrucci´on espacial de quinto orden para las ecuaciones de Hamilton-Jacobi En esta secci´on revisamos el m´etodo de reconstrucci´on Weighted Power-ENO de quinto orden de precisi´on espacial. Los m´etodos Weighted Power-ENO ([12]) se propusieron originalmente para leyes de conservaci´on hiperb´olicas. La idea fundamental de las reconstrucciones Power-ENO consiste en aplicar una clase extendida de limitadores ([11, 12]) a diferencias de segundo orden de los m´etodos de reconstrucci´on ENO cl´asicos de manera que la reconstrucci´on retiene m´as informaci´on de las escalas finas y mejora la resoluci´on alrededor de las discontinuidades de la soluci´on. Una estrategia de ponderaci´on basada en indicadores de suavidad apropiados ([12]) se usa para conseguir una aproximaci´on de quinto orden. El m´etodo de quinto orden que resulta es el llamado m´etodo Weighted Power-ENO. En [13] se adaptan los m´etodos de reconstrucci´on Weighted PowerENO para abordar la aproximaci´on de la soluci´on de viscosidad de las ecuaciones de las Hamilton-Jacobi y as´ı mejorar la resoluci´on alrededor de las esquinas donde el gradiente de la soluci´on es discontinua. 3 S. Serna Los limitadores powerp, a partir de los cuales se dise˜nan estos m´etodos, se definen a partir de las siguientes medias: dados dos n´umeros positivos x > 0 e y > 0, para un n´umero natural p, en [12] se define la media powerp(x, y) como: powerp(x, y) = (x+y) 21− x−y x+y p(7) Las siguientes desigualdades se satisfacen para todo x > 0 e y > 0: min(x, y)≤powerp(x, y)≤powerq(x, y)≤x+y 2 para 0 < p < q. Obs´ervese que l´ımp→∞ powerp(x, y) = x+y 2:= power∞(x, y). Dados los valores puntuales φ(xj), j = 1,2,··· de una funci´on en nodos equiespaciados xj, tal que xj+1 =xj+ ∆xdenotamos sus diferencias finitas divididas como zj+1 2=φj+1 −φj ∆x(8) Introducimos la siguiente notaci´on para las diferencias: dj=zj+1 2−zj−1 2(9) dj+1 2=dj+dj+1 2(10) Dj+1 2=dj+1 −dj(11) La reconstrucci´on Weighted PowerENO est´a basada en una combinaci´on convexa de las tres par´abolas siguientes asociadas a cada intervalo Ij= [xj, xj+1] pP j(x) = zj+1 2−Pj 24 + x−xj+1 2 ∆xdj+Pj 2+Pj 2x−xj+1 2 ∆x (12) pj+1 2(x) = zj+1 2− Dj+1 2 24 + x−xj+1 2 ∆x"dj+1 2+ Dj+1 2 2x−xj+1 2 ∆x#(13) pP j+1(x) = zj+1 2−Pj+1 24 + x−xj+1 2 ∆xdj+1 −Pj+1 2+Pj+1 2x−xj+1 2 ∆x (14) donde Pj:= powermodp(Dj−1 2 , Dj+1 2) y powermodp(x, y) = sign(x) + sign(y) 2powerp(|x|,|y|) (15) En particular, en x=xjtenemos 4 Aproximaci´on num´erica de las ecuaciones de Hamilton-Jacobi pP j(xj) = zj+1 2−1 2dj−1 6Pj(16) pj+1 2(xj) = zj+1 2−1 2dj+1 2−1 6Dj+1 2(17) pP j+1(xj) = zj+1 2−1 2dj+1 +1 3Pj+1 (18) Para obtener la precisi´on ´optima de u+(xj) en la interfase de Ijusamos las siguiente combinaci´on convexa u+(xj) = w0·pP j(xj) + w1·pj+1 2(xj) + w2·pP j+1(xj) (19) donde wk=αk α0+α1+α2 k= 0,1,2 (20) y αk=Ck (ǫ+ISk)2(21) aqu´ı C0= 0,6, C1= 0,2 y C2= 0,2 son los pesos ´optimos. Para calcular los indicadores de suavidad utilizamos la norma L2de las derivadas de los polinomios de manera que se alcanza el grado ´optimo de precisi´on. Obtenemos las expresiones siguientes: IS0=13 12 (Pj)2+1 42zj+1 2−2zj−1 2+Pj2 (22) IS1=13 12 zj−1 2−2zj+1 2+zj+3 22 +1 4zj−1 2−zj+3 22 (23) IS2=13 12 (Pj+1)2+1 42zj+3 2−2zj+1 2−Pj+12 (24) Para u−(xj) se obtiene una expresi´on similar para el polinomio asociado en Ij−1, u−(xj) = w0·pP j−1(xj) + w1·pj−1 2(xj) + w2·pP j(xj) (25) con C0= 0,2, C1= 0,2 y C2= 0,6 como pesos ´optimos. Una expresi´on reducida de estas par´abolas evaluadas en x=xjes, pP j−1(xj) = zj−1 2+1 2dj−1+1 3Pj−1(26) pj−1 2(xj) = zj−1 2+1 2dj−1 2+1 12Dj−1 2(27) pP j(xj) = zj−1 2+1 2dj−1 6Pj(28) 5 S. Serna El m´etodo que resulta es un m´etodo Weighted Power ENO de quinto orden de precisi´on para p≥3 como se demuestra en [12]. En [12] se demuestra que el valor ´optimo de ppara la aproximaci´on de leyes de conservaci´on hiperb´olicas donde se desarrollan discontinuidades en tiempo finito de manera que el esquema sea de variaci´on total local acotada ha de ser p= 3. En [13] se estudia aplicar los m´etodos Weighted PowerENO para la aproximaci´on de las soluciones de viscosidad de Hamilton-Jacobi utilizando limitadores powerpm´as d´ebiles con p > 3. En particular se propone la media aritm´etica (l´ımite infinito de las medias powerp) en lugar de los limitadores power3utilizados para leyes de conservaci´on. El m´etodo, al que se denomina Weighted Power∞ENO, es ´util para la aproximaci´on de soluciones de viscosidad de las ecuaciones de Hamilton-Jacobi ([13]) ya que estas soluciones son continuas con derivadas discontinuas y el procedimiento de reconstrucci´on se aplica a diferencias divididas de primer orden. 3. Reinicializaci´on de funciones de conjuntos de nivel Para ciertas curvas en R2no existe expresi´on anal´ıtica expl´ıcita de su funci´on distancia. En esta secci´on aproximamos con alta precisi´on la funci´on distancia euclidea signada de diferentes curvas sencillas. El c´alculo de la funci´on distancia se realiza aproximando la soluci´on de viscosidad de la ecuaci´on de reinicializaci´on tomando como dato inicial una perturbaci´on de dicha funci´on. La ecuaci´on de reinicializaci´on de los conjuntos de nivel es de la forma φt+ sign(φ0)hqφ2 x+φ2 y−1i= 0, φ(x, y, 0) = φ0(x, y) (29) cuyo dato inicial es una funci´on φ0(x, y) que tiene como curva de nivel cero la curva de R2para la que deseamos calcular la funci´on distancia. La funcion sign(φ0) representa la regularizaci´on de la funci´on signo sign(φ0) = φ0 pφ2 0+ǫ, ǫ =O(∆x) (30) El problema de valores iniciales (29) tiene como soluci´on estacionaria la funci´on distancia signada a la curva dada. Resolvemos (29) usando el Hamiltoniano num´erico de Osher-Sethian a primer y quinto orden (Weighted Power∞ENO ) para diferentes datos iniciales. Comparamos la precisi´on de los m´etodos calculando los errores globales en norma L∞con diferentes mallas computacionales. Los resultados se muestran en tablas donde se aprecia una mejora considerable en la reducci´on del error para los c´alculos realizados a quinto orden. 3.1. Experimento 1. Funci´on distancia regular Elegimos la funci´on φ0como la funci´on distancia signada a la circunferencia de centro el origen y radio 1 2perturbada en las direcciones radial y angular: φ0(x, y) = d+δ|d| ≤ ǫ den otro caso 6 Aproximaci´on num´erica de las ecuaciones de Hamilton-Jacobi N Primer orden Weighted Power∞ 100 3.85 10−32.5 10−5 200 2.43 10−34.51 10−7 Tabla 1: Errores globales en norma L∞para la funci´on distancia de la circunferencia N Primer orden Weighted Power∞ 100 8.10 10−12.17 10−2 200 1.03 10−17.61 10−3 Tabla 2: Errores globales en norma L∞para la funci´on distancia de la lemniscata donde d=px2+y2−0,5, θ= tan−1(y |x|), ǫ= 0,2 y δ=ǫ 16πsin 4πd sin 5θ ǫ Realizamos el c´alculo en el dominio [−1,1] ×[−1,1] utilizando las mallas computacionales 100 ×100, 200×200 con una restricci´on CFL de 0,6 y 256 y 1024 iteraciones respectivamente. Los errores en norma L∞para cada uno de los m´etodos se presentan en la Tabla 1. 3.2. Experimento 2. Funci´on distancia con discontinuidades en derivada Consideramos la restauraci´on de una funci´on distancia no regular. Elegimos φ0como la funcion distancia a la lemniscata a4= [(x−a)2+y2][(x+a)2+y2] (31) con a= 0,5 y a˜nadiendo una perturbaci´on en direcci´on radial y angular φ0(x, y) = d+ 3δ|d| ≤ ǫ den otro caso definida en [−1,1] ×[−1,1] donde θ= tan−1(y |x|), ǫ= 0,2 y δ=ǫ 16πsin 4πd sin 5θ ǫ. Calculamos la funci´on distancia dde la lemniscata aproximando la evoluci´on de la ecuaci´on (29) al estado estacionario tomando como dato incial la funci´on D(x, y) = p[(x−a)2+y2][(x+a)2+y2]−a2(32) Realizamos el c´alculo utilizando las mallas computacionales 100×100, 200×200 con una restricci´on CFL de 0,6 y 256 y 1024 iteraciones respectivamente. Se toma como referencia la soluci´on para la malla computacional 400 ×400. Los errores en norma L∞para cada uno de los m´etodos se presentan en la Tabla 2. 3.3. Experimento 3. Curva no regular con interior convexo Consideramos la curva frontera de la bola de centro cero de radio 1 2en la norma L1, {(x, y) : |x|+|y|=1 2}. Como dato inicial elegimos la funci´on distancia signada en la norma L1 φ0(x, y) = |x|+|y| − 1 2(33) 7 S. Serna N Primer orden Weighted Power∞ 100 2.03 10−21.35 10−3 200 8.12 10−34.50 10−4 Tabla 3: Errores globales en norma L∞para la funci´on distancia de la bola de radio 1 2en norma L1 Realizamos el c´alculo en el dominio [−1,1] ×[−1,1] utilizando las mallas computacionales 100 ×100, 200×200 con una restricci´on CFL de 0,6 y 256 y 1024 iteraciones respectivamente. Los errores en norma L∞para cada uno de los m´etodos se presentan en la Tabla 3. Los errores globales para las soluciones aproximadas de la lemniscata y la frontera de la bola en norma L1son significativamente mayores que los correspondientes a la circunferencia debido a la presencia de discontinuidades en derivada. Referencias [1] I. -L. Chern, J. Glimm, O. Mcbryan, B. Plohr and S. Yaniv , Front tracking for gas dynamics J. Comput. Phys., 62,(1986), pp. 83–110. [2] Crandall, M. and Lions, P., Viscosity solutions of Hamilton-Jacobi equations, Trans. Amer. Math. Soc., 277, (1983), pp. 1–42. [3] Crandall, M. and Lions, P., Two approximations of solutions of Hamilton-Jacobi equations, Math. Comput, 43, (1984), pp. 1–19. [4] Fedkiw, R, Aslam, T, Merriman, B. and Osher, S. , A Non-oscillatory Eulerian Approach to Interfaces in Multimaterial Flows (the Ghost Fluid Method) J. Comput. Phys., 152 ,(1999), pp. 457–492. [5] Hirt, C.W. and Nichols B.D. , Volume of fluid (VOF) method for the dynamics of free boundaries J. Comput. Phys., 39,(1981), pp. 201–225. [6] Jiang, G.S., and Peng, D.,Weighted ENO schemes for Hamilton-Jacobi equations, SIAM J. Sci. Comput., 21,(2000), pp. 2126–2143. [7] Lin, C. T. and Tadmor, E., High-resolution nonoscillatory central schemes for Hamilton-Jacobi Equations, SIAM J. Sci. Comput., 21, (2000), pp 2163-2186. [8] Mulder, W., Osher, S.J. and Sethian, J. , Computing interface motion in compressible gas dynamics J. Comput. Phys., 100, (1992), pp. 209-228. [9] Osher, S.J. and Sethian, J. , Fronts propagating with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulations, J. Comput. Phys., 79, (1988), pp. 12–49. [10] Osher, S.J. and Shu, C.W. , High-order essentially non-oscillatory schemes for Hamilton-Jacobi equations, SIAM J. Numer. Anal., 28, (1991), pp. 907–922. [11] Serna, S., A class of extended limiters applied to Piecewise Hyperbolic methods, SIAM Journal of Scientific Computing, v28 (1), 123-140, 2006. [12] Serna, S. and Marquina, A., Power ENO Methods: A fifth order accurate Weighted Power ENO method, J. Comput. Phys., 194, (2004), pp. 632–658. [13] Serna, S., Qian, J., Fifth Order Weighted Power-ENO schemes for Hamilton-Jacobi equations, Journal of Scientific Computing, v29, 57-81, 2006. [14] Souganidis, P. E., Approximation schemes for viscosity solutions of Hamilton-Jacobi equations, J. Diff. Eqn., 59, (1985), pp 1–43. [15] Sussman, M. Smereka, P. and Osher, S., A Level Set Approach for Computing Solutions to Incompressible Two-Phase Flow J. Comput. Phys., 114, (1994), pp. 146–159. [16] Spiteri, R. J. and Ruuth, S. J., A new class of optimal high order strong stability preserving time discretization methods, SIAM J. Num. Anal., 40, (2002), pp 469-491. 8