Approximation of the scalar convection-diffusion-reaction equation with stabilized finite element formulations of high order

==Resumen==

En este artículo presentamos formulaciones estabilizadas de elementos finitos para resolver la ecuación de convección-difusión-reacción escalar en el caso de difusión pequeña. Por un lado, resumimos la aplicación de las formulaciones ASGS y OSS basadas en el concepto de métodos variacionales multiescala. Por otro lado, discutimos aspectos de la utilización de elementos de alto orden, centrando nuestros experimentos numéricos en elementos cuadráticos, cúbicos y de cuarto orden. Asimismo, la aplicación del método OSS requiere de la introducción de una proyección, para lo cual introducimos una modificación del elemento simplicial de cuarto orden con una regla de integración numérica asociada.

'''Palabras clave''': Convección-difusión-reacción, métodos de elementos finitos estabilizados, convección dominante, alto orden, cuadraturas nodales

==Abstract==

In this paper we present stabilized finite element formulations of to solve the scalar convection-diffusion-reaction equation in the case of small diffusion. On the one hand, we summarize the application of the ASGS and OSS formulations based on the concept of  variational multiscale methods. On the other hand, we discuss aspects of the use of high order elements, focusing our numerical experiments on quadratic, cubic and fourth order elements. Likewise, the application of the OSS method requires the introduction of a projection, for which we introduce a modification of the fourth order simplicial element with an associated numerical integration rule. 

'''Keywords''': Convection-diffusion-reaction, stabilized finite element methods, dominant convection, high order, nodal quadratures

==1. Introducción==

El contar con equipos computacionales más robustos tanto en memoria RAM como en velocidad de los procesadores ha permitido entre otras cosas incentivar a numerosos investigadores al estudio de métodos numéricos complejos en fluidos que abordan modelos constitutivos cada vez más reales, buscando resolver las ecuaciones sin hacer hipótesis de simplificación por su complejidad. Dentro de este marco, la ecuación de convección-difusión-reacción (CDR) es un modelo matemático de simulación del transporte de una magnitud física que se utiliza entre otras cosas en la simulación del transporte de contaminantes.

Es bien conocido que el método estándar de elementos finitos de Galerkin aplicado a la ecuación de CDR presenta inestabilidades en la solución cuando el término convectivo es dominante frente al término difusivo, motivo por el cual distintos autores han desarrollado formulaciones estabilizadas. Dicha estabilización consiste en añadir términos que dependen de la malla de elementos finitos a los términos de Galerkin. Así, se pueden encontrar en la literatura diferentes métodos, como el método Streamline-upwind/Petrov-Galerkin (SUPG) <span id='citeF-1'></span><span id='citeF-2'></span><span id='citeF-3'></span>[[#cite-1|[1,2,3]]], el método de Galerkin-mínimos cuadrados (GLS) <span id='citeF-4'></span><span id='citeF-5'></span><span id='citeF-6'></span>[[#cite-4|[4,5,6]]], el método de las características-Galerkin (CG) <span id='citeF-7'></span>[[#cite-7|[7]]] y el método de Taylor-Galerkin (TG) <span id='citeF-8'></span>[[#cite-8|[8]]], entre otros. Una descripción y comparación de estos métodos se puede encontrar en <span id='citeF-9'></span>[[#cite-9|[9]]], en donde se expone que todos los métodos esencialmente consisten en la adición de un término de estabilización a la formulación original de Galerkin, y que fundamentalmente excepto por menores modificaciones este término se puede escribir como un parámetro numérico, denominado tiempo intrínseco, multiplicado al producto <math display="inline">L^{2}</math> entre el residuo de la ecuación diferencial a ser resuelta y un operador aplicado a la función de prueba.

En el presente trabajo nos centraremos en los métodos llamados variacionales multiescalas (VMS, por ''Variational Multi-Scale''), introducidos por Hughes y otros <span id='citeF-10'></span><span id='citeF-11'></span>[[#cite-10|[10,11]]] (véase <span id='citeF-12'></span>[[#cite-12|[12]]] para una revisión de esta teoría). En particular, nos centraremos en la versión más común, a la que nos referiremos como método ASGS (por ''Algebraic Sub-Grid Scale''), y al llamado método OSS (por ''Orthogonal Subscale Stabilization'') <span id='citeF-13'></span><span id='citeF-14'></span><span id='citeF-15'></span><span id='citeF-16'></span><span id='citeF-17'></span>[[#cite-13|[13,14,15,16,17]]].

Además de describir resumidamente los métodos ASGS y OSS, la contribución fundamental del presente trabajo es analizar su comportamiento con elementos finitos de alto orden, en particular cuadráticos, cúbicos y de cuarto orden, en todos los casos de clase <math display="inline">C^{0}</math>. Para estudiar dicho comportamiento, presentaremos resultados de convergencia en malla usando soluciones manufacturadas y analizaremos cómo son de robustos estos métodos en presencia de capas límite de la solución. Asimismo, en el caso del método OSS es necesario llevar a cabo una proyección en el espacio de elementos finitos, para lo cual resulta especialmente conveniente usar reglas de integración numérica con los puntos de integración en los nodos. Para elementos simpliciales de cuarto orden esto nos lleva a proponer una modificación del elemento estándar, que describimos en detalle para el caso del elemento triangular.

El artículo está organizado como sigue. En la siguiente sección se describe el problema a resolver y las aproximaciones y formulaciones de los métodos variacionales estabilizados ASGS y OSS. En la Sección 3 presentamos el elemento finito triangular de cuarto orden modificado. Las pruebas de convergencia en malla con soluciones analíticas conocidas y las pruebas de estabilidad con capas límite, así como un ejemplo práctico, se presentan en la Sección 4. Finalmente, en base a los resultados obtenidos  se extraen algunas conclusiones en la Sección&nbsp;5.

==2. Planteamiento del problema, aproximación y estabilización numérica==

En esta sección vamos a describir la aproximación y estabilización numérica de la ecuación escalar transitoria de CDR, con condiciones de contorno  y condiciones iniciales conocidas.

===2.1 Problema de contorno===

Sea <math display="inline">\Omega </math> un dominio acotado de <math display="inline">\mathbb{R}^d</math> <math display="inline">(d= 1,2,3)</math> y sea <math display="inline">\gamma =\left[0,T\right]</math>, con <math display="inline">T > 0</math>, el dominio temporal. El problema que nos planteamos consiste en encontrar <math display="inline">u:\Omega \times \gamma \rightarrow \mathbb{R}</math> tal que

<span id="eq-1"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\partial _{t}u+\mathcal{L}(u)=f,\qquad \left(\mathbf{x},t\right)\in \Omega \times \gamma , </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

con

<span id="eq-2"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{L}(u):=-k\partial _{i}\partial _{i}u+a_{i}\partial _{i}u+su, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

sujeta a condiciones de contorno de Dirichlet

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>u=u_{D},\qquad \left(\mathbf{x},t\right)\in \partial \Omega _{D}\times \mathbf{\gamma }, </math>
|}
|}

condiciones de contorno de Neumann

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>kn_{i}\partial _{i}u=\varphi ,\qquad \left(\mathbf{x},t\right)\in \partial \Omega _{N}\times \mathbf{\gamma ,} </math>
|}
|}

y condiciones iniciales

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>u=u_{0},\qquad \mathbf{x}\in \Omega ,\quad \mbox{ }t=0. </math>
|}
|}

En estas ecuaciones, <math display="inline">\partial _{t}</math> denota la derivada temporal y <math display="inline">\partial _{i}</math> la derivada con respecto a la <math display="inline">i</math>-ésima coordenada cartesiana <math display="inline">(i=1,\dots ,d)</math>. Hemos usado la notación de Eintein, de manera que los índices repetidos indican suma sobre todas las coordenadas espaciales. Los coeficientes en ([[#eq-2|2]]) son <math display="inline">k > 0</math> la difusión, <math display="inline">a_{i}</math> la componente en la dirección <math display="inline">i</math> del campo convectivo de velocidades <math display="inline">\mathbf{a}\in \mathbb{\mathbb{R}}^{d}</math> (para flujo incompresible <math display="inline">\nabla \cdot \mathbf{a}=0</math>) y <math display="inline">s\geq{0}</math> el coeficiente de reacción. <math display="inline">f</math> es el término fuente conocido y <math display="inline">n_{i}</math> es la componente <math display="inline">i</math>-ésima de la normal unitaria a la frontera, denotada por <math display="inline">\partial \Omega </math>, la cual está dada por <math display="inline">\partial \Omega =\partial \Omega _{D}\cup \partial \Omega _{N}.</math> Las condiciones de contorno de Dirichlet <math display="inline">u_{D}</math> y las condiciones de contorno de Neumann <math display="inline">\varphi </math> son conocidas en las fronteras <math display="inline">\partial \Omega _{D}</math> y <math display="inline">\partial \Omega _{N}</math>, respectivamente. Las condiciones iniciales y de contorno de Dirichlet verifican la siguiente condición de compatibilidad <math display="inline">u_{0}\left(\mathbf{x}\right)|_{\partial \Omega _{D}}=u_{D}\left(\mathbf{x},0\right).</math>

===2.2 Forma variacional o débil del problema===

Sea <math display="inline">\mathcal{W}</math> el espacio de funciones de <math display="inline">H^1(\Omega )</math> (funciones en <math display="inline">L^2(\Omega )</math> con derivadas en <math display="inline">L^2(\Omega )</math>) que sea anulan en <math display="inline">\partial \Omega _D</math>, y sea <math display="inline">L^{2}\left(0,T;\mathcal{W}\right)</math> el espacio de funciones en <math display="inline">\Omega \times \gamma </math> tales que la norma <math display="inline">H^1(\Omega )</math> es <math display="inline">L^2(\gamma )</math>. Para simplificar la exposición, supongamos que <math display="inline">u_D = 0</math>.

Multiplicando ([[#eq-1|1]]) por una función <math display="inline">v\in \mathcal{W}</math> e integrando el término de difusión por partes, se obtiene la forma débil o variacional del problema, la cual consiste en encontrar <math display="inline">u\in L^{2}\left(0,T;\mathcal{W}\right)</math> tal que

<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left(\partial _{t}u,v\right)+B(u,v)=\left\langle f,v\right\rangle + \left\langle \varphi ,v\right\rangle_{\partial \Omega _N} \qquad \forall v\in \mathcal{W}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

donde

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>B(u,v)  =\int _{\Omega }k\partial _{i}v\partial _{i}u+\int _{\Omega }va_{i}\partial _{i}u+\int _{\Omega }vsu,</math>
|-
| style="text-align: center;" | <math> \left(\partial _{t}u,v\right) = \int _{\Omega }v\partial _{t}u,</math>
|-
| style="text-align: center;" | <math> \left\langle f,v\right\rangle  =\int _{\Omega }vf,</math>
|-
| style="text-align: center;" | <math> \left\langle \varphi ,v\right\rangle_{\partial \Omega _N} =\int _{\partial \Omega _N}v \varphi , </math>
|}
|}

y cumpliéndose la condición inicial <math display="inline">u = u_0</math> en <math display="inline">t = 0</math> en <math display="inline">L^2(\Omega )</math>. En lo que sigue, tomaremos <math display="inline">\varphi = 0</math> para simplificar la notación.

===2.3 Aproximación numérica con el método estándar de Galerkin===

La ecuación <math display="inline">\left(\right. </math>[[#eq-3|3]]<math>\left. \right)</math> es la que aproximamos numéricamente en el espacio usando una formulación de elementos finitos. La discretización temporal la llevaremos a cabo mediante un método de diferencias finitas. En principio, podríamos discretizar primero en espacio o en tiempo indistintamente, aunque en nuestro caso llevaremos a cabo primero la discretización temporal.

====2.3.1 Discretización temporal en diferencias finitas====

Consideremos una partición uniforme del intervalo de tiempo <math display="inline">0=t^{0}<t^{1}<\dots{<}t^{N}=T</math>, con <math display="inline">\delta t=t^{n+1}-t^{n}</math>, siendo <math display="inline">\delta t</math> es el tamaño del paso de tiempo, considerado constante. Aquí y en adelante, usaremos el superíndice <math display="inline">n</math> para indicar el nivel de tiempo&nbsp;<math display="inline">n</math>.

Aunque la discretización en el tiempo puede realizarse mediante cualquier aproximación de la derivada temporal en diferencias finitas o elementos finitos, aquí utilizaremos la regla trapezoidal generalizada. Para una función genérica <math display="inline">f</math>, sea

<span id="eq-4"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>f^{n+\theta }=\theta f^{n+1}+\left(1-\theta \right)f^{n},\qquad{0}\leq \theta \leq{1.} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

Aproximaremos

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta _{t}^{n}f:=\frac{1}{\delta t}(f^{n+1}-f^{n})=\frac{1}{\theta \delta t}(f^{n+\theta }-f^{n})\approx (\partial _{t}f)^{n+\theta } . </math>
|}
|}

Esta aproximación es de segundo orden en <math display="inline">\delta t</math> en el caso en que <math display="inline">\theta=1/2</math> (método de Crank-Nicolson) y de primer orden si <math display="inline">\theta \not=1/2</math>. El caso <math display="inline">\theta=1</math> corresponde al método de Euler implícito o el método de diferencias hacia atrás de primer orden, a menudo abreviado como BDF1. Puesto que nuestro interés en este trabajo se centra en la aproximación espacial, en los ejemplos numéricos hemos usado este último método.

En cada paso de tiempo, el problema <math display="inline">\left(\right. </math>[[#eq-3|3]]<math>\left. \right)</math> discretizado en el tiempo, consiste en encontrar <math display="inline">u^{n+1}\in \mathcal{W}</math> tal que

<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left(\delta _{t}^{n}u,v\right)+B(u^{n+\theta },v)=\left\langle f^{n+\theta },v\right\rangle \qquad \forall v\in \mathcal{W}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

====2.3.2 Discretización de Galerkin en el espacio====

Para aproximar numéricamente la ecuación <math display="inline">(</math>[[#eq-5|5]]<math>)</math> en el espacio por el método de elementos finitos, consideremos una partición del dominio <math display="inline">\Omega </math>, de manera que <math display="inline">\overset{n_{el}}{\underset{e=1}{\cup }}\Omega ^{e}=\Omega </math> y <math display="inline">\Omega ^{e_{1}}\cap \Omega ^{e_{2}}=\emptyset </math> para cualesquiera <math display="inline">e_{1}</math>, <math display="inline">e_{2}</math>, <math display="inline">e_{1}\neq e_{2}</math>, donde <math display="inline">n_{el}</math> es el número de elementos de la partición y <math display="inline">\Omega ^{e}</math> es el dominio del elemento <math display="inline">e</math>.

Si <math display="inline">\mathcal{W}_{h}</math> es un espacio conforme de elementos finitos para aproximar <math display="inline">\mathcal{W}</math> (<math display="inline">\mathcal{W}_{h}\subset \mathcal{W}</math>), entonces el problema discreto que se conoce como el método estándar de Galerkin, que resuelve la ecuación escalar transitoria de CDR <math display="inline">(</math>[[#eq-1|1]]<math>)</math> con condiciones de Dirichlet y de Neumann homogéneas en la frontera, consiste en: dado <math display="inline">u_{h}^{n}\in \mathcal{W}_{h},</math> encontrar <math display="inline">u_{h}^{n+1}\in \mathcal{W}_{h}</math> tal que

<span id="eq-6"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left(\delta _{t}^{n}u_{h},v_{h}\right)+B(u_{h}^{n+\theta },v_{h})=\left\langle f^{n+\theta },v_{h}\right\rangle \qquad \forall v_{h}\in \mathcal{W}_{h}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

para <math display="inline">n = 0,1,\dots ,N-1</math>, siendo <math display="inline">u^0_h</math> conocido proyectando la condición inicial en <math display="inline">\mathcal{W}_h</math> (con el producto escalar de <math display="inline">L^2(\Omega )</math>).

===2.4 Métodos de estabilización numérica===

Como se mencionó anteriormente, el método estándar de Galerkin presenta inestabilidades numéricas cuando el término convectivo es dominante con respecto al término difusivo. Los métodos que vamos a utilizar para la estabilización numérica se pueden enmarcar en el método variacional multiescalas (VMS) <span id='citeF-10'></span><span id='citeF-11'></span>[[#cite-10|[10,11]]]. En particular, veremos dos opciones, correspondientes a dos elecciones del espacio de subescalas (véase también <span id='citeF-13'></span><span id='citeF-15'></span><span id='citeF-16'></span>[[#cite-13|[13,15,16]]]).

====2.4.1 El método variacional mutiescalas (VMS)====

Aunque el método VMS puede considerar un número arbitrario de escalas en la solución, en la mayoría de las ocasiones basta con considerar dos escalas para diseñar un método de estabilización. La idea básica es descomponer para cada instante de tiempo la variable continua desconocida <math display="inline">u(\cdot ,t)\in \mathcal{W}</math> en una parte resoluble en el espacio de elementos finitos <math display="inline">u_{h}(\cdot ,t)\in \mathcal{W}_{h}</math>, con <math display="inline">\mathcal{W}_{h}\subset \mathcal{W}</math>, y una parte en la escala de una submalla <math display="inline">\tilde{u}(\cdot ,t) \in \tilde{\mathcal{W}}</math>, con <math display="inline">\tilde{\mathcal{W}}\subset \mathcal{W}</math>,  la cual no puede ser capturada por la malla de elementos finitos, siendo <math display="inline">\mathcal{W}=\mathcal{W}_{h}\oplus \tilde{\mathcal{W}}</math>, donde <math display="inline">\mathcal{\tilde{W}}</math> es cualquier espacio para completar <math display="inline">\mathcal{W}_{h}</math> en <math display="inline">\mathcal{W}</math>. Para evitar tecnicismos, podemos pensar en <math display="inline">\mathcal{W}</math> y <math display="inline">\mathcal{\tilde{W}}</math> como espacios dimensionalmente finitos, con una dimensión grande. Puesto que <math display="inline">\tilde{\mathcal{W}}</math> representa la componente de <math display="inline">\mathcal{W}</math> que no es reproducida por el espacio de elementos finitos, la llamamos espacio de las subescalas. De este modo, la variable continua desconocida podemos escribirla como <math display="inline">u=u_{h}+\tilde{u}</math>, donde <math display="inline">u_{h}</math> es la componente de <math display="inline">u</math> en el espacio de elementos finitos y <math display="inline">\tilde{u}</math> es su componente en el complemento (con respecto a un cierto producto interno) de <math display="inline">\mathcal{W}</math>. La misma descomposición es aplicable a las funciones de test.

Aplicando la descomposición mencionada a la ecuación continua <math display="inline">\left(\right. </math>[[#eq-3|3]]<math>\left. \right)</math> se obtiene:

<span id="eq-7"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left(\partial _{t}\left(u_{h}+\tilde{u}\right),v_{h}+\tilde{v}\right)+B\left(u_{h}+\tilde{u},v_{h}+\tilde{v}\right)=\left\langle f,v_{h}+\tilde{v}\right\rangle \qquad \forall v{}_{h}\in \mathcal{W}_{h},\;\forall \tilde{v}\in \tilde{\mathcal{W}}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

Utilizando la linealidad de las formas involucradas, la ecuación <math display="inline">\left(\right. </math>[[#eq-7|7]]<math>\left. \right)</math> se puede descomponer en el siguiente sistema de ecuaciones:

<span id="eq-8"></span>
<span id="eq-9"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left(\partial _{t} u_{h},v_{h}\right)+\left(\partial _{t} \tilde{u},v_{h}\right)+B\left(u_{h},v_{h}\right)+B\left(\tilde{u},v_{h}\right) =\left\langle f,v_{h}\right\rangle \qquad \forall v_{h}\in \mathcal{W}_{h},</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|-
| style="text-align: center;" | <math> \left(\partial _{t} u_{h},\tilde{v}\right)+\left(\partial _{t}\tilde{u},\tilde{v}\right)+B\left(u_{h},\tilde{v}\right)+B\left(\tilde{u},\tilde{v}\right) =\left\langle f,\tilde{v}\right\rangle \qquad \forall \tilde{v}\in \tilde{\mathcal{W}}. </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}
|}

Para simplificar el método resultante, ''vamos a asumir que <math>\partial _t \tilde{u}=0</math>'',  lo cual implica que la variación temporal de las subescalas en <math display="inline">\left(\right. </math>[[#eq-8|8]]<math>\left. \right)</math> y <math display="inline">\left(\right. </math>[[#eq-9|9]]<math>\left. \right)</math> es despreciable comparada con el resto de términos. En <span id='citeF-17'></span>[[#cite-17|[17]]] las subescalas que satisfacen esta condición son llamadas cuasi-estáticas. La inclusión o no de la derivada temporal de las subescalas se analiza en detalle en <span id='citeF-18'></span>[[#cite-18|[18]]], donde se muestra que es esencial para evitar las inestabilidades que pueden producirse para pasos de tiempo pequeños, que en este trabajo no consideraremos.

Con esta suposición el sistema de ecuaciones <math display="inline">\left(\right. </math>[[#eq-8|8]]<math>\left. \right)</math> y <math display="inline">\left(\right. </math>[[#eq-9|9]]<math>\left. \right)</math> se escribe de la siguiente manera:

<span id="eq-10"></span>
<span id="eq-11"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left(\partial _{t} u_{h},v_{h}\right)+B\left(u_{h},v_{h}\right)+B\left(\tilde{u},v_{h}\right) =\left\langle f,v_{h}\right\rangle \qquad \forall v_{h}\in \mathcal{W}_{h},</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|-
| style="text-align: center;" | <math> \left(\partial _{t} u_{h},\tilde{v}\right)+B\left(u_{h},\tilde{v}\right)+B\left(\tilde{u},\tilde{v}\right) =\left\langle f,\tilde{v}\right\rangle \qquad \forall \tilde{v}\in \tilde{\mathcal{W}}. </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}
|}

La ecuación <math display="inline">\left(\right. </math>[[#eq-10|10]]<math>\left. \right)</math> corresponde a la escala resoluble <math display="inline">u_{h}</math> en el espacio de elementos finitos y tiene tres términos a su lado izquierdo, donde el primero y segundo términos son la contribución temporal y espacial de <math display="inline">u_{h}</math> del método estándar de Galerkin y el tercero toma en cuenta la influencia de la subescala en <math display="inline">u_{h}</math>.

De la ecuación <math display="inline">\left(\right. </math>[[#eq-11|11]]<math>\left. \right)</math> se obtiene <math display="inline">\tilde{u}</math>, que es la contribución de la subescala sobre la componente en el espacio de elementos finitos. Para evitar aproximar derivadas de <math display="inline">\tilde{u}</math>, podemos considerar la siguiente integración por partes en cada elemento:

<span id="eq-12"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>B\left(\tilde{u},v_{h}\right) =\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\Omega ^{e}}\left(-k\partial _{i}\left(\partial _{i}v_{h}\right)\tilde{u}-\partial _{i}v_{h}a_{i}\tilde{u}+v_{h}s\tilde{u}\right)</math>
|-
| style="text-align: center;" | <math>  +\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\partial \Omega ^{e}}k\partial _{i}v_{h}n_{i}\tilde{u}+\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\partial \Omega ^{e}}v_{h}n_{i}a_{i}\tilde{u}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

Sea <math display="inline">\mathcal{L}^{*}</math> el adjunto del operador <math display="inline">\mathcal{L}</math>,  dado por

<span id="eq-13"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{L}^{*}\left(v_{h}\right)=-k\partial _{i}\left(\partial _{i}v_{h}\right)-a_{i}\partial _{i}v_{h}+sv_{h}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

La ecuación <math display="inline">\left(\right. </math>[[#eq-12|12]]<math>\left. \right)</math> se puede escribir como

<span id="eq-14"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>B\left(\tilde{u},v_{h}\right)=\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\Omega ^{e}}\mathcal{L}^{*}\left(v_{h}\right)\tilde{u}+\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\partial \Omega ^{e}}k\partial _{i}v_{h}n_{i}\tilde{u}+\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\partial \Omega ^{e}}v_{h}n_{i}a_{i}\tilde{u}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

Reemplazando <math display="inline">\left(\right. </math>[[#eq-14|14]]<math>\left. \right)</math> en la ecuación <math display="inline">\left(\right. </math>[[#eq-10|10]]<math>\left. \right)</math>, ésta se puede escribir como

<span id="eq-15"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left(\partial _{t}u_{h},v_{h}\right)+B\left(u_{h},v_{h}\right)+\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\Omega ^{e}}\mathcal{L}^{*}\left(v_{h}\right)\tilde{u}+\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\partial \Omega ^{e}}k\partial _{i}v_{h}n_{i}\tilde{u}+\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\partial \Omega ^{e}}v_{h}n_{i}a_{i}\tilde{u}=\left\langle f,v_{h}\right\rangle   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

para todo <math display="inline">v_h\in \mathcal{W}_h</math>. Ahora bien integrando por partes los términos <math display="inline">B\left(u_{h},\tilde{v}\right)</math> y <math display="inline">B\left(\tilde{u},\tilde{v}\right)</math> de la ecuación <math display="inline">\left(\right. </math>[[#eq-11|11]]<math>\left. \right)</math>, ésta se puede escribir como

<span id="eq-16"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\partial \Omega ^{e}}k\tilde{v}n_{i}\partial _{i}\left(u_{h}+\tilde{u}\right)+\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\Omega ^{e}}\tilde{v}\mathcal{L}\left(\tilde{u}\right) =\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\Omega ^{e}}\tilde{v}\left[f-\mathcal{L}\left(u_{h}\right)\right]-\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\Omega ^{e}}\tilde{v}\partial _{t} u_{h}\qquad \forall \tilde{v}\in \tilde{\mathcal{W}}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

Observemos que el primer término de la ecuación <math display="inline">\left(\right. </math>[[#eq-16|16]]<math>\left. \right)</math> se anula, ya que en la frontera de los elementos las componentes normales de los flujos de <math display="inline">u</math> son iguales pero de signo contrario, debido a que los flujos difusivos deben ser continuos a través de los contornos ínter-elementales. Como resultado de esto, la ecuación <math display="inline">\left(\right. </math>[[#eq-16|16]]<math>\left. \right)</math> es equivalente a encontrar <math display="inline">\tilde{u}\in \tilde{\mathcal{W}}</math> tal que

<span id="eq-17"></span>
<span id="eq-18"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{L}\left(\tilde{u}\right) =f-\mathcal{L}\left(u_{h}\right)-\partial _{t} u_{h}+v_{h,{ort}},\qquad \mathrm{en}\;\Omega ^{e}\quad \forall v_{h,{ort}}\in \tilde{\mathcal{W}}^{\perp }, </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|-
| style="text-align: center;" | <math> \tilde{u}  =\tilde{u}_{esq},\qquad \mathrm{en}\;\partial \Omega ^{e}, </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}
|}

para <math display="inline">e=1,...,n_{el}</math> y para una cierta función <math display="inline">\tilde{u}_{esq}</math>, que llamamos el esqueleto de <math display="inline">u</math>. Es importante observar que la ecuación <math display="inline">\left(\right. </math>[[#eq-17|17]]<math>\left. \right)</math> se verifica para cualquier elemento <math display="inline">v_{h,{ort}}</math> ortogonal a <math display="inline">\tilde{\mathcal{W}}.</math> La presencia del elemento <math display="inline">v_{h,{ort}}</math> en la ecuación <math display="inline">\left(\right. </math>[[#eq-17|17]]<math>\left. \right)</math> garantiza que <math display="inline">\mathcal{L}\left(\tilde{u}\right)-\left[f-\mathcal{L}\left(u_{h}\right)\right]+\partial _{t} u_{h}</math> pertenezca a <math display="inline">\tilde{\mathcal{W}}^{\perp }</math>, que es lo que implica la ecuación <math display="inline">\left(\right. </math>[[#eq-16|16]]<math>\left. \right)</math> sin el primer término. En otras palabras, <math display="inline">v_{h,{ort}}</math> garantiza que <math display="inline">\tilde{u}</math> pertenezca al espacio de subescalas <math display="inline">\tilde{\mathcal{W}}</math> (y no a todo el espacio <math display="inline">\mathcal{W}</math>).

Las ecuaciones <math display="inline">\left(\right. </math>[[#eq-15|15]]<math>\left. \right),</math> <math display="inline">\left(\right. </math>[[#eq-17|17]]<math>\left. \right)</math> y <math display="inline">\left(\right. </math>[[#eq-18|18]]<math>\left. \right)</math> equivalen exactamente al sistema de ecuaciones <math display="inline">\left(\right. </math>[[#eq-10|10]]<math>\left. \right)</math> y <math display="inline">\left(\right. </math>[[#eq-11|11]]<math>\left. \right),</math> en donde hasta ahora no hemos hecho ninguna aproximación. El siguiente paso es hacer una selección de las funciones <math display="inline">\tilde{u}_{esq},</math> <math display="inline">v_{h,ort}</math> y una solución aproximada de la ecuación <math display="inline">\left(\right. </math>[[#eq-17|17]]<math>\left. \right).</math> La manera de escoger dichas funciones da lugar a las formulaciones ASGS y OSS que describimos más adelante. Los pasos comunes los describimos a continuación.

Para seleccionar <math display="inline">\tilde{u}_{esq}</math> vamos a usar el siguiente argumento. Para tener una aproximación discreta óptima de <math display="inline">u</math> que dé valores nodales exactos, se podría pedir que las subescalas se anulen en los contornos de los elementos. En problemas unidimensionales esto da condiciones de contorno homogéneas para el problema <math display="inline">\left(\right. </math>[[#eq-17|17]]<math>\left. \right)</math> y <math display="inline">\left(\right. </math>[[#eq-18|18]]<math>\left. \right)</math> en los extremos de cada elemento. Para los casos de más de una dimensión espacial, esto no es posible, puesto que la solución del problema continuo <math display="inline">u</math> tendría que ser polinómica en los contornos de los elementos para que pudiera coincidir con <math display="inline">u_h</math>. Sin embargo, imponer que <math display="inline">\tilde{u}_{esq}=0</math> puede pensarse como una aproximación.

Suponer que <math display="inline">\tilde{u}_{esq}=0</math> equivale a considerar el espacio de las subescalas <math display="inline">\tilde{\mathcal{W}}</math> como un espacio de funciones burbuja, es decir, un espacio de funciones que se anulan en los contornos de los elementos (véase <span id='citeF-19'></span><span id='citeF-6'></span>[[#cite-19|[19,6]]]). Con esta aproximación, los términos con las integrales en los contornos de los elementos de la ecuación <math display="inline">\left(\right. </math>[[#eq-15|15]]<math>\left. \right)</math> desaparecen, por lo cual el problema a ser resuelto consiste en encontrar para cada instante de tiempo <math display="inline">t</math> una función <math display="inline">u_h(\cdot ,t)\in \mathcal{W}_{h}</math> tal que

<span id="eq-19"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left(\partial _{t} u_{h},v_{h}\right)+B\left(u_{h},v_{h}\right)+\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\Omega ^{e}}\mathcal{L}^{*}\left(v_{h}\right)\tilde{u}=\left\langle f,v_{h}\right\rangle \qquad \forall v_{h}\in \mathcal{W}_{h}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

en donde todavía hace falta determinar <math display="inline">\tilde{u}</math>. Para ello, tenemos que resolver de forma aproximada el problema <math display="inline">\left(\right. </math>[[#eq-17|17]]<math>\left. \right)</math> y <math display="inline">\left(\right. </math>[[#eq-18|18]]<math>\left. \right)</math> en el espacio de las subescalas <math display="inline">\tilde{\mathcal{W}}</math> de funciones burbuja. Formalmente podemos escribir:

<span id="eq-20"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tilde{u}=\mathcal{L}^{-1}\left[f-\mathcal{L}\left(u_{h}\right)-\partial _{t}u_{h} + v_{h,ort}\right], </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}

en cada elemento, donde <math display="inline">\mathcal{L}^{-1}</math> es el operador inverso del operador <math display="inline">\mathcal{L}</math> con condiciones de Dirichlet homogéneas.

De la ecuación <math display="inline">\left(\right. </math>[[#eq-19|19]]<math>\left. \right)</math> vemos que solo se necesita la componente de <math display="inline">\tilde{u}</math> sobre el espacio <math display="inline">\mathcal{L}\left(\mathcal{W}_{h}\right)</math>, donde <math display="inline">\mathcal{L}\left(\mathcal{W}_{h}\right)</math> es el espacio de funciones de la forma <math display="inline">\mathcal{L}\left(v_{h}\right)</math>, con <math display="inline">v_{h}\in \mathcal{W}_{h}</math>. Esto sugiere aproximar <math display="inline">\left(\right. </math>[[#eq-20|20]]<math>\left. \right)</math> como

<span id="eq-21"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tilde{u}\approx \tau ^{e}\left[f-\mathcal{L}\left(u_{h}\right)-\partial _{t} u_{h} + v_{h,ort} \right]\quad \hbox{en}~\Omega ^e, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}

donde <math display="inline">\tau ^{e}</math> es el parámetro de estabilización y que se calcula dentro del dominio de cada elemento. Es una aproximación ''algebraica'' al operador <math display="inline">\mathcal{L}^{-1}</math>.

El diseño de <math display="inline">\tau ^{e}</math> es una de las piedras angulares en el desarrollo de los métodos de elementos finitos estabilizados. Algunas formulaciones algebraicas para el diseño del parámetro de estabilización para la ecuación escalar de CDR en 1D, basadas en el principio del máximo discreto, pueden encontrarse en <span id='citeF-10'></span><span id='citeF-11'></span><span id='citeF-9'></span>[[#cite-10|[10,11,9]]]. Un posible argumento para obtener <math display="inline">\tau ^e</math> es el análisis de Fourier que se presenta en <span id='citeF-17'></span><span id='citeF-20'></span><span id='citeF-21'></span>[[#cite-17|[17,20,21]]]. En este trabajo, para abordar la utilización de elementos finitos de alto orden, el diseño del <math display="inline">\tau </math> que adoptaremos es:

<span id="eq-22"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tau ^e=\left[c_{1}\frac{k}{\left(\frac{h^{e}}{p^{2}}\right)^{2}}+c_{2}\frac{|\boldsymbol{a}|}{\frac{h^{e}}{p}}+c_{3}s\right]^{-1}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|}

donde <math display="inline">c_{1},</math> <math display="inline">c_{2},</math> <math display="inline">c_{3}</math> son constantes algorítmicas,  que para los ejemplos numéricos hemos adoptado <math display="inline">c_{1}=12,</math> <math display="inline">c_{2}=2</math> y <math display="inline">c_{3}=1,</math> <math display="inline">k</math> es coeficiente de difusión, <math display="inline">|\mathbf{a}|</math> es la norma de la velocidad, <math display="inline">s</math> es el coeficiente de reacción, <math display="inline">h^{e}</math> es el diámetro del elemento y <math display="inline">p</math> es el orden polinomial de interpolación.

A falta de determinar la función <math display="inline">v_{h,ort}</math>, y por consiguiente el espacio de subescalas, el problema para la función <math display="inline">u_h</math> se obtiene reemplazando la ecuación ([[#eq-21|21]]) en ([[#eq-19|19]]) y usando la expresión ([[#eq-22|22]]).

====2.4.2 El método ASGS (Algebraic Sub-Grid Scales)====

La idea del método ASGS es tomar como espacio de subescalas el espacio de funciones que son residuos de funciones de elementos finitos, es decir, el espacio de funciones de la forma <math display="inline">f-\mathcal{L}\left(v_{h}\right)-\partial _{t} v_{h}</math>, con <math display="inline">v_h\in \mathcal{W}_h</math>. Esto nos lleva a que <math display="inline">v_{h,ort} = 0</math>.

El problema que finalmente obtenemos a partir de  ([[#eq-21|21]]) y ([[#eq-19|19]]) con <math display="inline">v_{h,ort} = 0</math> es:

<span id="eq-23"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left(\partial _{t}u_{h},v_{h}\right)+B\left(u_{h},v_{h}\right)+\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\Omega ^{e}}\left(-\mathcal{L}^{*}\left(v_{h}\right)\right)\tau ^{e}\left[\partial _{t}u_{h}+\mathcal{L}\left(u_{h}\right)-f\right]=\left\langle f,v_{h}\right\rangle \qquad \forall v_{h}\in \mathcal{W}_{h}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}

Esta ecuación variacional discreta, a la cual hay que añadirle las condiciones iniciales, la podemos ahora discretizar en el tiempo usando por ejemplo la regla trapezoidal descrita anteriormente. Para ello, habrá que sustituir <math display="inline">\partial _{t}u_{h}</math> por <math display="inline">\delta ^n_{t}u_{h}</math> y evaluar el resto de funciones dependientes del tiempo que aparecen en <math display="inline">t^{n+\theta }</math>.

El parámetro de estabilización <math display="inline">\tau ^{e}</math> para elementos finitos de alto orden está dado por la ecuación <math display="inline">\left(\right. </math>[[#eq-22|22]]<math>\left. \right)</math>, el operador <math display="inline">\mathcal{L}</math> por <math display="inline">\left(\right. </math>[[#eq-2|2]]<math>\left. \right)</math> y el operador <math display="inline">\mathcal{L}^{*}</math> por la ecuación <math display="inline">\left(\right. </math>[[#eq-13|13]]<math>\left. \right)</math>.

====2.4.3 El método OSS (Orthogonal Subscales Stabilization)====

En este caso, el espacio de subescalas no será el espacio de residuos de elementos finitos como en el caso anterior, sino que se toma como espacio <math display="inline">\mathcal{\tilde{W}}</math> el espacio ortogonal en el sentido de <math display="inline">L^{2}</math> al espacio de elementos finitos, es decir,

<span id="eq-24"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{\tilde{W}=\mathcal{W}}_{h}^{\perp }, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|}

que es una legítima elección que cumple con <math display="inline">\mathcal{W}=\mathcal{W}_{h}\oplus \mathcal{\tilde{W}}</math>.

Tomando en cuenta ([[#eq-24|24]]), de ([[#eq-17|17]]) se sigue que

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> v_{h,{ort}}\in \mathcal{W}_{h},\quad \tilde{u}\in \mathcal{W}_{h}^{\perp }, </math>
|}
|}

lo que significa que <math display="inline">v_{h,{ort}}</math> es una función del espacio de elementos finitos y por lo tanto numéricamente calculable, mientras que <math display="inline">\tilde{u}</math> estará en el espacio ortogonal a <math display="inline">\mathcal{W}_{h}</math>. A este espacio <math display="inline">\mathcal{W}_{h}^{\perp }</math> escogido para las subescalas <math display="inline">\tilde{u}</math> lo denominamos espacio de subescalas ortogonales.

Puesto que <math display="inline">\tilde{u}\in \mathcal{W}_{h}^{\perp }</math>, de ([[#eq-21|21]]) se tiene que

<span id="eq-25"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>0=\int _{\Omega }\tilde{u}v_{h}\approx \overset{n_{{el}}}{\underset{e=1}{\sum }}\int _{\Omega ^{e}}\tau ^{e}\left[f-\mathcal{L}\left(u_{h}\right)-\partial _{t} u_{h}\right]v_{h}+\overset{n_{{el}}}{\underset{e=1}{\sum }}\int _{\Omega ^{e}}\tau ^{e}\left(v_{h,{ort}}\right)v_{h}\qquad \forall v_{h}\in \mathcal{W}_{h}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|}

Si llamamos <math display="inline">{\Pi }_{\tau }</math> la proyección <math display="inline">L^{2}</math> sobre el espacio de elementos finitos pesada con <math display="inline">\tau ^{e}</math> dentro de cada elemento, de ([[#eq-25|25]]) se tiene que

<span id="eq-26"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>v_{h,{ort}}=-{\Pi }_{\tau }\left[f-\mathcal{L}\left(u_{h}\right)-\partial _{t}u_{h}\right]. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|}

Reemplazando ([[#eq-26|26]]) en ([[#eq-21|21]]) tenemos

<span id="eq-27"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tilde{u}\approx \tau ^{e}\left[f-\mathcal{L}\left(u_{h}\right)-\partial _{t} u_{h}-{\Pi }_{\tau }\left[f-\mathcal{L}\left(u_{h}\right)-\partial _{t}u_{h}\right]\right]. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
|}

Sea ahora

<span id="eq-28"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\Pi }_{\tau }^{\perp }={I}-{{\Pi }_{\tau }}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|}

donde <math display="inline">{I}</math> es la identidad en <math display="inline">\mathcal{W}</math> y <math display="inline">{\Pi }_{\tau }^{\perp }</math> es la proyección sobre el espacio ortogonal al espacio de elementos finitos. La ecuación ([[#eq-27|27]]) se puede escribir entonces como

<span id="eq-29"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tilde{u}\approx \tau ^{e}{\Pi }_{\tau }^{\perp }\left[f-\mathcal{L}\left(u_{h}\right)-\partial _{t}u_{h}\right]. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
|}

Podemos reemplazar esta expresión en ([[#eq-15|15]]) y usar el hecho de que las integrales sobre los contornos de los elementos son cero, con lo cual obtenemos la formulación estabilizada de elementos finitos OSS que resuelve la ecuación escalar de CDR transitoria, la cual consiste en encontrar <math display="inline">u_{h}(\cdot ,t)\in \mathcal{W}_{h}</math> en cada instante de tiempo <math display="inline">t</math> tal que

<span id="eq-30"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left(\partial _{t}u_{h},v_{h}\right)+B\left(u_{h},v_{h}\right)+\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\Omega ^{e}}\left(-\mathcal{L}^{*}\left(v_{h}\right)\right)\tau ^{e}{\Pi }_{\tau }^{\perp }\left[\partial _{t} u_{h}+\mathcal{L}\left(u_{h}\right)-f\right]=\left\langle f,v_{h}\right\rangle \qquad \forall v_{h}\in \mathcal{W}_{h}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
|}

Podemos aplicar aquí los mismos comentarios que para <math display="inline">\left(\right. </math>[[#eq-23|23]]<math>\left. \right)</math>. Sin embargo, hay una consideración importante a tener en cuenta, y es que puesto que <math display="inline">\partial _{t} u_{h}</math> es una función que pertenece al espacio de elementos finitos <math display="inline">\mathcal{W}_{h}</math>, su componente en el espacio ortogonal es nula, esto es, la ecuación variacional discreta anterior se reduce a

<span id="eq-31"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left(\partial _{t} u_{h},v_{h}\right)+B\left(u_{h},v_{h}\right)+\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\Omega ^{e}}\left(-\mathcal{L}^{*}\left(v_{h}\right)\right)\tau ^{e}{\Pi }_{\tau }^{\perp }\left[\mathcal{L}\left(u_{h}\right)-f\right]=\left\langle f,v_{h}\right\rangle \qquad \forall v_{h}\in \mathcal{W}_{h}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|}

====2.4.4 Propiedades de estabilidad y convergencia====

Aunque los métodos ASGS y OSS difieren en la definición del espacio de subescalas, ambos tienen las mismas propiedades de estabilidad y convergencia. No es nuestro objetivo presentar el análisis numérico de estos métodos, el cual puede encontrarse en <span id='citeF-22'></span>[[#cite-22|[22]]] (véase también <span id='citeF-23'></span><span id='citeF-24'></span>[[#cite-23|[23,24]]] y <span id='citeF-25'></span>[[#cite-25|[25]]] para la teoría general), sino mostrar en qué sentido dichos métodos son estables y convergen, justificar la expresión del parámetro de estabilización (y en particular de su dependencia en <math display="inline">p</math>) e indicar cómo debe ser la integración numérica para preservar la convergencia asociada a la integración exacta. Explotaremos este último punto en el apartado siguiente.

Para simplificar la exposición, consideremos el problema estacionario, con coeficientes constantes, <math display="inline">s = 0</math> y malla de elementos finitos uniforme, de tamaño <math display="inline">h</math>. En este caso, el parámetro de estabilización dado por ([[#eq-22|22]]) será el mismo en todos los elementos.

Los métodos ASGS y OSS se puede demostrar que son estables y convergentes en la norma:

<span id="eq-32"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\vert \!\vert \!\vert }v_{h}{\vert \!\vert \!\vert }^2 :=  k \Vert \nabla v_{h}\Vert ^2 +  \tau \Vert \mathbf{a}\cdot \nabla v_{h}\Vert ^2, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

donde <math display="inline">\Vert \cdot \Vert </math> es la norma de <math display="inline">L^2(\Omega )</math>. El punto clave para demostrar esta estabilidad es hacer uso del estimador inverso (ver por ejemplo <span id='citeF-25'></span>[[#cite-25|[25]]])

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\Vert \nabla v_{h}\Vert \leq C_{inv}\frac{p^2}{h}\Vert v_{h}\Vert },  </math>
|}
|}

válido para cualquier función de elementos finitos <math display="inline">v_h</math>, siendo <math display="inline">C_{inv}</math> una constante que depende de la forma de la malla, pero no de su tamaño <math display="inline">h</math> ni del orden de interpolación <math display="inline">p</math> . La potencia 2 de <math display="inline">p^2</math> es precisamente quien obliga a dividir <math display="inline">h</math> por <math display="inline">p^2</math> en la contribución del término difusivo en el parámetro de estabilización dado por&nbsp;([[#eq-22|22]]). Con ello tenemos garantizada la estabilidad.

Para analizar la convergencia es necesario hacer uso de las propiedades de interpolación. Si <math display="inline">\Vert \cdot \Vert _{H^r}</math> es la norma del espacio de Sobolev <math display="inline">H^r(\Omega )</math>, la solución <math display="inline">u</math> del problema continuo es suficientemente regular y <math display="inline">\tilde{u}_h</math> es su mejor aproximación en el espacio de elementos finitos, se tiene que

<span id="eq-33"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\Vert u - \tilde{u}_h \Vert _{H^r} \leq C \frac{h^{p+1-r}}{p^{p+1-r}} \Vert u\Vert _{H^{p+1}}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}

donde los casos de interés corresponden a <math display="inline">r=0</math> y <math display="inline">r=1</math>. Aquí y en lo que sigue, <math display="inline">C</math> es una constante que no depende de la malla. Con este estimador de interpolación, es posible demostrar que la función de error de las formulaciones que consideramos, dada por

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>E(h)  =  k^{1/2} \Vert \nabla u - \nabla  \tilde{u}_h\Vert + \tau ^{1/2} \sum _K  \Vert -k\Delta (u - \tilde{u}_h) + \mathbf{a}\cdot \nabla (u - \tilde{u}_h)\Vert _K  + \tau ^{-1/2}\Vert u - \tilde{u}_h\Vert  </math>
|}
|}

se comporta como

<span id="eq-34"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>E(h)  = \left[k + \frac{\vert \mathbf{a}\vert h}{p} + \min \left(p^2, \frac{\vert \mathbf{a}\vert ^2 h^2}{p^2 k^2}\right)k\right]^{1/2}   \frac{h^p}{p^p} \Vert u\Vert _{H^{p+1}}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|}

Que <math display="inline">E(h)</math> sea la función de error de los métodos ASGS y OSS en la norma ([[#eq-32|32]]) quiere decir que

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\vert \!\vert \!\vert }u - u_{h}{\vert \!\vert \!\vert }\leq C E(h).  </math>
|}
|}

En virtud de la expresión ([[#eq-32|32]]) de la norma de trabajo y de la función de error ([[#eq-34|34]]), se observa que el estimador de error presentado es óptimo tanto para advección dominante como para difusión dominante. Hay que notar sin embargo que si se considera la convergencia en <math display="inline">p</math>, hay un rango de parámetros en los que se pierde la optimalidad (véase <span id='citeF-22'></span>[[#cite-22|[22]]]).

El objetivo de haber presentado los resultados de estabilidad y convergencia en este trabajo es doble. Por un lado, justifica la validez de la expresión del parámetro de estabilización ([[#eq-22|22]]) para elementos de alto orden, que es el caso que nos ocupa. Por otro lado, nos permite observar que si la integración es numérica y no exacta (como obviamente siempre es el caso), hay un límite de exactitud mínimo para no perder orden de convergencia de la aproximación. En el apartado siguiente presentaremos una modificación del elemento triangular de cuarto orden con una cuadratura numérica asociada. En este caso, <math display="inline">p = 4</math>, con lo que de acuerdo con ([[#eq-34|34]]) la cuadratura debe ser como mínimo de cuarto orden en el caso de difusión dominante y de orden <math display="inline">9/2</math> para convección dominante. De hecho, en presencia de término reactivo (<math display="inline">s>0</math>) y con reacción dominante se obtiene que la cuadratura mínima debe ser de quinto orden. Este hecho lo usamos a continuación.

==3. Elemento finito triangular modificado de cuarto orden==

===3.1 Introducción===

En el método OSS, para calcular la proyección ortogonal del residuo en la ecuación <math display="inline">\left(\right. </math>[[#eq-31|31]]<math>\left. \right)</math>, lo que hacemos es calcular la proyección del residuo en el espacio de elementos finitos. Suponiendo para simplificar que los parámetros de estabilización son iguales para todos los elementos, debemos resolver el sistema

<span id="eq-35"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _\Omega v_h  \Pi _\tau = \int _\Omega v_h \left[\mathcal{L}(u_h) - f\right], \quad \forall v_h \in \mathcal{W}_h, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|}

siendo <math display="inline">\Pi _\tau = \Pi _\tau ( \mathcal{L}(u_h) - f)</math>. Cuando expresamos <math display="inline">\Pi _\tau </math> en función de los valores nodales <math display="inline">\mathbf{\Pi }</math> y tomamos <math display="inline">v_h</math> como las funciones de forma, ([[#eq-35|35]]) da lugar al sistema lineal

<span id="eq-36"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{M} \mathbf{\Pi } = \mathbf{R}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
|}

donde <math display="inline">\mathbf{R}</math> es conocido y <math display="inline">\mathbf{M}</math> es la matriz de masa de la partición de elementos finitos. Esta matriz de masa es en principio llena, por lo que es necesario resolver el sistema ([[#eq-36|36]]). Sin embargo, es bien sabido que <math display="inline">\mathbf{M}</math> se puede aproximar por una matriz diagonal <math display="inline">\mathbf{M}_d</math> usando una regla de integración numérica con los puntos de integración en los nodos (cuadratura cerrada) en vez de la integral exacta que aparece en ([[#eq-35|35]]). Con ello la resolución del sistema que aproxima a ([[#eq-36|36]]) es trivial. El punto clave es que la exactitud de la regla de integración numérica debe preservar la de la aproximación de elementos finitos que se tendría con la integración exacta.

En ocasiones las matrices de masa diagonales obtenidas a través de una cuadratura cerrada presentan en su diagonal valores nulos, lo cual hace que no sean invertibles e imposibilita su uso. En elementos lineales como el <math display="inline">P_1</math> (simplicial, es decir, triangular para <math display="inline">d=2</math> y tetraédrico para <math display="inline">d=3</math>), esta matriz se obtiene directamente utilizando integración cerrada, con todos sus pesos distintos de cero. Lo mismo sucede para elementos cúbicos, mientras que los elementos cuadráticos tienen que subdividirse en elementos lineales y usar cuadradura nodal para la división resultante; con este procedimiento puede obtenerse una matriz de masa diagonal que, si bien reduciría el orden de convergencia en caso de que se reemplazara por ella la matriz de masa original, puede usarse como precondicionador para resolver el sistema ([[#eq-36|36]]).

Atención especial merece el elemento de cuarto orden <math display="inline">P_4</math> simplicial. Veremos a continuación que es posible modificar el elemento convencional, en el sentido de alterar la posición de los nodos centrales, y plantear una regla de integración numérica con los puntos de integración en los nodos y los pesos distintos de cero, lo cual nos permitirá aproximar la matriz de masa por una matriz diagonal invertible. Además, dicha regla de integración veremos que es capaz de integrar de manera exacta hasta un polinomio completo de quinto grado, con lo que de acuerdo con la discusión del apartado anterior no se disminuirá el orden de aproximación de la formulación final.

Aunque es posible aplicar la modificación que proponemos tanto para <math display="inline">d= 2</math> como para <math display="inline">d= 3</math>, haremos los desarrollos en el primer caso. Plantearemos el elemento <math display="inline">P_4</math> modificado considerando móviles los tres nodos interiores del elemento <math display="inline">P_4</math> original y parametrizando la posición de los mismos con una variable cuyo valor será incógnita, al igual que los pesos de la regla de integración nodal asociada; la posición de los 12 nodos restantes es la misma que la del elemento <math display="inline">P_4</math> original.

===3.2 <span id='lb-3.2'></span>Descripción del elemento P₄ modificado===

En la Figura [[#img-1|1]] se representa el elemento en el dominio de referencia, dado por el triángulo <math display="inline">0\leq x \leq 1</math>, <math display="inline">0\leq y \leq 1</math>, <math display="inline">x + y \leq 1</math>. Las posiciones de los nodos sobre los lados es la del elemento estándar, pero hemos introducido una variable <math display="inline">z</math> para parametrizar  la posición de los tres nodos interiores, obteniéndose así el elemento que llamaremos <math display="inline">P_4</math> modificado, con los tres nodos interiores móviles en función del parámetro <math display="inline">z</math>. Obsérvese que con <math display="inline">z=1/4</math> tendríamos el elemento <math display="inline">P_4</math> original.

<div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 50%;"
|-
|[[Image:Review_640781029593_1838_test-FIG_1.png|526px|Elemento P4-modificado con los 3 nodos interiores móviles]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figura 1'''. Elemento P4-modificado con los 3 nodos interiores móviles
|}

===3.3 Integral exacta===

Consideremos un polinomio de grado <math display="inline">n</math> de la forma:

<span id="eq-37"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p\left(x,y\right)=\overset{n+1}{\underset{i=1}{\sum }}\overset{i}{\underset{j=1}{\sum }}a_{ij}x^{i-j}y^{j-1}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}

donde los coeficientes <math display="inline">a_{ij}</math> son constantes. La integral exacta <math display="inline">I_{E}</math> de <math display="inline">p(x,y)</math> sobre la región triangular de la Figura [[#img-1|1]] está dada por:

<span id="eq-38"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>I_{E}=\sum _{i=1}^{n+1}\sum _{j=1}^{i}a_{ij}\int _{0}^{1}\int _{0}^{1-x}x^{i-j}y^{j-1}{d}y{d}x. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
|}

Sea

<span id="eq-39"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\alpha _{ij}=\int _{0}^{1}\int _{0}^{1-x}x^{i-j}y^{j-1}{d}y{d}x\quad \,\,\,\,\forall i=1,2,...,n+1,\,\,\forall j=1,2,...,i.\qquad  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
|}

Reemplazando ([[#eq-39|39]]) en ([[#eq-38|38]]) tenemos:

<span id="eq-40"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>I_{E}=\sum _{i=1}^{n\hbox{+1}}\sum _{j=1}^{i}a_{ij}\alpha _{ij}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|}

===3.4 Integral numérica===

La integral exacta <math display="inline">I_{E}</math> vamos a aproximarla mediante una cuadratura cerrada sobre una región triangular con 15 puntos distribuidos de la siguiente manera: 3 puntos (1, 2, 3) en los vértices del triángulo, 9 puntos numerados del 4 al 12 distribuidos simétricamente en los lados del triángulo y los 3 puntos móviles interiores restantes (13, 14, 15) distribuidos simétricamente como se ve en la Figura [[#img-2|2]].

<div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 65%;"
|-
|[[Image:Review_640781029593_5330_test-FIG_2.png|526px|Posición de los 15 puntos de integración numérica]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figura 2'''. Posición de los 15 puntos de integración numérica
|}

Los pesos <math display="inline">w_{1},w_{2},w_{3}</math> y <math display="inline">w_{4}</math> están asociados a los puntos detallados en la Figura [[#img-3|3]] para los 15 puntos de integración. Estos pesos, junto con el parámetro <math display="inline">z</math>, son las incógnitas a determinar.

<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 65%;"
|-
|[[Image:Review_640781029593_1721_test-FIG_3.png|526px|Pesos para los 15 puntos de integración numérica]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figura 3'''. Pesos para los 15 puntos de integración numérica
|}

Con las posiciones y los pesos de los 15 puntos de integración que se detallan en las Figuras [[#img-2|2]] y [[#img-3|3]], planteamos la integral numérica, <math display="inline">I_{N}</math>, como sigue:

<span id="eq-41"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>I_{N}  =  \sum _{k=1}^{3}w_{1}p\left(x_{k},y_{k}\right)+\sum _{k=1}^{3}w_{2}p\left(x_{3k+2},y_{3k+2}\right)+\sum _{k=1}^{3}w_{3}p\left(x_{3k+1},y_{3k+1}\right)</math>
|-
| style="text-align: center;" | <math>     +\sum _{k=1}^{3}w_{3}p\left(x_{3k+3},y_{3k+3}\right)+\sum _{k=1}^{3}w_{4}p\left(x_{12+k},y_{12+k}\right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
|}

De la ecuación ([[#eq-41|41]]), reemplazando <math display="inline">p(x,y)</math> dado en la ecuación ([[#eq-37|37]]) y agrupando términos, obtenemos:

<span id="eq-42"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>I_{N}  =  \sum _{k=1}^{3}(\sum _{i=1}^{n+1}\sum _{j=1}^{i}a_{ij}(w_{1}x_{k}^{i-j} y_{k}^{j-1} +w_{2}x_{3k+2}^{i-j} y_{3k+2}^{j-1}+w_{3}x_{3k+1}^{i-j}y_{3k+1}^{j-1}</math>
|-
| style="text-align: center;" | <math>    \quad + w_{3}x_{3k+3}^{i-j}y_{3k+3}^{j-1}+w_{4}x_{12+k}^{i-j}y_{12+k}^{j-1}))</math>
|-
| style="text-align: center;" | <math>  =  \sum _{i=1}^{n+1}\sum _{j=1}^{i}a_{ij}\sum _{k=1}^{3}(w_{1}x_{k}^{i-j}y_{k}^{j-1}+w_{2}x_{3k+2}^{i-j}y_{3k+2}^{j-1}+w_{3}x_{3k+1}^{i-j}y_{3k+1}^{j-1}</math>
|-
| style="text-align: center;" | <math>   \quad +w_{3}x_{3k+3}^{i-j}y_{3k+3}^{j-1}+w_{4}x_{12+k}^{i-j}y_{12+k}^{j-1}). </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
|}

Puesto queremos que la integral exacta y la integral numérica sean iguales, esto es <math display="inline">I_{E}=I_{N}</math>, de las ecuaciones ([[#eq-40|40]]) y ([[#eq-42|42]]) se tiene que:

<span id="eq-43"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\alpha _{ij}  =  \sum _{k=1}^{3}(w_{1}x_{k}^{i-j}y_{k}^{j-1}+w_{2}x_{3k+2}^{i-j}y_{3k+2}^{j-1}+w_{3}x_{3k+1}^{i-j}y_{3k+1}^{j-1}</math>
|-
| style="text-align: center;" | <math>     \quad +w_{3}x_{3k+3}^{i-j}y_{3k+3}^{j-1}+w_{4}x_{12+k}^{i-j}y_{12+k}^{j-1}), </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
|}

para todo <math display="inline">i=1,2,...,n+1</math>, <math display="inline">j=1,2,...,i</math>, donde <math display="inline">n</math> es el grado del polinomio completo cuya integral exacta sobre una región triangular se quiere calcular numéricamente.

Las ecuaciones ([[#eq-39|39]]) y ([[#eq-43|43]]) nos permiten obtener un sistema de ecuaciones no lineal para determinar <math display="inline">w_{1}</math>, <math display="inline">w_{2}</math>, <math display="inline">w_{3}</math>, <math display="inline">w_4</math> y <math display="inline">z</math>. Recordemos que los puntos <math display="inline">(x_{k},y_{k}),</math> <math display="inline">(x_{3k+1},y_{3k+1}),</math> <math display="inline">(x_{3k+2},y_{3k+2}),</math> <math display="inline">(x_{3k+3},y_{3k+3})</math>, <math display="inline">k=1,2,3</math>, son las coordenadas de los 12 puntos ubicados en los vértices y lados del triángulo respectivamente, y los puntos <math display="inline">(x_{12+k},y_{12+k})</math>, <math display="inline"> k=1,2,3</math>, son las coordenadas de los 3 puntos interiores cuya ubicación está parametrizada con la variable <math display="inline">z</math>.

===3.5 Cuadratura cerrada para integración exacta de un polinomio de cuarto grado===

Utilizando las ecuaciones ([[#eq-39|39]]) y ([[#eq-43|43]]) para <math display="inline">n=4,</math> fácilmente se puede plantear el sistema no lineal de 15 ecuaciones para integrar exactamente un polinomio completo de cuarto grado. De las 15 ecuaciones, obtenemos tres ecuaciones linealmente independientes en términos del parámetro <math display="inline">z</math>, que son:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{w_{2}}{16}+\frac{15}{128}w_{3}+(-9z^{4}+13z^{3}-6z^{2}+z)w_{4}  =\frac{1}{120},</math>
|-
| style="text-align: center;" | <math> \frac{w_{2}}{16}+\frac{9}{128}w_{3}+(9z^{4}-8z^{3}+2z^{2})w_{4}  =\frac{1}{180},</math>
|-
| style="text-align: center;" | <math> \frac{w_{2}}{8}+\frac{3}{16}w_{3}+(3z^{3}-3z^{2}+z)w_{4}  =\frac{1}{60}. </math>
|}
|}

Resolviendo este sistema de ecuaciones en términos del parámetro <math display="inline">z</math>, se tiene:

<span id="eq-44"></span>
<span id="eq-45"></span>
<span id="eq-46"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>w_{2} =\frac{72z^{3}-77z^{2}+27z-3}{45z(2z-1)},</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
|-
| style="text-align: center;" | <math> w_{3} =\frac{-8(18z^{3}-23z^{2}+9z-1}{135z(2z-1)},</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
|-
| style="text-align: center;" | <math> w_{4} =\frac{-1}{360z^{2}(2z-1)}. </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
|}
|}

Para obtener la ecuación parametrizada para <math display="inline">w_{1}</math>, podemos usar la primera ecuación del sistema de 15 ecuaciones mencionado anteriormente. Reemplazando <math display="inline">i=1,</math> <math display="inline">j=1,</math> en las ecuaciones ([[#eq-39|39]]) y ([[#eq-43|43]]), tenemos respectivamente:

<span id="eq-47"></span>
<span id="eq-48"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\alpha _{11}  =\int _{0}^{1}\int _{0}^{1-x}{d}y {d}x = \frac{1}{2},</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (47)
|-
| style="text-align: center;" | <math> \alpha _{11}   =3w_{1}+3w_{2}+6w_{3}+3w_{4}. </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (48)
|}
|}

Igualando las ecuaciones ([[#eq-47|47]]) y ([[#eq-48|48]])  y reemplazando las ecuaciones ([[#eq-44|44]]), ([[#eq-45|45]]) y ([[#eq-46|46]]), tenemos:

<span id="eq-49"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>w_{1}=\frac{288z^{3}-224z^{2}+50z-3}{1080z^{2}}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (49)
|}

Las ecuaciones parametrizadas ([[#eq-44|44]]), ([[#eq-45|45]]), ([[#eq-46|46]]) y ([[#eq-49|49]]), verifican las 11 ecuaciones restantes del sistema no lineal de 15 ecuaciones planteadas para <math display="inline">n=4.</math> Por lo tanto, con las ecuaciones parametrizadas para <math display="inline">w_{1},w_{2},w_{3}</math> y <math display="inline">w_{4}</math> hemos obtenido un conjunto infinito de soluciones de parámetro <math display="inline">z</math>. Por ejemplo para <math display="inline">z=\frac{1}{6}</math> tenemos: <math display="inline">w_{1}=\frac{2}{135},\,w_{2}=\frac{11}{180},\,w_{3}=\frac{-4}{135},\,w_{4}=\frac{3}{20}</math>, donde observamos que todos los pesos son distintos de cero. Como observación, notemos que para <math display="inline">z=\frac{1}{4}</math>, nos da, <math display="inline">w_{1}=0,\,w_{2}=-\frac{1}{90},\,w_{3}=\frac{2}{45},\,w_{4}=\frac{4}{45},</math> que es la cuadratura cerrada con los puntos de integración numérica coincidiendo con los nodos del elemento <math display="inline">P_4</math> original, donde no todos los pesos son distintos de cero.

===3.6 Cuadratura cerrada para integración exacta de un polinomio de quinto grado===

En el apartado anterior hemos obtenido un conjunto infinito de cuadraturas cerradas que integran en forma exacta hasta un polinomio completo de cuarto grado sobre una región triangular con puntos de integración coincidentes con los nodos del elemento <math display="inline">P_4</math> modificado, con sus tres nodos interiores móviles dependientes del parámetro <math display="inline">z</math>, y con valores de <math display="inline">z</math> para los cuales todos los pesos son distintos de cero. Ahora nos proponemos obtener de este conjunto infinito de cuadraturas, al menos una cuadratura que integre exactamente un polinomio completo de quinto grado, también con todos sus pesos distintos de cero.

Utilizando las ecuaciones ([[#eq-39|39]]) y ([[#eq-43|43]]) para <math display="inline">n=5,</math> obtenemos 6 ecuaciones no lineales adicionales al sistema no lineal anterior de 15 ecuaciones. De estas 6 ecuaciones adicionales, solo 1 es linealmente independiente con el sistema anterior, la misma que se obtiene reemplazando <math display="inline">i=6</math> y <math display="inline">j=1</math> en las ecuaciones ([[#eq-39|39]]) y ([[#eq-43|43]]):

<span id="eq-50"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>w_{1}+\frac{w_{2}}{16}+\frac{61}{128}w_{3}+(2z^{5}+(1-2z)^{5}w_{4})=\frac{1}{42}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (50)
|}

Reemplazando en la ecuación ([[#eq-50|50]]) las soluciones parametrizadas de <math display="inline">w_{1},w_{2},w_{3}</math> y <math display="inline">w_{4}</math> dadas por las ecuaciones ([[#eq-49|49]]), ([[#eq-44|44]]), ([[#eq-45|45]]) y ([[#eq-46|46]]) respectivamente, se obtiene la ecuación:

<span id="eq-51"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{z^{2}}{2}-\frac{z}{3}+\frac{1}{21}=0, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (51)
|}

cuya solución es:

<span id="eq-52"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>z=\frac{7+\sqrt{7}}{21},\qquad z=\frac{7-\sqrt{7}}{21}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (52)
|}

Reemplazando las soluciones parametrizadas de <math display="inline">w_{1},w_{2},w_{3}</math> y <math display="inline">w_{4}</math> en las 5 ecuaciones restantes, es decir, para <math display="inline">i=6</math> y <math display="inline">j=2,3,4,5,6</math> en la ecuaciones ([[#eq-39|39]]) y ([[#eq-43|43]]), se obtiene la misma ecuación ([[#eq-51|51]]), con lo cual la solución es única para todas las ecuaciones del sistema.

De las dos soluciones de <math display="inline">z</math>, vamos a encontrar la cuadratura cerrada para:

<span id="eq-53"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>z=\frac{7-\sqrt{7}}{21}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (53)
|}

Reemplazando el valor <math display="inline">z</math> de la ecuación ([[#eq-53|53]]) en las ecuaciones parametrizadas ([[#eq-49|49]]), ([[#eq-44|44]]), ([[#eq-45|45]]) y ([[#eq-46|46]]) se obtienen los pesos <math display="inline">w_{1},w_{2},w_{3}</math> y <math display="inline">w_{4}</math>, todos distintos de cero:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>w_{1}=\frac{11\sqrt{7}}{15120}+\frac{1}{216},\qquad \qquad \qquad  w_{2}=\frac{11\sqrt{7}}{630}-\frac{1}{30},</math>
|-
| style="text-align: center;" | <math>  w_{3}=\frac{4}{135}-\frac{4\sqrt{7}}{945},\qquad \qquad \qquad  w_{4}=\frac{49}{360}-\frac{7\sqrt{7}}{720}. </math>
|}
|}

Con todo esto, hemos obtenido un elemento <math display="inline">P_4</math> modificado, con las posiciones de los nodos interiores dadas por <math display="inline">z=\frac{7-\sqrt{7}}{21}</math>, y  cuya cuadratura cerrada sobre una región triangular con puntos de integración coincidentes con los nodos del elemento <math display="inline">P_4</math> modificado integra de forma exacta hasta un polinomio completo de quinto grado, siendo todos los pesos distintos de cero.

La matriz de masa aproximada que se obtiene con la cuadratura cerrada sobre el elemento <math display="inline">P_4</math> modificado es una matriz diagonal con todos sus elementos distintos de cero. Este elemento <math display="inline">P_4</math> modificado es el que utilizamos en todos nuestros ejemplos numéricos  para el cálculo de la proyección ortogonal del residuo en el método OSS.

En la Tabla [[#table-1|1]] presentamos la cuadratura cerrada del elemento <math display="inline">P_4</math> modificado que acabamos de obtener.

<math>a=\frac{11\sqrt{7}}{15120}+\frac{1}{216},\qquad </math>

<math>b=\frac{11\sqrt{7}}{630}-\frac{1}{30},</math>  

<math>c=\frac{4}{135}-\frac{4\sqrt{7}}{945},\qquad \quad \,d=\frac{49}{360}-\frac{7\sqrt{7}}{720}.</math>

<div class="center" style="font-size: 75%;">'''Tabla 1'''. Cuadratura cerrada del elemento <math>P_4</math> modificado</div>

{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;font-size:85%;"
|- style="border-top: 1px solid;"
| style="border-top: 2px solid;border-bottom: 2px solid;" | '''Puntos '''
| style="border-top: 2px solid;border-bottom: 2px solid;" | <math>x</math>
| style="border-top: 2px solid;border-bottom: 2px solid;"  | <math>y</math>
| style="border-top: 2px solid;border-bottom: 2px solid;"  | '''Pesos'''
|- style="border-top: 1px solid;"
| style="border: 1px solid;"  |   1 
| style="border: 1px solid;"  | 0 
| style="border: 1px solid;"  | 0 
| style="border: 1px solid;"  | <math>a</math>
|- style="border-top: 1px solid;"
| style="border: 1px solid;"  |  2 
| style="border: 1px solid;"  | 1 
| style="border: 1px solid;"  | 0 
| style="border: 1px solid;"  | <math>a</math>
|- style="border-top: 1px solid;"
| style="border: 1px solid;" |  3 
| style="border: 1px solid;" | 0 
| style="border: 1px solid;"  | 1 
| style="border: 1px solid;"  | <math>a</math>
|- style="border-top: 1px solid;"
| style="border: 1px solid;"  |  4 
| style="border: 1px solid;"  | 1/4 
| style="border: 1px solid;" | 0 
| style="border: 1px solid;"  | <math>c</math>
|- style="border-top: 1px solid;"
| style="border: 1px solid;"  |  5 
| style="border: 1px solid;"  | 1/2 
| style="border: 1px solid;"  | 0 
| style="border: 1px solid;"  | <math>b</math>
|- style="border-top: 1px solid;"
| style="border: 1px solid;"  |  6 
| style="border: 1px solid;"  | 3/4 
| style="border: 1px solid;"  | 0 
| style="border: 1px solid;"  | <math>c</math>
|- style="border-top: 1px solid;"
| style="border: 1px solid;" |  7 
| style="border: 1px solid;"  | 3/4 
| style="border: 1px solid;"  | 1/4 
| style="border: 1px solid;"  | <math>c</math>
|- style="border-top: 1px solid;"
| style="border: 1px solid;"  |  8 
| style="border: 1px solid;"  | 1/2 
| style="border: 1px solid;"  | 1/2 
| style="border: 1px solid;" | <math>b</math>
|- style="border-top: 1px solid;"
| style="border: 1px solid;"  |  9 
| style="border: 1px solid;"  | 1/4 
| style="border: 1px solid;" | 3/4 
| style="border: 1px solid;"  | <math>c</math>
|- style="border-top: 1px solid;"
| style="border: 1px solid;"  |  10 
| style="border: 1px solid;"  | 0 
| style="border: 1px solid;"  | 3/4 
| style="border: 1px solid;"  | <math>c</math>
|- style="border-top: 1px solid;"
| style="border: 1px solid;"  |  11 
| style="border: 1px solid;"  | 0 
| style="border: 1px solid;" | 1/2 
| style="border: 1px solid;"  | <math>b</math>
|- style="border-top: 1px solid;"
| style="border: 1px solid;" |  12 
| style="border: 1px solid;" | 0 
| style="border: 1px solid;"  | 1/4 
| style="border: 1px solid;"  | <math>c</math>
|- style="border-top: 1px solid;"
| style="border: 1px solid;"  |  13 
| style="border: 1px solid;"  | <math>\frac{7-\sqrt{7}}{21}</math>
| style="border: 1px solid;"  | <math>\frac{7-\sqrt{7}}{21}</math>
| style="border: 1px solid;"  | <math>d</math>
|- style="border-top: 1px solid;"
| style="border: 1px solid;" |  14 
| style="border: 1px solid;"  | <math>\frac{7+2\sqrt{7}}{21}</math>
| style="border: 1px solid;"  | <math>\frac{7-\sqrt{7}}{21}</math>
| colspan='1' style="border: 1px solid;"  | <math>d</math>
|- style="border-top: 1px solid;border-bottom: 1px solid;"
| style="border: 1px solid;" |  15 
| style="border: 1px solid;"  | <math>\frac{7-\sqrt{7}}{21}</math>
| style="border: 1px solid;"  | <math>\frac{7+2\sqrt{7}}{21}</math>
| style="border: 1px solid;"  | <math>d</math>
|}

==4. Resultados de los experimentos numéricos==

===4.1 Pruebas de convergencia en malla===

Las pruebas de convergencia en malla que hemos llevado a cabo consisten en calcular el error en norma <math display="inline">L^{2}</math>, es decir, la norma <math display="inline">L^2</math> de la diferencia entre entre la solución exacta y su aproximación numérica. Considerando solamente refinamiento en <math display="inline">h</math>, si el error del método se comporta como el error de interpolación, situación que puede considerarse óptima, este error debe comportarse como

<span id="eq-54"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>e = \left\Vert u-u_{h}\right\Vert _{L^{2}}\leq C h^{p+1}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (54)
|}

que corresponde al mismo comportamiento que el error de interpolación ([[#eq-33|33]]) con <math display="inline">r= 0</math> (y absorbiendo en la constante <math display="inline">C</math> la dependencia con <math display="inline">p</math>).  Consideraremos los casos <math display="inline">1 \leq p\leq 4</math>.

Si graficamos la ecuación [[#eq-54|(54)]] (en el caso de la igualdad) en el plano, con <math display="inline">\log (e)</math> en las ordenadas y <math display="inline">\log (h)</math> en las abscisas, tenemos una línea recta en la que <math display="inline">p+1</math> es la pendiente teórica óptima de convergencia. Las pruebas de convergencia en malla consisten entonces en verificar que la pendiente calculada sea precisamente <math display="inline">p+1</math>. Para ello, seleccionamos como solución exacta <math display="inline">u</math> a la función polinómica espacio-temporal

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> u\left(x,y,t\right)=x^{6}y^{6}\left(1-x\right)^{6}\left(1-y\right)^{6}t, </math>
|}
|}

e imponemos el término fuente <math display="inline">f</math> en la ecuación escalar de CDR transitoria ([[#eq-1|1]]) para que esta sea precisamente la solución exacta del problema:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> f=\partial _{t}u\left(x,y,t\right)-k\partial _{i}\partial _{i}u\left(x,y,t\right)+a_{i}\partial _{i}u\left(x,y,t\right)+su\left(x,y,t\right), </math>
|}
|}

donde las constantes de difusión <math display="inline">k</math>, las componentes de la velocidad de convección <math display="inline">a_{i}</math> y la reacción <math display="inline">s</math> están dadas en cada uno de los ejemplos, tomando <math display="inline">k=10^{-3}</math> para todos los casos. Elegimos también las condiciones de contorno (todas de Dirichlet) e iniciales para que la solución sea la que hemos elegido.

En las Figuras [[#img-4|4]], [[#img-5|5]] y [[#img-6|6]] presentamos los resultados de los experimentos numéricos de las pruebas de convergencia en malla con solución analítica conocida, tanto para el método ASGS como para el OSS. En cada gráfica se puede apreciar con línea continua la pendiente teórica de convergencia y con línea con apéndices sobre ella la pendiente calculada. Presentamos los resultados para los distintos grados polinómicos de las funciones de forma estudiados, es decir, elementos lineales, cuadráticos, cúbicos y de cuarto orden triangulares <math display="inline">P_{1},\;P_{2},\;P_{3},\;P_{4}</math> y elementos cuadrangulares <math display="inline">Q_{1},\;Q_{2},\;Q_{3},\;Q_{4}</math>, respectivamente.

El dominio computacional es el cuadrado <math display="inline">\left[0,1\right]\times \left[0,1\right]</math> y el intervalo de tiempo es <math display="inline">\left[0,1\right]</math>. El tamaño del paso de tiempo lo hemos tomado como <math display="inline">\delta t=0.2</math>. Como se verá a continuación, en este caso el error está dominado por la aproximación espacial, no se ve afectado por la aproximación de la discretización temporal. La malla de elementos finitos consiste en triángulos o cuadrados formando una malla regular. Para cada grado polinómico de las funciones de forma <math display="inline">P_{1},\;P_{2},\;P_{3},\;P_{4}</math> o <math display="inline">Q_{1},\;Q_{2},\;Q_{3},\;Q_{4},</math> hemos calculado los errores en norma <math display="inline">L^{2}</math> para ocho mallas de tamaño <math display="inline">h</math>, siendo <math display="inline">\frac{1}{h}=15,\;20,\;25,\;30,\;35,\;40,\;45,\;50,</math> el número de partes en que se ha dividido cada lado del cuadrado. En la Tabla [[#table-2|2]] se muestra el número de elementos y el número de nodos para cada refinamiento, tanto para elementos triangulares como para elementos cuadrangulares. El integrador temporal utilizado es BDF1.

<div class="center" style="font-size: 75%;">'''Tabla 2'''. Refinamiento para elementos triangulares y cuadrangulares</div>

{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;font-size:85%;"
|- style="border-top: 2px solid;"
| colspan='2' | Triángulos
| colspan='4' | Número de nodos
|- style="border-top: 2px solid;"
| <math display="inline">1/h</math> 
| elem. 
| <math>P1</math>
| <math>P2</math>
| <math>P3</math>
| <math>P4</math>
|- style="border-top: 2px solid;"
|  15 
| 450 
| 256 
| 961 
| 2116 
| 3721
|-
| 20 
| 800 
| 441 
| 1681 
| 3721 
| 6561
|-
| 25 
| 1250 
| 676 
| 2601 
| 5776 
| 10201
|-
| 30 
| 1800 
| 961 
| 3721 
| 8281 
| 14641
|-
| 35 
| 2450 
| 1296 
| 5041 
| 11236 
| 19881
|-
| 40 
| 3200 
| 1681 
| 6561 
| 14641 
| 25921
|-
| 45 
| 4050 
| 2116 
| 8281 
| 18496 
| 32761
|- style="border-bottom: 2px solid;"
| 50 
| 5000 
| 2601 
| 10201 
| 22801 
| 40401
|}

{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;font-size:85%;"
|- style="border-top: 2px solid;"
| colspan='2' | Cuadrados
| colspan='4' | Número de nodos
|- style="border-top: 2px solid;"
| <math display="inline">1/h</math> 
| elem. 
| <math>Q1</math>
| <math>Q2</math>
| <math>Q3</math>
| <math>Q4</math>
|- style="border-top: 2px solid;"
|  15 
| 225 
| 256 
| 961 
| 2116 
| 3721
|-
| 20 
| 400 
| 441 
| 1681 
| 3721 
| 6561
|-
| 25 
| 625 
| 676 
| 2601 
| 5776 
| 10201
|-
| 30 
| 900 
| 961 
| 3721 
| 8281 
| 14641
|-
| 35 
| 1250 
| 1296 
| 5041 
| 11236 
| 19881
|-
| 40 
| 1600 
| 1681 
| 6561 
| 14641 
| 25921
|-
| 45 
| 2025 
| 2116 
| 8281 
| 18496 
| 32761
|- style="border-bottom: 2px solid;"
| 50 
| 2500 
| 2601 
| 10201 
| 22801 
| 40401

|}


En el encabezado de las gráficas de las pruebas de convergencia en malla se adjuntan dos filas con los valores de las pendientes calculadas; la primera fila corresponde a los valores de las pendientes de las rectas que pasan por los primeros 5 puntos y la segunda fila son las pendientes de las rectas que pasan por los últimos 5 puntos de los 8 puntos del refinamiento de malla evaluados. De esta manera podemos visualizar numéricamente la tendencia de la convergencia de las pendientes calculadas hacia el valor teórico. Para el cálculo de dichas pendientes se ha utilizado el método de los mínimos cuadrados. La simbología <math display="inline">m(P1),</math> <math display="inline">m(P2),</math> <math display="inline">m(P3),</math> <math display="inline">m(P4),</math> corresponde a las pendientes para elementos triangulares lineales, cuadráticos, cúbicos y de cuarto orden, respectivamente, mientras que para elementos cuadrangulares lineales, cuadráticos, cúbicos y de cuarto orden la simbología es <math display="inline">m(Q1),</math> <math display="inline">m(Q2),</math> <math display="inline">m(Q3),</math> <math display="inline">m(Q4)</math>. Los valores de las pendientes teóricas se corresponden con <math display="inline">m(P1)=m(Q1)=2</math>, <math display="inline">m(P2)=m(Q2)=3</math>, <math display="inline">m(P3)=m(Q3)=4</math>, <math display="inline">m(P4)=m(Q4)=5</math>.

Los experimentos numéricos que presentamos son los siguientes: en la Figura [[#img-4|4]] convección dominante <math display="inline">\left|\boldsymbol{a}\right|=1</math> y una pequeña reacción <math display="inline">s=10^{-3}</math>, en la Figura [[#img-5|5]] convección y reacción del mismo orden, con <math display="inline">\left|\boldsymbol{a}\right|=1</math> y <math display="inline">s=1</math>, y en la Figura [[#img-6|6]] una pequeña convección <math display="inline">\left|\boldsymbol{a}\right|=10^{-3}</math> y reacción dominante <math display="inline">s=1</math>. Para todos los ejemplos tomamos <math display="inline">k=10^{-3}</math>, como se ha indicado anteriormente.


<div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;"
|-
|[[Image:Review_640781029593-Ej1_triangulos_ASGS.png|600px|]]
|[[Image:Review_640781029593-Ej1_triangulos_OSS.png|600px|]]
|-
|[[Image:Review_640781029593-Ej1_cuadrados_ASGS.png|600px|]]
|[[Image:Review_640781029593-Ej1_cuadrados_OSS.png|600px|Convergencia ASGS-OSS, elementos triangulares y cuadrangulares. k=10⁻³, a=\left[\cos60<sup>∘</sup>,sin60<sup>∘</sup>\right], s=10⁻³.]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" style="padding-bottom:10px;"| '''Figura 4'''. Convergencia ASGS-OSS, elementos triangulares y cuadrangulares. <math>k=10^{-3},</math> <math>\mathbf{a}=\left[\cos{60}^{\circ },\mathrm{sin60^{\circ }}\right]</math>, <math>s=10^{-3}</math>
|}

<div id='img-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;"
|-
|[[Image:Review_640781029593-Ej2_cuadrados_OSS.png|600px|Convergencia ASGS-OSS, elementos triangulares y cuadrangulares. k=10⁻³, a=\left[\cos60<sup>∘</sup>,sin60<sup>∘</sup>\right], s=1.]]
|[[Image:Review_640781029593-Ej2_triangulos_ASGS.png|600px|]]
|-
|[[Image:Review_640781029593-Ej2_triangulos_OSS.png|600px|]]
|[[Image:Review_640781029593-Ej2_cuadrados_ASGS.png|600px|]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" style="padding-bottom:10px;"| '''Figura 5'''. Convergencia ASGS-OSS, elementos triangulares y cuadrangulares. <math>k=10^{-3}</math>, <math>\mathbf{a}=\left[\cos{60}^{\circ },\mathrm{sin60^{\circ }}\right]</math>, <math>s=1</math>
|}

<div id='img-6'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Review_640781029593-Ej3_triangulos_ASGS.png|600px|]]
|[[Image:Review_640781029593-Ej3_triangulos_OSS.png|600px|]]
|-
|[[Image:Review_640781029593-Ej3_cuadrados_ASGS.png|600px|]]
|[[Image:Review_640781029593-Ej3_cuadrados_OSS.png|600px|Convergencia ASGS-OSS, elementos triangulares y cuadrangulares. k=10⁻³, a=\left[\cos60<sup>∘</sup>,sin60<sup>∘</sup>\right]×10⁻³, s=1.]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" style="padding-bottom:10px;"| '''Figura 6'''. Convergencia ASGS-OSS, elementos triangulares y cuadrangulares. <math>k=10^{-3}</math>, <math>\mathbf{a}=\left[\cos{60}^{\circ },\mathrm{sin60^{\circ }}\right]\times 10^{-3}</math>, <math>s=1</math>
|}

Para los dos métodos de estabilización ASGS y OSS, los valores de las constantes algorítmicas <math display="inline">c_{i},</math> <math display="inline">i=1,2,3</math> se han calibrado a <math display="inline">c_{1}=12,</math> <math display="inline">c_{2}=2,</math> <math display="inline">c_{3}=1</math>. También es necesario indicar que para el método OSS con elementos triangulares de cuarto orden <math display="inline">P_4</math>, hemos usado en el mallado el elemento <math display="inline">P_4</math> modificado mencionado en la Sección 3.

En las gráficas de las figuras de convergencia en malla se observa tanto gráfica como numéricamente que las pendientes calculadas son mayores o tienden al valor de la pendiente teórica <math display="inline">p+1</math>. Los resultados  muestran claramente que las formulaciones ASGS y OSS son capaces de aproximar correctamente el problema de CDR escalar transitorio con todas las combinaciones de convección y reacción dominantes.

===4.2 Pruebas con capas límite===

A continuación presentamos un grupo de pruebas cualitativas para observar la estabilidad de los métodos de elementos finitos ASGS y OSS en la ecuación escalar transitoria de CDR. Se analizan cuatro casos en el dominio <math display="inline">\Omega =[0,1]\times [0,1]</math>, <math display="inline">\gamma =\left[0,1\right],</math> combinando los valores de término fuente <math display="inline">f=1</math>, <math display="inline">f=0</math> y con el coeficiente de reacción <math display="inline">s=0,</math> <math display="inline">s=1</math>, y además variando la dirección de la velocidad, considerando un valor de <math display="inline">\left|\mathbf{a}\right|=1</math> para todos los casos. El valor del coeficiente de difusión es <math display="inline">k=10^{-5}</math> para todos los ejemplos, con el fin de tener convección dominante y poder probar la estabilidad de los métodos. El tamaño del paso de tiempo considerado constante es <math display="inline">\delta t=0.2</math> y el tamaño del elemento para todas las figuras es <math display="inline">1/h=20</math>. En todos los casos excepto el caso de la Figura [[#img-8|8]], se imponen condiciones de Dirichlet homogéneas, observando inestabilidades locales en el contorno, como era de esperar, pero siendo las soluciones globalmente estables. Para el caso de la Figura [[#img-8|8]], las condiciones de Dirichlet en el contorno son: <math display="inline">u(x,0,t)=1\;\mathrm{en\;0\leq x\leq{0.3}}</math>, <math display="inline">u(x,0,t)=0\;\mathrm{en} 0.3<x\leq{1}</math>, <math display="inline">u(0,y,t)=1,\;u(1,y,t)=0,\;u(1,x,t)=0</math>. La condición inicial para todos los casos es <math display="inline">u\left(\mathbf{x},0\right)=0</math> en todo <math display="inline">\Omega </math>.

En cada figura se presentan ocho gráficos que corresponden a los cuatro grados polinómicos de las funciones de forma analizadas para cada método estabilizado de elementos finitos ASGS y OSS. La solución se presenta en <math display="inline">t= 1</math>, tiempo en el cual se ha alcanzado el estado estacionario usando el integrador temporal BDF1.

En la Figura [[#img-7|7]] presentamos el caso con término fuente <math display="inline">f=1</math>,  velocidad de convección en la dirección del eje <math display="inline">y</math> <math display="inline">\mathbf{a}=[0,1]</math>, coeficiente de reacción <math display="inline">s=0</math> y condiciones de Dirichlet homogéneas. En la Figura [[#img-8|8]] se encuentra el caso con término fuente <math display="inline">f=0</math>, velocidad de convección inclinada <math display="inline">\mathbf{a}=\left[\sin 60^{\circ },\cos 60^{\circ }\right],</math> coeficiente de reacción <math display="inline">s=0,</math> y las condiciones de Dirichlet en el contorno dadas anteriormente. En la Figura [[#img-9|9]] encontramos el caso con término fuente <math display="inline">f=1,</math> velocidad de convección inclinada <math display="inline">\mathbf{a}=\left[\sin 60^{\circ },\cos 60^{\circ }\right],</math> coeficiente de reacción <math display="inline">s=0</math> y condiciones de Dirichlet homogéneas. Y, finalmente, en la Figura [[#img-10|10]] tenemos el caso con término fuente <math display="inline">f=1,</math> velocidad de convección inclinada <math display="inline">\mathbf{a}=\left[\sin 30^{\circ },\cos 30^{\circ }\right],</math> coeficiente de reacción <math display="inline">s=1</math> y condiciones de Dirichlet homogéneas.

<div id='img-7'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|
{|  style="text-align: center; margin: 1em auto;min-width:50%;width:100%;"
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |   
| style="border-left: 2px solid;border-right: 2px solid;" | <math>ASGS</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>OSS</math>
|- style="border-top: 2px solid;"
| style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>Q1</math>
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej1_capas_ASGS_q1.png|300px|Capas_Ej1/Ej1_capas_ASGS_q1]]
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej1_capas_OSS_q1.png|300px|Capas_Ej1/Ej1_capas_OSS_q1]]
|- style="border-top: 2px solid;"
| style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>Q2</math>
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej1_capas_ASGS_q2.png|300px|Capas_Ej1/Ej1_capas_ASGS_q2]]
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej1_capas_OSS_q2.png|300px|Capas_Ej1/Ej1_capas_OSS_q2]]
|- style="border-top: 2px solid;"
| style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>Q3</math>
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej1_capas_ASGS_q3.png|300px|Capas_Ej1/Ej1_capas_ASGS_q3]]
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej1_capas_OSS_q3.png|300px|Capas_Ej1/Ej1_capas_OSS_q3]]
|- style="border-top: 2px solid;border-bottom: 2px solid;"
| style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>Q4</math>
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej1_capas_ASGS_q4.png|300px|Capas_Ej1/Ej1_capas_ASGS_q4]]
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej1_capas_OSS_q4.png|300px|Capas_Ej1/Ej1_capas_OSS_q4]]
|}
|-
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figura 7'''. Capas límite con ASGS y OSS, elementos cuadrangulares. <math>k=10^{-5}</math>, <math>\mathbf{a}=\left[0,1\right]</math>, <math>s=0</math>, <math>f=1</math> y condiciones de Dirichlet homogéneas en todo el contorno
|}


<div id='img-8'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|
{|  style="text-align: center; margin: 1em auto;min-width:50%;width:100%;"
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |   
| style="border-left: 2px solid;border-right: 2px solid;" | <math>ASGS</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>OSS</math>
|- style="border-top: 2px solid;"
| style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>Q1</math>
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej2_capas_ASGS_q1.png|300px|Capas_Ej2/Ej2_capas_ASGS_q1]]
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej2_capas_OSS_q1.png|300px|Capas_Ej2/Ej2_capas_OSS_q1]]
|- style="border-top: 2px solid;"
| style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>Q2</math>
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej2_capas_ASGS_q2.png|300px|Capas_Ej2/Ej2_capas_ASGS_q2]]
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej2_capas_OSS_q2.png|300px|Capas_Ej2/Ej2_capas_OSS_q2]]
|- style="border-top: 2px solid;"
| style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>Q3</math>
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej2_capas_ASGS_q3.png|300px|Capas_Ej2/Ej2_capas_ASGS_q3]]
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej2_capas_OSS_q3.png|300px|Capas_Ej2/Ej2_capas_OSS_q3]]
|- style="border-top: 2px solid;border-bottom: 2px solid;"
| style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>Q4</math>
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej2_capas_ASGS_q4.png|300px|Capas_Ej2/Ej2_capas_ASGS_q4]]
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej2_capas_OSS_q4.png|300px|Capas_Ej2/Ej2_capas_OSS_q4]]
|}
|-
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figura 8'''. Capas límite con ASGS y OSS, elementos cuadrangulares. <math>k=10^{-5}</math>, <math>\mathbf{a}=\left[\sin{60}^{\circ },\cos 60^{\circ }\right]</math>, <math>s=0</math>, <math>f=0</math> y condiciones de Dirichlet homogéneas en todo el contorno 
|}


<div id='img-9'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|
{|  style="text-align: center; margin: 1em auto;min-width:50%;width:100%;"
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |   
| style="border-left: 2px solid;border-right: 2px solid;" | <math>ASGS</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>OSS</math>
|- style="border-top: 2px solid;"
| style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>Q1</math>
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej3_capas_ASGS_q1.png|300px|Capas_Ej3/Ej3_capas_ASGS_q1]]
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej3_capas_OSS_q1.png|300px|Capas_Ej3/Ej3_capas_OSS_q1]]
|- style="border-top: 2px solid;"
| style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>Q2</math>
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej3_capas_ASGS_q2.png|300px|Capas_Ej3/Ej3_capas_ASGS_q2]]
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej3_capas_OSS_q2.png|300px|Capas_Ej3/Ej3_capas_OSS_q2]]
|- style="border-top: 2px solid;"
| style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>Q3</math>
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej3_capas_ASGS_q3.png|300px|Capas_Ej3/Ej3_capas_ASGS_q3]]
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej3_capas_OSS_q3.png|300px|Capas_Ej3/Ej3_capas_OSS_q3]]
|- style="border-top: 2px solid;border-bottom: 2px solid;"
| style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>Q4</math>
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej3_capas_ASGS_q4.png|300px|Capas_Ej3/Ej3_capas_ASGS_q4]]
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej3_capas_OSS_q4.png|300px|Capas_Ej3/Ej3_capas_OSS_q4]]
|}
|-
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figura 9'''. Capas límite con ASGS y OSS, elementos cuadrangulares. <math>k=10^{-5}</math>, <math>\mathbf{a}=\left[\sin 60^{\circ },\cos{60}^{\circ }\right]</math>, <math>s=0,\quad f=1</math> y condiciones de Dirichlet homogéneas en todo el contorno
|}


<div id='img-10'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|
{|  style="text-align: center; margin: 1em auto;min-width:50%;width:100%;"
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |   
| style="border-left: 2px solid;border-right: 2px solid;" | <math>ASGS</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>OSS</math>
|- style="border-top: 2px solid;"
| style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>Q1</math>
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej4_capas_ASGS_q1.png|300px|Capas_Ej4/Ej4_capas_ASGS_q1]]
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej4_capas_OSS_q1.png|300px|Capas_Ej4/Ej4_capas_OSS_q1]]
|- style="border-top: 2px solid;"
| style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>Q2</math>
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej4_capas_ASGS_q2.png|300px|Capas_Ej4/Ej4_capas_ASGS_q2]]
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej4_capas_OSS_q2.png|300px|Capas_Ej4/Ej4_capas_OSS_q2]]
|- style="border-top: 2px solid;"
| style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>Q3</math>
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej4_capas_ASGS_q3.png|300px|Capas_Ej4/Ej4_capas_ASGS_q3]]
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej4_capas_OSS_q3.png|300px|Capas_Ej4/Ej4_capas_OSS_q3]]
|- style="border-top: 2px solid;border-bottom: 2px solid;"
| style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>Q4</math>
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej4_capas_ASGS_q4.png|300px|Capas_Ej4/Ej4_capas_ASGS_q4]]
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Ej4_capas_OSS_q4.png|300px|Capas_Ej4/Ej4_capas_OSS_q4]]
|}
|-
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figura 10'''. Capas límite con ASGS y OSS, elementos cuadrangulares. <math>k=10^{-5}</math>, <math>\mathbf{a}=\left[\sin 30^{\circ },\cos 30^{\circ }\right]</math>, <math>s=1</math>, <math>f=1</math> y condiciones de Dirichlet homogéneas en todo el contorno
|}


En todos las figuras se observa que el comportamiento es el esperado en concordancia con los datos de cada figura. Además observamos que en todas la figuras existen inestabilidades en los bordes correspondientes a las capas límite, las cuales es importante notar que se atenúan al incrementar el refinamiento polinómico.

En la Figura [[#img-11|11]] encontramos una comparación de los métodos ASGS y OSS, realizando cortes sobre las gráficas de capas límite para los dos métodos ASGS y OSS del ejemplo de la Figura [[#img-7|7]], y graficando los dos cortes sobre un mismo plano cartesiano. Los cortes que encontramos en la Figura [[#img-11|11]] son cortes longitudinales y cortes trasversales. Como era de esperar, aparecen oscilaciones localizadas en las capas límites que no se propagan al interior del dominio. Para eliminarlas habría que introducir algún método de captura de discontinuidades. Se observa también que los picos espúreos del método OSS son mayores que los del método ASGS. Este comportamiento es conocido, como también lo es la mayor exactitud del método OSS, y los resultados presentados aquí corroboran este hecho para elementos de alto orden.

<div id='img-11'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|
{|  style="text-align: center; margin: 1em auto;min-width:50%;width:100%;"
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |   
| style="border-left: 2px solid;border-right: 2px solid;" | Cortes longitudinales paralelos al eje <math display="inline">y</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | Cortes transversales paralelos al eje <math display="inline">x</math>
|- style="border-top: 2px solid;"
| style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>Q1</math>
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Long_ASGS_OSS_Q1.png|300px|Capas_Ej1/Cortes_longitudinales/Long_ASGS_OSS_Q1]]
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-ASGS_OSS_Q5.png|300px|Capas_Ej1/Cortes_transversales/ASGS_OSS_Q5]]
|- style="border-top: 2px solid;"
| style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>Q2</math>
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Long_ASGS_OSS_Q2.png|300px|Capas_Ej1/Cortes_longitudinales/Long_ASGS_OSS_Q2]]
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-ASGS_OSS_Q2.png|300px|Capas_Ej1/Cortes_transversales/ASGS_OSS_Q2]]
|- style="border-top: 2px solid;"
| style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>Q3</math>
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Long_ASGS_OSS_Q3.png|300px|Capas_Ej1/Cortes_longitudinales/Long_ASGS_OSS_Q3]]
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-ASGS_OSS_Q3.png|300px|Capas_Ej1/Cortes_transversales/ASGS_OSS_Q3]]
|- style="border-top: 2px solid;border-bottom: 2px solid;"
| style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>Q4</math>
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-Long_ASGS_OSS_Q4.png|300px|Capas_Ej1/Cortes_longitudinales/Long_ASGS_OSS_Q4]]
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-ASGS_OSS_Q4.png|300px|Capas_Ej1/Cortes_transversales/ASGS_OSS_Q4]]
|}
|-
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figura 11'''. Comparación ASGS y OSS, elementos cuadrangulares. Cortes longitudinales y transversales del ejemplo de la Figura [[#img-7|7]]
|}

===4.3 Transporte de una concentración===

Presentamos a continuación los resultados numéricos para un ejemplo menos académico, consistente en el transporte y difusión de una concentración en parte de la superficie de un ducto de directriz <math display="inline">x</math>. Suponemos que tenemos un flujo constante con velocidad unidireccional <math display="inline">\mathbf{a}=\left[100,0,0\right]</math> en un ducto, como se ilustra en la Figura [[#img-12|12]]. También asumimos que sobre una región de la pared superior del ducto tenemos una concentración constante <math display="inline">u=1</math> de una variable <math display="inline">u</math>, y sobre el resto de la periferia del ducto, incluyendo la sección transversal a la entrada del flujo, asumimos el valor <math display="inline">u=0</math>. Queremos estudiar como se propaga esta concentración dentro del flujo considerado.

<div id='img-12'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 60%;;"
|-
<!-- [[Image:Review_640781029593_2608_test-Ducto en 3D.png|600px|Ducto en 3D]]--> 
|[[File:Villota_Codina_2018a_1462_Fig12.png|400px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;" | '''Figura 12'''. Ducto en 3D
|}


Considerando el ancho mucho mayor que la altura y debido a que el flujo es constante en toda la sección transversal del ducto, podemos simplificar el planteamiento del problema a dos dimensiones como vemos en la Figura [[#img-13|13]].

<div id='img-13'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 60%;"
|-
|[[File:Villota_Codina_2018a_2163_Fig13.png|400px]]
<!-- [[Image:Review_640781029593-Ducto en 2D.png|600px|Ducto en 2D]]--> 
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;" | '''Figura 13'''. Ducto en 2D
|}


El modelo matemático que resuelve el problema planteado es:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\partial _{t}u-k\partial _{i}\partial _{i}u+a_{i}\partial _{i}u  =f,\quad  \hbox{ en}~\Omega \times \gamma </math>
|-
| style="text-align: center;" | <math> u  =u_{D},\quad  \hbox{ en}~ \partial \Omega _D\times \mathbf{\gamma }, </math>
|-
| style="text-align: center;" | <math> \frac{\partial{u}}{\partial x}  = 0,\quad  \hbox{ en}~ \partial \Omega _N\times \mathbf{\gamma },  </math>
|-
| style="text-align: center;" | <math> u  =0,\quad  \hbox{ en}~\Omega \times \{ 0\} . </math>
|}
|}

Para nuestro problema, <math display="inline">k=10^{-3},</math> <math display="inline">\mathbf{a}=\left[1,0\right],</math> <math display="inline">f=0,</math> el dominio <math display="inline">\Omega </math> es el rectángulo de dimensiones <math display="inline">l_{x}=1</math> y <math display="inline">l_{y}=0,1</math> y el intervalo de tiempo es <math display="inline">\gamma =\left[0,1\right]</math>.

Las condiciones de Dirichlet en el contorno <math display="inline">\partial \Omega _D</math> donde se imponen son:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>u_D(x,y,t) =  \begin{cases}0 & x\in [0,1], ~ y = 0, ~ t \in [0,1]\\ 0 & x = 0, ~ y \in [0, 0.1], ~ t \in [0,1] \\ 0 & x \in [0, 0.05) \cup (0.5, 1],  ~ y  = 0.1,  ~ t \in [0,1]\\ 1 & x \in [0.05 , 0.5],  ~ y  = 0.1,  ~ t \in [0,1] \end{cases} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (55)
|}

En el extremo derecho (<math display="inline">\partial \Omega _N</math>) no se impone ninguna condición de Dirichlet, sino que se considera libre.

El intervalo de tiempo <math display="inline">[0,1]</math> se discretiza mediante una partición uniforme de tamaño <math display="inline">\delta t=0.2</math>. Como en los ejemplos anteriores, el integrador temporal utilizado es BDF1. La malla de elementos finitos utilizada consiste en 850 elementos cuadrangulares formando una malla regular, con 918 nodos para el elemento <math display="inline">Q_1,</math> 3535 nodos para el <math display="inline">Q_2,</math> 7852 nodos para el <math display="inline">Q_3</math> y 13869 nodos para el <math display="inline">Q_4.</math>

En la Figura [[#img-14|14]] se muestran los resultados de la simulación numérica de este problema en <math display="inline">t = 1</math> para los dos métodos ASGS y OSS y para los cuatro grados polinómicos de las funciones de forma, <math display="inline">Q_1,</math> <math display="inline">Q_2,</math> <math display="inline">Q_3,</math> y <math display="inline">Q_4,</math> observando en todos los casos que los resultados son los esperados en concordancia con los datos del problema. Se visualiza el efecto de la convección dominante, ya que se observa que la concentración <math display="inline">u</math> que se prescribe como condición de contorno se transporta al interior de dominio con poca difusión. En ninguno de los casos considerados se aprecian oscilaciones, que es el objetivo de los métodos de estabilización.

<div id='img-14'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|
{|  style="text-align: center; margin: 1em auto;width:80%;"
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |   
| style="border-left: 2px solid;border-right: 2px solid;" | <math>ASGS</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>OSS</math>
|- style="border-top: 2px solid;"
| style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>Q1</math>
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-ASGS_Q1.png|350px|Ejemplo_practico/ASGS_Q1]]
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-OSS_Q1.png|350px|Ejemplo_practico/OSS_Q1]]
|- style="border-top: 2px solid;"
| style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>Q2</math>
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-ASGS_Q2.png|350px|Ejemplo_practico/ASGS_Q2]]
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-OSS_Q2.png|350px|Ejemplo_practico/OSS_Q2]]
|- style="border-top: 2px solid;"
| style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>Q3</math>
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-ASGS_Q3.png|350px|Ejemplo_practico/ASGS_Q3]]
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-OSS_Q3.png|350px|Ejemplo_practico/OSS_Q3]]
|- style="border-top: 2px solid;border-bottom: 2px solid;"
| style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>Q4</math>
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-ASGS_Q4.png|350px|Ejemplo_practico/ASGS_Q4]]
| style="border-left: 2px solid;border-right: 2px solid;" | [[Image:Review_640781029593-OSS_Q4.png|350px|Ejemplo_practico/OSS_Q4]]
|}
|-
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figura 14'''. Transporte de una concentración
|}

==5. Conclusiones==

En este artículo hemos presentado una revisión de los métodos estabilizados ASGS y OSS y su estudio con elementos finitos de alto orden, cuadráticos, cúbicos y de cuarto orden, para resolver la ecuación escalar de CDR en problemas con convección y reacción dominantes. En el cálculo del parámetro de estabilización <math display="inline">\tau </math> hemos incluido  el grado del polinomio, para de esta manera poder utilizar elementos finitos de alto orden y mantener el orden de aproximación óptimo tanto en <math display="inline">h</math> como en <math display="inline">p</math>. Para poder usar de forma eficiente el método OSS con el elemento simplicial de cuarto orden, hemos presentado el procedimiento para construir un elemento <math display="inline">P_4</math> modificado, el cual permite aproximar la matriz de masa por una matriz diagonal con todos sus elementos distintos de cero; hemos usado este elemento modificado en todos los ejemplos con elementos de cuarto orden con el método OSS.

Los experimentos numéricos de las pruebas de convergencia en malla para elementos lineales y de alto orden revelan que la convergencia de los dos métodos estabilizados es óptima, para todos los casos de convección y reacción dominantes, tanto para elementos triangulares como para elementos cuadrangulares, y para todos los grados polinómicos de las funciones de forma que hemos considerado. Hemos observado que las pendientes de convergencia del método OSS son ligeramente mayores que las del método ASGS, lo cual indica mayor precisión del método OSS. Además, también se concluye que para un número de nodos dado, los errores son menores cuando mayor es el grado del polinomio de forma utilizado, lo cual indica que para obtener un menor coste computacional para una precisión dada podría ser conveniente utilizar elementos finitos de alto orden. Sin embargo, no ha sido nuestro objetivo estudiar con detalle el coste computacional de los métodos presentados.

En términos generales, los resultados de las pruebas de capas límite han sido los esperados y están en concordancia con las condiciones impuestas en cada caso. Asimismo, hemos observado que el método OSS produce más oscilaciones localizadas en las capas límites, como ya era conocido para el caso de elementos lineales y cuadráticos.

Finalmente, los resultados del ejemplo de la simulación numérica del transporte de una concentración han sido satisfactorios, de acuerdo con el fenómeno físico planteado. En este caso, las diferencias de los resultados entre los métodos ASGS y OSS son insignificantes.

==Agradecimientos==

R. Codina agradece la ayuda recibida a través del programa ICREA Academia, del Gobierno de Catalunya.

==Referencias==

<div class="auto" style="width: auto; margin-left: auto; margin-right: auto;font-size: 85%;">

<div id="cite-1"></div>
[[#citeF-1|[1]]]  Brooks A.N.,  Hughes T.J.R. Streamline upwind/Petrov-Galerkin formulations for      convection dominated flows with particular emphasis on the     incompressible Navier-Stokes equation.  Computer Methods in Applied Mechanics and Engineering, 32:199-259, 1982.

<div id="cite-2"></div>
[[#citeF-2|[2]]]  Kelly D.W., Nakazawa S.,   Zienkiewicz O.C.,   Heinrich J.C.  A note on upwinding and anisotropic balancing dissipation in finite element approximations to convective diffusion problems. International Journal for Numerical Methods in Engineering,  15:1705-1711, 1980.

<div id="cite-3"></div>
[[#citeF-3|[3]]] Johnson C.,  Nävert U.,  Pitkäranta J. Finite element methods for linear hyperbolic equations. Computer Methods in Applied Mechanics and Engineering, 45:285-312, 1984.

<div id="cite-4"></div>
[[#citeF-4|[4]]] Hughes T.J.R.,  Franca L.P.,  Hulbert G.M. A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective-diffusive equations. Computer Methods in Applied Mechanics and Engineering, 73:173-189, 1989.

<div id="cite-5"></div>
[[#citeF-5|[5]]] Franca L.P.,  Frey S.L.,   Hughes T.J.R.  Stabilized finite element methods: I. Application to the advective-diffusive model.  Computer Methods in Applied Mechanics and Engineering, 95:253-276, 1992.

<div id="cite-6"></div>
[[#citeF-6|[6]]] Brezzi F.,   Franca L.P.,    Hughes T.J.R.,  Russo A.  <math>b=\int g</math>. Computer Methods in Applied Mechanics and Engineering, 145:329-339, 1997.

<div id="cite-7"></div>
[[#citeF-7|[7]]] Douglas J., Russel T. Numerical methods for convection dominated problems based on combining the method of characteristics with finite elements or finite difference procedures. SIAM Journal on Numerical Analysis, 9:871-885, 1982.

<div id="cite-8"></div>
[[#citeF-8|[8]]] Donea J.  A Taylor-Galerkin method for convection transport problems. International Journal for Numerical Methods in Engineering, 20:101-119, 1984.
<div id="cite-9"></div>
[[#citeF-9|[9]]] Codina R.  Comparison of some finite element methods for solving the diffusion-convection-reaction equation.  Computer Methods in Applied Mechanics and Engineering, 156:185-210, 1998.

<div id="cite-10"></div>
[[#citeF-10|[10]]]  Huges T.J.R.  Multiscale phenomena: Green's functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods. Computer Methods in Applied Mechanics and Engineering, 127:387-401, 1995.

<div id="cite-11"></div>
[[#citeF-11|[11]]] Huges T.J.R.,  Feijóo G.O.,  Mazzei L.,  Quincy J.B.  The variational multiscale method-a paradign for computational mechanics. Computer Methods in Applied Mechanics and Engineering, 166:3-24, 1998.

<div id="cite-12"></div>
[[#citeF-12|[12]]] Codina R., Badia S.,  Baiges J., Principe J. Variational multiscale methods in computational fluid dynamics.  In Encyclopedia of Computational Mechanics, Part 1. Fluids, John Wiley & Sons Ltd, 2017.

<div id="cite-13"></div>
[[#citeF-13|[13]]] Codina R.  Finite element approximation of the convection-diffusion equation: Subgrid-scale spaces, local instabilities and anisotropic space-time discretizations. In Bail 2010 - Boundary and Interior Layers, Computational and Asymptotic Methods, Part of the Lecture Notes in Computational Science and Engineering, C. Clavero, J.L. Gracia and F.J. Lisbona (Eds.), 81:85-97, Springer, 2011.

<div id="cite-14"></div>
[[#citeF-14|[14]]] Codina R.  Stabilization of incompressibility and convection through orthogonal sub-scales in finite element methods. Computer Methods in Applied Mechanics and Engineering, 190:1579-1599, 2000.

<div id="cite-15"></div>
[[#citeF-15|[15]]] Codina R. On stabilized finite element methods for linear systems of convection-diffusion-reaction equations. Computer Methods in Applied Mechanics and Engineering, 188:61-82, 2000.

<div id="cite-16"></div>
[[#citeF-16|[16]]] Codina R. A stabilized finite element method for generalized stationary incompressible flows. Computer Methods in Applied Mechanics and Engineering, 190:2681-2706, 2001.

<div id="cite-17"></div>
[[#citeF-17|[17]]] Codina R. Stabilized finite element approximation of transient incompressible flows using orthogonal subscales. Computer Methods in Applied Mechanics and Engineering, 191:4295-4321, 2002.

<div id="cite-18"></div>
[[#citeF-18|[18]]] Codina R.,  Principe J., Guasch O., Badia S.  Time dependent subscales in the stabilized finite element approximation of incompressible flow problems. Computer Methods in Applied Mechanics and Engineering, 196:2413-2430, 2007.

<div id="cite-19"></div>
[[#citeF-19|[19]]] Baiocchi C., Brezzi F.,  Franca L.P. Virtual bubbles and Galerkin/least-squares type methods (Ga.L.S). Computer Methods in Applied Mechanics and Engineering, 105:125-141, 1993.

<div id="cite-20"></div>
[[#citeF-20|[20]]] Codina R., Blasco J.  Analysis of a stabilized finite element approximation of the transient convection-diffusion-reaction equation using orthogonal subscales. Computing and Visualization in Science, 4:167-174, 2002.

<div id="cite-21"></div>
[[#citeF-21|[21]]] Principe J., Codina R.  On the stabilization parameter in the subgrid scale approximation of scalar convection-diffusion-reaction equations on distorted meshes. Computer Methods in Applied Mechanics and Engineering, 199:1386-1402, 2010.

<div id="cite-22"></div>
[[#citeF-22|[22]]] Codina R.  On <math display="inline">hp</math> convergence of stabilized finite element methods for the convection-diffusion equation. SeMA Journal, 75:591-606, 2018.

<div id="cite-23"></div>
[[#citeF-23|[23]]] Houston P., Süli E.  Stabilised <math display="inline">hp</math>-finite element approximation of partial differential equations with nonnegative characteristic form. Computing, 66:99-119, 2001.

<div id="cite-24"></div>
[[#citeF-24|[24]]] Lube G., Rapin G.  Residual-based stabilized higher-order FEM for advection-dominated problems. Computer Methods in Applied Mechanics and Engineering, 195:4124-4138, 2006.

<div id="cite-25"></div>
[[#citeF-25|[25]]] Schwab Ch.  <math>p</math>- and <math display="inline">hp</math>-Finite element methods. Theory and application to solid and fluid mechanics. Oxford University Press, 388 pp., 1998.

</div>