Métodos linealmente implícitos de tipo Runge-Kutta de Pasos Fraccionarios aplicados a problemas parabólicos semilineales: reducción de orden y técnicas para evitarla

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) M´etodos linealmente impl´ıcitos de tipo Runge-Kutta de Pasos Fraccionarios aplicados a problemas parab´olicos semilineales: reducci´on de orden y t´ecnicas para evitarla B. Bujanda1, J.C. Jorge 2, M. J. Moreta3 1,2Dpto. Ingenier´ıa Matem´atica e Inform´atica, Campus Arrosad´ıa s/n, Universidad P´ublica de Navarra, Pamplona (Navarra). E-mails: [email protected], [email protected]. 3Dpto. de Fundamentos del An´alisis Econ´omico I. Universidad Complutense de Madrid. E-mail: [email protected] Palabras clave: M´etodos Runge-Kutta de Pasos Fraccionarios, Problemas semilineales, Reducci´on de orden Resumen A linearly implicit variant of fractional-step Runge-Kutta methods was introduced in [2] in order to integrate efficiently semi-linear multidimensional parabolic problems. The computational cost reduction of these methods is remarkable, compared to classical implicit methods but they present the classical drawback of order reduction which has its maximum relevance when they are applied to problems with boundary conditions depending on the time. In this paper we show a technique which permits us to avoid this reduction in some cases. 1. Introducci´on En la naturaleza existen diversos fen´omenos que son modelados mediante ecuaciones diferenciales parab´olicas semilineales; por ejemplo, en [7] y sus referencias podemos encontrar modelos de poluci´on de aire descritos en la forma: ∂c ∂t +∇ · (u c) = ∇ · (K∇c) + R(c) donde ces un vector de concentraciones de especies qu´ımicas, ues el vector velocidad, K es la matriz de difusi´on y Runa funci´on no lineal. En esta comunicaci´on vamos a mostrar m´etodos num´ericos que integran eficientemente problemas parab´olicos semilineales que, en forma abstracta, admiten ser enunciados como sigue: 1 B. Bujanda, J.C. Jorge, M.J. Moreta Encontrar y: [t0, T ]→ H soluci´on de    y0(t) = L y(t) + f(t) + g(t, y),∀t∈[t0, T ], y(t0) = y0, B y(t) = β(t)∈ Hb,∀t∈[t0, T], (1) siendo HyHbdos espacios de Hilbert de funciones definidas sobre Ω ⊆Rsy sobre su frontera Γ respectivamente; Bun operador frontera definido entre HyHb,Lun operador lineal generalmente no acotado, definido sobre dominio denso en Hcon valores en H,f(t) el t´ermino fuente, g(t, y) una funci´on regular que marca el aporte no lineal y β(t) una condici´on de contorno dependiente del tiempo. Consideraremos que se tiene suficiente regularidad y compatibilidad en los datos para asegurar que la soluci´on de (1) es suficientemente diferenciable. Supondremos adem´as que el operador −L0:D ⊆ H −→ H que es la restricci´on del operador Lal Ker(B) es lineal, en general no acotado, maximal y mon´otono, es decir, ∀u∈ H,∃v∈ D,tal que (I−L0)v=uy ((−L0v, v )) ≥0,1∀v∈ D respectivamente. La obtenci´on de soluciones num´ericas para problemas de la forma (1) normalmente requiere realizar un doble proceso de discretizaci´on espacial y temporal. Para la discretizaci´on temporal se puede optar por utilizar m´etodos expl´ıcitos o impl´ıcitos. Si se utilizan m´etodos expl´ıcitos, el costo computacional del m´etodo ser´a elevado en el caso de utilizar mallas finas para la variable espacial, ya que tendremos que utilizar pasos de tiempo peque˜nos para asegurar la estabilidad num´erica del algoritmo. Si nos decidimos por m´etodos impl´ıcitos cl´asicos nos encontraremos con el problema de tener que resolver en cada etapa sistemas no lineales que pueden llegar a tener una dimensi´on muy elevada; al resolver estos sistemas utilizando m´etodos iterativos se suele observar habitualmente una velocidad de convergencia lenta. Adem´as en el caso de problemas no lineales se han utilizado habitualmente m´etodos num´ericos con propiedades de tipo B-estabilidad. Este tipo de propiedades es muy restrictiva e implica la utilizaci´on de m´etodos de orden bajo o, en el caso de optar por m´etodos de orden alto, m´etodos totalmente impl´ıcitos. Con el fin de evitar estos inconvenientes, en [2] los autores presentaron un tipo de m´etodos que integra de manera muy eficiente esta clase de problemas. La base de los m´etodos indicados es la misma que la de los m´etodos de Pasos Fraccionarios cuando son aplicados a problemas parab´olicos lineales. As´ı los m´etodos se construyen combinando discretizaciones espaciales adecuadas con discretizaciones temporales de tipo Runge-Kutta Aditivo en las que la aportaci´on de la parte lineal de la funci´on derivada, L y(t) + f(t), se define mediante un m´etodo Runge-Kutta de Pasos Fraccionarios (RKPF) (ver [6]) y de la parte no lineal, g(t, y), mediante un m´etodo RK expl´ıcito adecuado. Los algoritmos que se obtienen siguiendo este proceso presentan la caracter´ıstica de ser incondicionalmente convergentes imponiendo solamente propiedades de tipo estabilidad absoluta lineal (ver [3]); adem´as, si la descomposici´on del operador el´ıptico y del t´ermino fuente se realiza de manera adecuada, el costo computacional es bajo, comparado con los algoritmos obtenidos con m´etodos impl´ıcitos cl´asicos, incluyendo algunos otros m´etodos linealmente impl´ıcitos, ya que, en este caso, en la resoluci´on de las etapas intermedias aparecen sistemas lineales reducibles a una serie de subsistemas (muy simples en algunos casos) cuya resoluci´on puede hacerse en paralelo. 1((·,·)) denota un producto escalar en H; posteriormente usaremos k·k para indicar su norma asociada. 2 Reducci´on de orden en m´etodos RKPF aplicados a problemas parab´olicos semilineales Al integrar este tipo de problemas mediante m´etodos de un paso suele presentarse el fen´omeno conocido como reducci´on de orden, que alcanza su m´axima expresi´on cuando las condiciones de contorno var´ıan en el tiempo (ver [4, 5]). Mostraremos una t´ecnica para evitar esta reducci´on de orden similar a la propuesta en [1]. Este proceso tiene la ventaja del bajo costo computacional adicional que a˜nade al algoritmo, puesto que s´olo afecta a la evaluaci´on de ciertos datos en la frontera del dominio Ω. 2. Un m´etodo de discretizaci´on temporal El m´etodo de integraci´on temporal linealmente impl´ıcito introducido en [2] proporciona aproximaciones num´ericas Ymde la soluci´on exacta y(tm) mediante el siguiente esquema:                    Y0=y0∈ H, Ym,i =Ym+τ i X j=1 akj ij LkjYm,j +fkj(tm,j)+τ i−1 X j=1 an+1 ij g(tm,j, Y m,j), para i= 1, . . . , s, Ym+1 =Ym+τ s X i=1 bki iLkiYm,i +fki(tm,i)+τ s X i=1 bn+1 ig(tm,i, Y m,i). (2) donde τindica el paso en tiempo (que, por simplicidad, consideraremos constante), tm= t0+m τ ytm,i =t0+ (m+ci)τ,∀m= 0,1,2, . . .,nes el n´umero de niveles del m´etodo RKPF involucrado, kj∈ {1, . . . , n}para j= 1,2, . . . , s y las aproximaciones intermedias Ym,i para i= 1, . . . , s son denominadas etapas intermedias del m´etodo, las cuales pueden ser consideradas aproximaciones num´ericas de la soluci´on exacta y(t) en tm,i. Para poder aplicar el esquema anterior, en el caso de considerar condiciones de contorno homog´eneas, es suficiente asumir que el operador lineal Ladmite una descomposici´on en nsumandos de la forma L=Pn i=1 Li(t) con Li:Di→ H maximales, mon´otonos y tales que D=Tn i=1 Di. Definimos los operadores Bi:Di→ Hb i,i= 1, . . . , n tales que Ker(B) = ∩n i=1Ker(Bi) y Ker(Bi) = Di(ver [2, 3]). Para el caso general, es decir, condiciones de contorno no homog´eneas, supondremos adem´as que los operadores −Li,0: Di⊆ H −→ H que son la restricci´on de los operadores −Lial Ker(Bi) son lineales, en general no acotados, maximales y mon´otonos. Para que la expresi´on de los algoritmos, as´ı como de los resultados, sea lo m´as sim´etrica posible, se ha descompuesto adem´as el t´ermino fuente en nsumandos suficientemente regulares f(t) = Pn i=1 fi(t). Los coeficientes de estos m´etodos se pueden organizar en una tabla de la forma: C A1. . . AnAn+1 (b1)T. . . (bn)T(bn+1)T donde C=diag(c1, . . . , cs)∈Rs×s. Para obtener estas tablas hemos completado las ma3 B. Bujanda, J.C. Jorge, M.J. Moreta trices con algunos coeficientes nulos en la forma: Ak= (ak ij)∈Rs×sdonde ak ij =     an+1 ij si k=n+ 1y i > j, akj ij si k=kjyi≥j, 0 otro caso, bk= (bk j)∈Rsdonde bk j=     bn+1 jsi k=n+ 1, bkj jsi k=kj, 0 otro caso. En (2) puede observarse como el costo computacional de estos esquemas puede llegar a ser muy bajo, ya que en cada etapa tenemos que resolver problemas en los que aparece impl´ıcitamente tan solo uno de los operadores en los que se ha fraccionado el operador el´ıptico L, y el aporte de la parte no lineal g(t, y) viene dado de forma expl´ıcita. As´ı, los sistemas que se deben resolver en cada una de las etapas son lineales. No obstante para que cada una de estas etapas tenga una ´unica soluci´on en el caso general, es decir, condiciones de contorno dependientes del tiempo, es necesario a˜nadir en el esquema las condiciones frontera de las etapas. As´ı,para cada i= 1, . . . , s tenemos que resolver: (I−τ aki ii Lki)Ym,i =Ym+τ i−1 X j=1 akj ij LkjYm,j +fkj(tm,j)+τ i−1 X j=1 an+1 ij g(tm,j, Y m,j), BkiYm,i =βm,i donde se suele tomar como valor de contorno βm,i ≡Bkiu(tm,i) = β(tm,i).(3) Asumiremos que los problemas de contorno de la forma (I−d lki)u=f, Bkiu=β, admiten soluci´on ´unica para todo d∈(0, D] cumpliendo adem´as que kuk ≤ C(kfk+ kβkb), siendo Cuna constante independiente de d. Las condiciones de contorno elegidas en la forma (3) provocan una reducci´on en el orden del m´etodo, tal y como veremos posteriormente. El esquema (2), con las condiciones de contorno indicadas para las etapas puede ser escrito en forma abreviada utilizando la siguiente notaci´on tensorial: IHdenota la matriz identidad en H; dados M≡(mij)∈Rs×syv≡(vi)∈Rs, denotamos por ¯ M≡(mijIH)∈ L(Hs,Hs) y ¯v≡(viIH)∈ L(Hs,H). Agrupamos las etapas Ym,i, las evaluaciones de fi(t), g(t, y) y Lipara todo i= 1,· · · , n, y para todo m= 1,2,· · · , en la forma Ym=Ym,1,· · · , Y m,sT∈ Hs, Fm i= (fi(tm,1),· · · , fi(tm,s))T∈ Hs, Gm(Ym) = g(tm,1, Y m,1),· · · , g(tm,s, Y m,s)T∈ Hs, ˆ Li=diag (Li,· · · , Li)∈ L(Ds i,H), B Ym= (Bk1Ym,1,· · · , BksYm,s)T∈ Hb k1× · · · × Hb ks, βm= (βm,1,· · · , βm,s)T∈ Hb k1× · · · × Hb ks. 4 Reducci´on de orden en m´etodos RKPF aplicados a problemas parab´olicos semilineales As´ı, se tiene que las etapas vienen dadas por     ¯ I−τ n X i=1 Aiˆ LiYm= ¯e Y m+τ n X i=1 AiFm i+τAn+1 Gm(Ym), B Ym=βmsiendo βm,i =Bkiy(tm,i), y para avanzar un paso aplicamos: Ym+1 =Ym+τ n X i=1 biTˆ LiYm+Fm i+τbn+1TGm(Ym). Recordemos que un m´etodo RKPF se dice consistente con orden cl´asico psi el error local verifica, para funciones y(t) suficientemente regulares, que kem+1k=ky(tm+1)− ¯ Ym+1k=O(τp+1) donde ¯ Ym+1 es la soluci´on num´erica obtenida a partir del esquema:                    ¯ I−τ n X i=1 Aiˆ Li¯ Ym= ¯e y(tm) + τ n X i=1 AiFm i+τAn+1 Gm(¯ Ym), B¯ Ym=βm,siendo βm,i =Bkiy(tm,i), ¯ Ym+1 =y(tm) + τ n X i=1 biTˆ Li¯ Ym+Fm i+τbn+1TGm(¯ Ym). El orden observado finalmente en el esquema tiene que ver con el orden de las etapas intermedias; un m´etodo RKPF se dice que tiene orden de etapas qcuando q= m´ın{p, ˜q} siendo pel orden cl´asico y ˜qel m´aximo valor tal que, (C)ke−kAj(C)k−1e= 0, ∀j= 1, . . . , n+1, k= 1, . . . , ˜q. Notar que, en este tipo de m´etodos si se considera que la primera etapa no es expl´ıcita, el orden por etapas es cero, ya que An+1 e6=Akeen la primera componente para alg´un k= 1, . . . , n; por ello recuperar el orden perdido es especialmente importante en estos casos. Para ver que se produce una p´erdida de orden descomponemos el error local en la forma2: em+1,[0] =y(tm+1)−¯ Ym+1,[0] = (y(tm+1)−e Ym+1,[0])+(e Ym+1,[0] −¯ Ym+1,[0]), donde hemos introducido e Ym+1,[0] que se obtiene a partir de e Ym+1,[0] =y(tm) + τ n X i=1 biTˆ Lie Ym,[0] +Fm i+τbn+1TGm(e Ym,[0]), e Ym,[0] = [Ym,1,[0], . . . , Y m,s,[0]]T,siendo Ym,i,[0] =y(tm,i). Esta soluci´on est´a relacionada con los errores que aparecen en las f´ormulas de cuadratura de las etapas intermedias que llamaremos, δm,[0], y que definimos como sigue     ¯ I−τ n X i=1 Aiˆ Lie Ym,[0] = ¯e y(tm) + τ n X i=1 AiFm i+τAn+1 Gm(e Ym,[0]) + δm,[0], Be Ym,[0] =βm,[0]. 2El super´ındice [0] indica que no se ha realizado ninguna modificaci´on, por ejemplo, em+1,[0] indica em+1. 5 B. Bujanda, J.C. Jorge, M.J. Moreta Para acotar el primero de los sumandos (y(tm+1)−e Ym+1,[0]) realizamos adecuados desarrollos de Taylor para las expresiones anteriores y aplicamos las condiciones de orden (bk1)T(C)ρe=1 ρ+ 1, 0 ≤ρ≤p−1, obteniendo ky(tm+1)−e Ym+1,[0]k=O(τp+1). Para acotar el segundo de los sumandos realizamos tambi´en en este caso adecuados desarrollos de Taylor para la expresi´on de δm,[0] y, teniendo en cuenta el orden de las etapas intermedias, se tiene: ke Ym+1,[0] −¯ Ym+1,[0]k=O(τm´ın{˜q+1,p+1}) = O(τq+1). De esta forma se obtiene que kem+1,[0]k ≤ O(τq+1) + O(τp+1) = O(τm´ın{q+1,p+1}). Siguiendo las ideas utilizadas en [1] para evitar la reducci´on de orden que puede aparecer en la integraci´on num´erica de problemas no lineales (los autores se refieren a m´etodos de tipo Rosenbrock) se definen los valores de contorno de las etapas intermedias en la forma: Ym,i,[0] =y(tm,i), βm,i,[0] =BkiYm,i,[0] =Bkiy(tm,i), i = 1, . . . , s,      Ym,i,[1] =y(tm) + τ i X j=1 akj ij LkjYm,kj,[0] +fkj(tm,j)+τ i−1 X j=1 an+1 ij g(tm,j, Y m,j,[0]), βm,i,[1] =BkiYm,i,[1], i = 1, . . . , s. Utilizando notaci´on tensorial queda:      Ym,[1] = ¯e y(tm) + τ n X i=1 AiYm,[0] +τ n X i=1 AiFm i+τAn+1 Gm(Ym,[0]), B Ym,[1] =βm,[1],siendo βm,i,[1] =BkiYm,i,[1]. En este caso, el error local mejorado se calcula en la forma: em+1,[1] =y(tm+1)−¯ Ym+1,[1] = (y(tm+1)−e Ym+1,[1])+(e Ym+1,[1] −¯ Ym+1,[1]), expresi´on que nos permite deducir que kem+1,[1]k ≤ O(τq+2) + O(τp+1) = O(τm´ın{q+2,p+1}), y por lo tanto afirmar que de forma muy sencilla, se puede recuperar un orden siempre que haya una p´erdida de orden debida a que el orden de las etapas es inferior al orden del m´etodo, lo que es habitual. 3. Discretizaci´on total Para obtener la soluci´on num´erica del problema inicial debemos completar el proceso de discretizaci´on temporal con una segunda etapa de discretizaci´on en espacio. La formulaci´on que vamos a indicar permite incluir tanto el caso de Diferencias Finitas como el de Elementos Finitos. Para discretizar el dominio espacial Ω consideramos un par´ametro h∈(0, h0] destinado a tender a cero y para cada htomamos un espacio de Hilbert Vh 6 Reducci´on de orden en m´etodos RKPF aplicados a problemas parab´olicos semilineales finito dimensional con norma k · kh; para la discretizaci´on de la frontera consideramos Vb h, Vb i,h y las normas k·kb hyk·kb i,h. Para conectar los espacios consideramos las siguientes aplicaciones: rh:D ⊂ H −→ Vhcumpliendo que l´ım h→0krhykh=kyk,∀y∈ D, ri,h :Di⊂ H −→ Vhcumpliendo que l´ım h→0kri,h ykh=kyk,∀y∈ Di, πh:H −→ Vhcumpliendo que l´ım h→0kπhgkh=kgk,∀g∈ H, πb h:Hb−→ Vb hcumpliendo que l´ım h→0kπb hgkb h→ kgk,∀g∈ Hb πb i,h :Hb i−→ Vb i,h cumpliendo que l´ım h→0kπb i,hgkb i,h → kgk,∀g∈ Hb i. En estos espacios consideramos aproximaciones discretas de los operadores que aparecen en el problema continuo. Estos operadores discretos heredar´an algunas de las caracter´ısticas de los operadores continuos; as´ı, tomamos Lh:Vh→Vhcomo aproximaci´on discreta de L;Li,h :Vh→Vhcomo aproximante de Lipara i= 1, . . . , n; para los operadores frontera B:D → HbyBi:Di→ Hb iconsideramos las aproximaciones discretas Bh:Vh→Vb hyBi,h :Vi,h →Vb i,h respectivamente. Con el fin de realizar el estudio de la convergencia del esquema totalmente discreto se introduce el concepto de error local de truncatura asociado a los operadores LyLien la forma: τL h(v)≡Lhrhv−πhL v, ∀v∈ D yτLi h(v)≡Li,hri,hv−πhLiv, ∀v∈ Di, y para los operadores del contorno ByBien la forma: τB h(v)≡Bhrhv−πb hB v, ∀v∈ D yτBi h(v)≡Bi,hri,hv−πb i,hLiv, ∀v∈ Di. En este contexto diremos que las aproximaciones obtenidas son consistentes de orden rsi para funciones suficientemente regulares se tiene que: kτL h(v)k=O(hr),kτLi h(v)k=O(hr),kτB h(v)k=O(hr) y kτBi h(v)k=O(hr). El esquema totalmente discreto que se obtiene viene dado por:                    ¯ Ih−τ n X i=1 Aiˆ Li,hYm,[j] h=Ym,[j] h+τ n X i=1 AiFm i,h +τAn+1 Gm h(Ym,[j] h), BhYm,[j] h=βm,[j] h,siendo βh,m,i,[j]=πb ki,hBkiYm,i,[j], Ym+1,[j] h=Ym,[j] h+τ n X i=1 biTˆ Li,hYm,[j] h+Fm i,h+τbn+1TGm h(Ym,[j] h), donde j= 0, corresponde a considerar condiciones de contorno cl´asicas y j= 1, al caso de considerar condiciones de contorno mejoradas. Para realizar un an´alisis t´ıpico de la convergencia de este esquema combinando propiedades de estabilidad y de consistencia debemos recordar que la discretizaci´on del operador (I−τ c L) se dice estable si y s´olo si k(Ih−τ c Lh)−1kh≤C, ∀h∈(0, h0]. Estas propiedades se tienen como consecuencia 7 B. Bujanda, J.C. Jorge, M.J. Moreta de preservar algunas de las propiedades relacionadas con el buen planteamiento de los problemas (4). Para estudiar la convergencia asumiremos que los siguientes problemas (relacionados con la resolubilidad de las etapas discretizadas) tienen soluci´on ´unica: (Ih−c τ Ai,h)vh=wh∈Vh,(c > 0), Bi,hvh=vb,h ∈Vb i,h, verificando adem´as que kvkh≤C(kwhkh+kvb,hkb,h), con Cindependiente de τ∈(0, τ0] y de h. Recordemos que el esquema totalmente discreto se dice que es convergente de orden pen tiempo y ren espacio si el error global (definido en el instante tmcomo Em,[j] h= krhy(tm)−Ym,[j] hkhpara j= 0,1) verifica que Em,[j] h=O(τp+hr). En el contexto que nos ocupa se puede demostrar que si la integraci´on temporal se realiza utilizando un adecuado m´etodo RKPF linealmente impl´ıcito A-estable (ver [3]) con orden cl´asico py orden de etapas q(ver [4]) y la discretizaci´on espacial es estable, consistente de orden r, las aplicaciones de conexi´on verifican las propiedades de compatibilidad entre normas, las restricciones de los operadores discretos Li,h al Ker(Bi,h) son mon´otonos y conmutativos y adem´as se cumple que k(πh−rh)ukh≤Chr, u ∈ H, entonces se tiene la siguiente cota para el error global Em,[j] h=O(hr+τm´ın{p,q+j}). Agradecimientos Este trabajo ha sido realizado con ayuda de los proyectos: MTM 2004-08012, MTM 2004-05521 y la Red Tem´atica 05/R-8. Referencias [1] I. Alonso-Mallo & B. Cano Efficient time integration of nonlinear partial differential equations by means of Rosenbrock methods, Applied Mathematics Reports, Universidad de Valladolid (2006). [2] B. Bujanda & J.C. Jorge, Efficient linearly implicit methods for nonlinear multidimensional parabolic problems, J. Comput. Appl. Math. 164/165 (2004), 159–174. [3] B. Bujanda & J.C. Jorge, Stability results for linearly implicit Fractional Step discretizations of non-linear time dependent parabolic problems, Appl. Numer. Math. 56 (2006), no. 8, 1061–1076. [4] B. Bujanda & J.C. Jorge, Order conditions for linearly implicit Fractional Step Runge-Kutta methods, IMAJNA (en prensa). [5] Sanz-Serna, J. M.; Verwer, J. G. & Hundsdorfer, W. H., Convergence and order reduction of RungeKutta schemes applied to evolutionary problems in partial differential equations, Numer. Math. 50 (1987), no. 4, 405–418. [6] Yanenko, N.N., “The method of fractional steps”, Springer, 1971. [7] Zlatev, Z. “Using efficient numerical methods in large-scale air pollution modelling”, Problems in Modern Applied Mathematics, 60–65. 8