Métodos conservativos de direcciones alternadas para problemas parabólicos semilineales sobre mallados rectangulares lógicos

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) M´etodos conservativos de direcciones alternadas para problemas parab´olicos semilineales sobre mallados rectangulares l´ogicos A. Arrar´ as, L. Portero, J.C. Jorge Dpto. de Ingenier´ıa Matem´atica e Inform´atica, Universidad P´ublica de Navarra, Campus de Arrosad´ıa s/n, 31006 Pamplona. E-mails: {andres.arraras,laura.portero,jcjorge}@unavarra.es. Palabras clave: Diferencias Finitas, Direcciones Alternadas, Elementos Finitos Mixtos, Mallado Rectangular L´ogico, Problema Parab´olico Semilineal. Resumen En este trabajo se propone y analiza un m´etodo num´erico eficiente para la resoluci´on de problemas parab´olicos semilineales que contienen derivadas mixtas y se plantean sobre dominios que pueden ser no rectangulares. El proceso de semidiscretizaci´on en espacio se efect´ua mediante un esquema de diferencias finitas centrado en celdas. Esta primera etapa de discretizaci´on es localmente conservativa y, tras ser combinada con una integraci´on temporal de tipo direcciones alternadas, nos permite obtener un esquema totalmente discreto que involucra conjuntos de sistemas tridiagonales desacoplados cuya resoluci´on es f´acilmente paralelizable. El art´ıculo incluye algunos resultados te´oricos de convergencia incondicional (de segundo orden en espacio y primer orden en tiempo) y se completa con un experimento num´erico que ilustra el comportamiento del m´etodo propuesto. 1. Introducci´on El presente art´ıculo versa sobre la resoluci´on num´erica eficiente de problemas parab´olicos de car´acter semilineal en cuya formulaci´on aparecen derivadas mixtas. En particular, consideramos un modelo evolutivo de flujo en medios porosos basado en una versi´on simplificada de la ecuaci´on de Richards (cf. [4, 7]). Dicho modelo, que supone la existencia de direcciones de flujo preferente debido a la naturaleza anis´otropa del medio, admite la 1 A. Arrar´as, L. Portero, J.C. Jorge siguiente formulaci´on: Hallar ˆ ψ:ˆ Ω×[0, T]→Rtal que          c∂ˆ ψ(ˆ x, t) ∂t +ˆ ∇ · (ˆ K(ˆ x)ˆ ∇ˆ ψ(ˆ x, t)) = ˆ f(ˆ x, t) + ˆg(ˆ ψ),(ˆ x, t)∈ˆ Ω×(0, T], ˆ ψ(ˆ x,0) = ˆ ψ0(ˆ x),ˆ x∈ˆ Ω, ˆ ψ(ˆ x, t) = ˆ ψD(ˆ x, t),(ˆ x, t)∈∂ˆ Ω×(0, T], (1) donde ˆ ψ≡ˆ ψ(ˆ x, t) es la presi´on capilar, c > 0 denota la constante de humedad espec´ıfica y ˆ K≡ˆ K(ˆ x)=(ˆ Kij(ˆ x))i,j∈{1,2}es un tensor sim´etrico y definido-positivo que representa la conductividad hidr´aulica del medio. Por su parte, ˆ f(ˆ x, t) es un t´ermino fuente/sumidero y ˆg(ˆ ψ) constituye una funci´on no lineal suave que modeliza, por ejemplo, la absorci´on radicular en perfiles de suelo (cf. [9]). Suponemos que (1) est´a definido sobre un dominio abierto y acotado ˆ Ω⊆R2de geometr´ıa no rectangular, cuya frontera se denota por ∂ˆ Ω. Finalmente, las condiciones inicial y de contorno vienen dadas por las funciones ˆ ψ0(ˆ x) y ˆ ψD(ˆ x, t), respectivamente. El m´etodo num´erico desarrollado combina los procedimientos de semidiscretizaci´on espacial e integraci´on temporal que se describen en las secciones 2 y 3, respectivamente. Su contenido se completa con algunos resultados te´oricos de convergencia incondicional para los esquemas semidiscreto y totalmente discreto. Por ´ultimo, la secci´on 4 incluye un experimento num´erico que ilustra el comportamiento del algoritmo en la simulaci´on de modelos de flujo como el formulado previamente. 2. Semidiscretizaci´on en espacio Para la resoluci´on num´erica del problema (1) sobre un dominio de flujo ˆ Ω de geometr´ıa no rectangular, proponemos una transformaci´on previa del mismo que permita plantearlo sobre un dominio computacional Ω de geometr´ıa rectangular. Sea F:R2→R2una transformaci´on inversible y suficientemente regular (al menos C2) que satisface F(x) = ˆ x para x∈Ω y ˆ x∈ˆ Ω. Denotemos por J≡J(x) a la matriz jacobiana de F, siendo d≡d(x) = |det(J)|. Si {Th}h>0representa una malla de elementos finitos rectangulares cuasi-uniforme definida sobre Ω, donde hdesigna el di´ametro m´aximo de los elementos que la componen, entonces {ˆ Th}h>0≡F(Th) ser´a una malla rectangular l´ogica definida sobre ˆ Ω. Es decir, Ftransforma la estructura rectangular de un elemento E∈ Then la estructura curvil´ınea de su correspondiente ˆ E∈ˆ Th, manteniendo la proporci´on entre Eyˆ E. Para toda funci´on escalar ϕ(x) definida sobre E, tenemos que ˆϕ(ˆ x) = F(ϕ)(ˆ x)≡ϕ◦F−1(ˆ x). Por su parte, para las funciones vectoriales aplicamos la transformaci´on de Piola (cf. [10]): dada q(x) definida sobre E, tenemos que ˆ q(ˆ x) = G(q)(ˆ x)≡¡d−1Jq¢◦F−1(ˆ x). Por tanto, es posible reformular (1) a trav´es de FyGsobre el dominio rectangular Ω del modo siguiente: Hallar ψ: Ω ×[0, T]→Rtal que        c∂ψ(x, t) ∂t +1 d(x)∇ · (K(x)∇ψ(x, t)) = f(x, t) + g(ψ),(x, t)∈Ω×(0, T], ψ(x,0) = ψ0(x),x∈Ω, ψ(x, t) = ψD(x, t),(x, t)∈∂Ω×(0, T], (2) donde K≡K(x) = dJ−1ˆ K(J−1)T. N´otese que, en general, el tensor Kes no diagonal 2 M´etodos conservativos de direcciones alternadas para problemas parab´olicos ••• ••• ••• xi− 3 2 xi− 1 2 xi+1 2 xi+3 2 xi−1xixi+1 yj− 3 2 yj− 1 2 yj+1 2 yj+3 2 yj−1 yj yj+1 - - -¾ ¾ ¾ hx i−1hx ihx i+1 6 6 6 ? ? ? hy j−1 hy j hy j+1 Figura 1: Estructura de la malla que discretiza el dominio computacional Ω. incluso en el caso de que ˆ Ksea diagonal, lo que da lugar a la aparici´on de derivadas mixtas en la nueva formulaci´on. La semidiscretizaci´on espacial del problema (2) se basa en una variante del m´etodo de elementos finitos mixtos que permite aproximar expl´ıcitamente tres variables: la presi´on ψ, su gradiente cambiado de signo ˜ uy el flujo u=K˜ u(cf. [1]). Considerando los espacios de Raviart-Thomas-N´ed´elec de orden m´as bajo, si en la formulaci´on variacional del problema semidiscreto se aproximan los productos escalares entre vectores por determinadas reglas de cuadratura (e.g. v´eanse las propuestas en [8]), es inmediato obtener un esquema de diferencias finitas centrado en celdas para la variable ψ. Dicho esquema preserva las propiedades esenciales de los m´etodos mixtos est´andar (precisi´on asint´otica y conservatividad local), reduciendo notablemente su coste computacional. A continuaci´on se realiza una descripci´on detallada del mismo. Si denotamos los nodos de la malla por (xi+1/2, yj+1/2) para i= 0,1, . . . , Nxyj= 0,1, . . . , Ny, es posible definir xi= (xi−1/2+xi+1/2)/2, hx i=xi+1/2−xi−1/2eyj= (yj−1/2+yj+1/2)/2, hy j=yj+1/2−yj−1/2para i= 1,2, . . . , Nxyj= 1,2, . . . , Ny. La Fig. 1 muestra la estructura de la malla para el dominio computacional Ω. En lo que sigue, ψh(t) es el vector formado por las aproximaciones de la funci´on escalar ψ(x, t) en los centros de las celdas, as´ı como por las evaluaciones de ψD(x, t) en los puntos medios de los segmentos de contorno. Por su parte, uh(t)≡(ux h(t), uy h(t)) y ˜ uh(t)≡ (˜ux h(t),˜uy h(t)) aproximan a las funciones vectoriales u(x, t)≡(ux(x, t), uy(x, t)) y ˜ u(x, t)≡ (˜ux(x, t),˜uy(x, t)) en los puntos medios de los segmentos entre nodos; en concreto, ux h(t) y ˜ux h(t) se localizan en los segmentos verticales, mientras que uy h(t) y ˜uy h(t) se sit´uan en los horizontales. Finalmente, Khcontiene las evaluaciones nodales de las componentes del tensor K(x). La discretizaci´on del operador lineal el´ıptico A≡ −∇ · K∇se realiza en tres fases sucesivas. La primera de ellas permite obtener una aproximaci´on de las componentes de ˜ u=−∇ψen los puntos interiores (i+1 2, j) e (i, j +1 2) de la forma ˜ux h,i+1 2,j =−ψh,i+1,j −ψh,i,j 1 2¡hx i+hx i+1¢,˜uy h,i,j+1 2 =−ψh,i,j+1 −ψh,i,j 1 2³hy j+hy j+1´,(3) 3 A. Arrar´as, L. Portero, J.C. Jorge • •× (a) × × × × × ¤ ◦ ◦ (c) • • × (b) • • • • • • • • • × × × × × × × × × × × × ¤ ¤ ¤ ¤ ◦ ◦ ◦ ◦ (d) Figura 2: Stencil asociado a: (a) ˜ux h(t), (b) ˜uy h(t), (c) ux h(t) y (d) ∇h·uh. Simbolog´ıa: •, ψh(t); ×, ˜ux h(t) y ˜uy h(t); ¤,ux h(t) y uy h(t); ◦,Kh. donde ψh,i,j ≈ψ(xi, yj, t), ˜ux h,i+1 2,j ≈˜ux(xi+1 2, yj, t) y ˜uy h,i,j+1 2 ≈˜uy(xi, yj+1 2, t). En los puntos de frontera, tenemos las siguientes expresiones: ˜ux h,1 2,j =−ψh,1,j −ψh, 1 2,j 1 2hx 1 ,˜ux h,Nx+1 2,j =−ψh,Nx+1 2,j −ψh,Nx,j 1 2hx Nx , ˜uy h,i,1 2 =−ψh,i,1−ψh,i,1 2 1 2hy 1 ,˜uy h,i,Ny+1 2 =−ψh,i,Ny+1 2−ψh,i,Ny 1 2hy Ny , (4) donde ψh,1 2,j,ψh,Nx+1 2,j,ψh,i, 1 2yψh,i,Ny+1 2se deducen de las condiciones de contorno de tipo Dirichlet. Los stencils correspondientes a las componentes xeyse muestran en las Figs. 2(a) y 2(b). En la segunda fase, obtenemos una aproximaci´on de las componentes de u=K˜ ua partir de las discretizaciones previas. As´ı, la expresi´on de ux h(t) en el punto interior (i+1 2, j) viene dada por ux h,i+1 2,j =1 2³(K11)h,i+1 2,j−1 2+ (K11)h,i+1 2,j+1 2´˜ux h,i+1 2,j +µµ(K12)h,i+1 2,j−1 2˜uy h,i+1,j−1 2 + (K12)h,i+1 2,j+1 2˜uy h,i+1,j+1 2¶hx i+1 +µ(K12)h,i+1 2,j−1 2˜uy h,i,j−1 2 + (K12)h,i+1 2,j+1 2˜uy h,i,j+1 2¶hx i¶1 2(hx i+hx i+1), (5) donde (K1`)h,i+1 2,j±1 2≈K1`(xi+1 2, yj±1 2), `= 1,2. En el contorno, se tiene ux h,i+1 2,j =1 2³³(K11)h,i+1 2,j−1 2+ (K11)h,i+1 2,j+1 2´˜ux h,i+1 2,j +µ(K12)h,i+1 2,j−1 2˜uy h,ˇı+1,j−1 2 + (K12)h,i+1 2,j+1 2˜uy h,ˇı+1,j+1 2¶¶, (6) para i= ˇı= 0 e i= ˇı+ 1 = Nx. La Fig. 2(c) muestra el stencil asociado a esta discretizaci´on. An´alogamente, es inmediato obtener las expresiones correspondientes a uy h,i,j+1/2en el interior de Ω y en su frontera. 4 M´etodos conservativos de direcciones alternadas para problemas parab´olicos La tercera fase permite calcular la aproximaci´on de ∇ · u=−∇ · (K∇ψ)≡Aψ, que denotamos por ∇h·uh, como (∇h·uh)i,j =  ux h,i+1 2,j −ux h,i−1 2,j hx i + uy h,i,j+1 2 −uy h,i,j−1 2 hy j  .(7) Si combinamos las ecuaciones (3)-(7), obtendremos una aproximaci´on de Aψque contendr´a ´unicamente a la variable ψh(t) y, por tanto, podr´a representarse en forma compacta como un stencil de nueve puntos (v´ease Fig. 2(d)). Considerando la semidiscretizaci´on descrita, el problema (2) puede aproximarse por la siguiente familia de problemas de valor inicial stiff:    cdΨh(t) dt =1 dh (AhΨh(t) + Bh(t)) + Fh(t) + Gh(Ψh), t ∈(0, T], Ψh(0) = rhψ0(x) = Ψ0 h, (8) donde Ψh(t) es un vector que contiene las aproximaciones de ψ(x, t) en los centros de las celdas. Con esta notaci´on, AhΨh(t)+Bh(t) representa la discretizaci´on del t´ermino difusivo Aψbasada en (3)-(7). N´otese que, a diferencia de ψh(t), Ψh(t) no incluye la contribuci´on de ψD(x, t) que, en el esquema (8), est´a contenida en Bh(t). Por su parte, Gh(Ψh) designa las aproximaciones en los centros de las celdas del t´ermino no lineal g(ψ). Finalmente, si denotamos por rhla restricci´on a dichos centros, tenemos que Fh(t) = rhf(x, t) y dh=rhd(x). Sea Hhun espacio finito-dimensional de funciones discretas que toman sus valores en los centros de las celdas y denotemos por k·khla norma L2discreta asociada a Hh. Asumiendo condiciones suficientes de suavidad y compatibilidad en los datos f(x, t), ψ0(x) y ψD(x, t) y suponiendo que la soluci´on ψ(x, t) es suficientemente regular, se puede demostrar que krhψ(x, t)−ψh(t)kh⩽Ch2∀t∈[0, T],donde Ces una constante independiente de h (v´ease [2]). 3. Integraci´on en tiempo La aproximaci´on en espacio da lugar al sistema (8), cuya integraci´on temporal se realiza mediante un m´etodo de direcciones alternadas linealmente impl´ıcito. Para ello, el t´ermino AhΨh(t)+ Bh(t) es descompuesto en tres sumandos: dos de ellos agrupan las derivadas de segundo orden discretas con respecto a cada una de las variables espaciales, mientras que el tercero contiene las discretizaciones de las derivadas mixtas. En particular, consideramos sendas descomposiciones de la matriz Ahy del vector Bh(t) de la forma Ah=A0h+A1h+A2hyBh(t) = B0h(t) + B1h(t) + B2h(t), de modo que A0hΨh(t)+B0h(t)≈(∂x(K12∂yψ)+∂y(K12∂xψ)), A1hΨh(t)+B1h(t)≈(∂x(K11∂xψ)) yA2hΨh(t) + B2h(t)≈(∂y(K22∂yψ)). Adem´as, particionamos Fh(t) en dos sumandos suficientemente regulares, i.e. Fh(t) = F1h(t) + F2h(t). N´otese que los vectores Bih(t), i= 0,1,2, contienen las contribuciones de las condiciones de contorno a cada uno de los operadores discretos Aih,i= 0,1,2. Si el vector de inc´ognitas Ψh(t) se ordena siguiendo las filas de la malla, la matriz A1hpresenta una estructura tridiagonal, mientras que si la ordenaci´on se realiza por columnas, dicha estructura recae en A2h. 5 A. Arrar´as, L. Portero, J.C. Jorge El integrador temporal elegido para discretizar (8) es una variante del esquema de Euler impl´ıcito fraccionario que admite la siguiente formulaci´on:          ³I−∆t dhcA1h´Ψn,1 h= Ψn h+∆t dhc¡dhFn+1 1h+dhGh(Ψn h) + A0hΨn h+Bn 0h+Bn+1 1h¢, ³I−∆t dhcA2h´Ψn,2 h= Ψn,1 h+∆t dhc¡dhFn+1 2h+Bn+1 2h¢, Ψn+1 h= Ψn,2 h,para n= 0,1, . . . , NT−1¡NT≡£T ∆t¤¢, (9) donde ∆tdenota el paso en tiempo y Ψn hes una aproximaci´on de Ψh(tn), siendo tn=n∆t. Los vectores Ψn,1 hy Ψn,2 hcontienen la soluci´on de las etapas internas, mientras que Fn ih = Fih(tn), i= 1,2, y Bn ih =Bih(tn), i= 0,1,2. Al igual que en los esquemas cl´asicos de direcciones alternadas, los operadores discretos A1hyA2hact´uan en (9) de modo impl´ıcito; sin embargo, tanto A0hcomo el t´ermino no lineal discreto Gh(Ψn h) contribuyen al esquema expl´ıcitamente. El tratamiento expl´ıcito de los t´erminos que discretizan las derivadas mixtas ha sido estudiado en [5, 6] para algunos esquemas de direcciones alternadas aplicados a problemas parab´olicos lineales. Por su parte, en [3] se analiza la contribuci´on expl´ıcita de t´erminos no lineales para una familia de m´etodos Runge-Kutta de pasos fraccionarios linealmente impl´ıcitos. El esquema totalmente discreto as´ı obtenido est´a compuesto por un conjunto de sistemas lineales (uno por cada etapa interna), cada uno de los cuales puede ser descompuesto en subsistemas tridiagonales desacoplados f´acilmente paralelizables. El an´alisis cl´asico de convergencia del esquema (9) implica un estudio previo de las propiedades de consistencia y estabilidad. Para analizar la consistencia, introducimos el concepto de error local en el instante de tiempo tn+1 como sigue: en+1 h= Ψh(tn+1)−˘ Ψn+1 h, donde ˘ Ψn+1 hes la soluci´on num´erica obtenida con (9) considerando Ψn h≡Ψh(tn). Bajo hip´otesis adecuadas de regularidad y compatibilidad en los datos, es f´acil probar que ken+1 hkh⩽C(∆t)2∀∆t∈(0,∆t0] y n= 0,1, . . . , NT−1,donde Ces una constante independiente de hy ∆t. Por tanto, la variante del m´etodo de Euler impl´ıcito fraccionario formulada en (9) es uniformemente consistente de orden uno. Para estudiar la estabilidad, definimos en primer lugar el operador de transici´on como Rh³∆t dhcA0h,∆t dhcA1h,∆t dhcA2h´=³I−∆t dhcA2h´−1³I−∆t dhcA1h´−1³I+∆t dhcA0h´. Tal como se describe en [2], es posible demostrar la estabilidad del integrador temporal aplicado a un problema parab´olico lineal de coeficientes constantes mediante un an´alisis de von Neumann, obteniendo que ° ° °Rh³∆t dhcA0h,∆t dhcA1h,∆t dhcA2h´° ° °h ⩽1.(10) En lo que sigue, asumiremos que esta propiedad es extensible a un problema de tipo (8). Introducimos a continuaci´on el esquema perturbado          ³I−∆t dhcA1h´˜ Ψn,1 h=˜ Ψn h+∆t dhc³dhFn+1 1h+dhGh(˜ Ψn h) + A0h˜ Ψn h+Bn 0h+Bn+1 1h´+δn,1 h, ³I−∆t dhcA2h´˜ Ψn,2 h=˜ Ψn,1 h+∆t dhc¡dhFn+1 2h+Bn+1 2h¢+δn,2 h, ˜ Ψn+1 h=˜ Ψn,2 h,para n= 0,1, . . . , NT−1,siendo ˜ Ψ0 h= Ψ0 h+δ0 h. 6 M´etodos conservativos de direcciones alternadas para problemas parab´olicos Suponiendo que la funci´on no lineal g(ψ) satisface una condici´on de Lipschitz, i.e. kg(ψ1)− g(ψ2)k⩽Lkψ1−ψ2k ∀ ψ1, ψ2: Ω ×[0, T ]→R, y que el operador de transici´on satisface (10), es inmediato demostrar la siguiente cota de estabilidad (cf. [2]): k˜ Ψn+1 h−Ψn+1 hkh⩽C kδ0 hkh+ n X j=0 ³kδj,1 hkh+kδj,2 hkh´ ,para n= 0,1, . . . , NT−1. Finalmente, si definimos el error global en el instante de tiempo tn+1 como En+1 h= rh(ψ(x, tn+1)) −Ψn+1 hy combinamos los resultados de convergencia del esquema semidiscreto con los de consistencia uniforme y estabilidad del integrador temporal, es sencillo demostrar que kEn+1 hkh≤C(h2+ ∆t)∀∆t∈(0,∆t0] y n= 0,1, . . . , NT−1. 4. Experimento num´erico En esta secci´on aplicamos el m´etodo descrito a un problema parab´olico de tipo (1) en el que (ˆ x, t)∈ˆ Ω×[0,1], donde ˆ Ω≡ {(ˆx, ˆy) : 1 <ˆx2+ ˆy2<4,ˆx > 0,ˆy > 0}. Denotamos por ˆ K≡(ˆ Kij)i,j∈{1,2}el tensor de conductividad hidr´aulica cuyos coeficientes son ˆ K11 =ˆ K22 = 2 y ˆ K12 =ˆ K21 = 1, siendo c= 1 la constante de humedad espec´ıfica. El t´ermino no lineal ˆg(ˆ ψ) = 1/(1 + ˆ ψ3) reproduce la funci´on de absorci´on radicular asociada al modelo de infiltraci´on en suelos de van Genuchten. Finalmente, el t´ermino fuente ˆ f(ˆ x, t) y las condiciones inicial ˆ ψ0(ˆ x) y de contorno ˆ ψD(ˆ x, t) han sido elegidos de tal modo que ˆ ψ(ˆ x, t) = e−t 8¡ˆxˆy¡ˆx2+ ˆy2−1¢¡ˆx2+ ˆy2−4¢¢2es la soluci´on exacta del problema. Aplicando la transformaci´on F(x)=(xcos y, x sin y), definida del dominio computacional Ω ≡ {(x, y)∈R2: 1 <x<2,0< y < π/2}en el dominio de flujo ˆ Ω, obtenemos un nuevo tensor K(x)≡(Kij(x))i,j∈{1,2}con coeficientes variables K11(x) = x(2 + sin (2 y)), K12(x) = K21(x) = cos (2 y) y K22(x) = x−1(2 −sin (2 y)), as´ı como las correspondientes funciones transformadas g(ψ), f(x, t), ψ0(x) y ψD(x, t). A continuaci´on hacemos uso del esquema totalmente discreto (9) para aproximar num´ericamente la soluci´on del problema dado. Con el fin de estudiar el comportamiento asint´otico del m´etodo, calcularemos los errores globales en norma L∞(0,1; L2(Ω)) discreta, definidos como E∆t Nx,Ny= m´ax n∈{1,2,...,NT}  Nx X i=1 Ny X j=1 hxhy¡ψ(xi, yj, tn)−ψn ij¢2  1/2 , asumiendo que los par´ametros hx iyhy json constantes e iguales a hx= 1/Nxyhy= π/(2Ny), respectivamente. El t´ermino ψ(xi, yj, tn) denota la soluci´on exacta evaluada en (xi, yj) para tn=n∆t, mientras que ψn ij designa la soluci´on num´erica obtenida en dicha posici´on espacial para el mismo instante de tiempo. La Tabla 1 contiene los errores globales E∆t Nx,Ny(fila superior) y ´ordenes num´ericos de convergencia p(fila inferior) para diferentes mallados en los que Nx=Ny. Los ´ordenes de convergencia han sido obtenidos como p= log2(E∆t Nx,Ny/E∆t/4 2Nx,2Ny). Al objeto de preservar el mismo orden de magnitud en los errores en espacio y tiempo, hemos impuesto la relaci´on ∆t NxNy= 10. De los resultados expuestos, se puede deducir que el algoritmo es convergente de segundo orden en espacio y primer orden en tiempo. 7 A. Arrar´as, L. Portero, J.C. Jorge Finalmente, la Tabla 2 describe la influencia de la discretizaci´on temporal en el esquema num´erico, de nuevo en t´erminos de errores globales y ´ordenes num´ericos de convergencia. En este caso, se considera una malla suficientemente fina (Nx=Ny= 300), con el fin de que los errores asociados a la discretizaci´on espacial puedan asumirse despreciables. La estructura de esta tabla es an´aloga a la de la Tabla 1, siendo p= log2(E∆t Nx,Ny/E∆t/2 Nx,Ny). Los resultados mostrados ilustran la estabilidad incondicional del esquema, as´ı como su convergencia de primer orden en tiempo. Nx=Ny20 40 80 160 320 ∆t∆t0∆t0/4 ∆t0/16 ∆t0/64 ∆t0/256 E∆t Nx,Ny6,833E-1 1,860E-1 4,760E-2 1,197E-2 2,998E-3 p1,8772 1,9663 1,9915 1,9974 − Tabla 1: Errores globales y ´ordenes num´ericos de convergencia para ∆t0= 2,5·10−2. ∆t∆t0∆t0/2 ∆t0/4 ∆t0/8 ∆t0/16 ∆t0/32 E∆t Nx,Ny3,562E-1 1,937E-1 1,012E-1 5,190E-2 2,629E-2 1,326E-2 p0,8789 0,9366 0,9634 0,9812 0,9874 − Tabla 2: Errores globales y ´ordenes num´ericos de convergencia para Nx=Ny= 300 y ∆t0= 2,0·10−2. Referencias [1] T. Arbogast, M.F. Wheeler, I. Yotov. Mixed finite elements for elliptic problems with tensor coefficients as cell-centered finite differences. SIAM J. Numer. Anal., 34 (1997), 828–852. [2] A. Arrar´as, L. Portero, J.C. Jorge. A combined mixed finite element ADI scheme for solving Richards’ equation with mixed derivatives on irregular grids. Preprint. [3] B. Bujanda, J.C. Jorge. Stability results for linearly implicit fractional step discretizations of nonlinear time dependent parabolic problems. Appl. Numer. Math., 56 (2006), 1061–1076. [4] M.A. Celia, E.T. Bouloutas, R.L. Zarba. A general mass-conservative numerical solution for the unsaturated flow equation. Water Resour. Res., 26 (1990), 1483–1496. [5] K.J. in’t Hout, B.D. Welfert. Stability of ADI schemes applied to convection-diffusion equations with mixed derivative terms. Appl. Numer. Math., 57 (2007), 19–35. [6] S. McKee, D.P. Wall, S.K. Wilson. An alternating direction implicit scheme for parabolic equations with mixed derivative and convective terms. J. Comput. Phys., 126 (1996), 64–76. [7] L.A. Richards. Capillary conduction of liquids through porous mediums. Physics, 1 (1931), 318–333. [8] T.F. Russell, M.F. Wheeler. Finite element and finite difference methods for continuous flows in porous media. In: The Mathematics of Reservoir Simulation, Frontiers in Applied Mathematics 1, SIAM, Philadelphia, PA, 1983, 35–106. [9] J. ˇ Sim˚unek, J.W. Hopmans, J.A. Vrugt, M.T. van Wijt. One-, two- and three-dimensional root water uptake functions for transient modeling. Water Resour. Res., 37 (2001), 2457–2470. [10] J.M. Thomas, Sur l’analyse num´erique des m´ethodes d’´el´ements finis hybrides et mixtes, Th`ese d’Etat, Universit´e Pierre et Marie Curie, Paris, 1977. 8