1. No category

2. FUNDAMENTOS TEÓRICOS
21
Capítulo 2
FUNDAMENTOS TEÓRICOS DEL
CÁLCULO VECTORIAL DISCRETO
2.1 LA RETÍCULA TIPO.
En este capítulo se va a presentar el cálculo vectorial discreto que, siguiendo la filosofía
de los métodos miméticos, nos proporcionará una herramienta de cálculo para la resolución numérica de edp’s. Esto exige que se establezcan antes ciertas definiciones y proposiciones que nos ayuden a cimentar la estructura que vamos a edificar, lo que nos
permitirá dotar de coherencia y completitud al discurso matemático que se va a presentar.
El primer paso consiste en establecer el “soporte” de nuestra herramienta de cálculo.
Este soporte hace referencia al medio en que nos vamos a encontrar. Por supuesto será
un medio discreto, esto es, un conjunto finito o infinito numerable de puntos en el espacio n-dimensional. ¿Pero es suficiente con esta definición?.
Si el objetivo es establecer un análogo al caso continuo vamos a tener que buscar alguna
“relación” que represente las propiedades topológicas que se dan en él. Así pues se llega
de forma natural a la definición de retícula.
Definición 2.1: Retícula
Sea V un conjunto discreto de puntos contenido en n .
Sea χ una aplicación de conexión entre los elementos de V .
χ se puede interpretar como la aplicación:
χ : V × V 
→
→ χ ( x, y )
( x, y ) 
que cumple las siguientes propiedades:
i)
χ ( x , x) = 0
no hay autolazos (no es una aplicación reflexiva)
ii)
χ ( x, y) = 1
si x e y están conectados o son adyacentes
χ ( x , y) = 0
en otro caso
iii)
χ ( x , y ) = χ ( y , x ) propiedad de simetría
Entonces el par (V , χ ) recibe el nombre de retícula y se denota por Γ = (V , χ ) .
2. FUNDAMENTOS TEÓRICOS
22
En estas condiciones al conjunto V se le llama conjunto de nodos de la retícula y a χ
se le llama relación de adyacencia de la retícula.
Pues bien, imaginemos que ahora representamos geométricamente la conexión entre los
nodos mediante segmentos rectilíneos o ramas. Sea E el conjunto formado por dichas
ramas. Entonces llamamos gráfico de la retícula o simplemente grafo al par (V , E ) y lo
denotamos por Gf = (V , E ) . Este grafo simplemente tiene la misión de informar visualmente de la conectividad de la retícula, sin tener mayor importancia para el cálculo
que aquí se desarrolla. En la figura 2.1 puede observarse el grafo de una retícula irregular en un espacio tridimensional.
V
E
Figura 2.1 Retícula irregular en un espacio tridimensional.
Fuente: elaboración propia.
De la retícula también nos interesará conocer su tamaño, pues nos indicará la envergadura del sistema lineal asociado que deberemos resolver. Así, llamaremos tamaño de la
retícula al número de elementos contenidos en V y lo denotaremos por V . En general,
de ahora en adelante, la notación A denotará el cardinal de un conjunto finito Α dado.
Después de estas definiciones básicas el objetivo ahora es desarrollar la teoría con vistas
a establecer con facilidad un cálculo vectorial operativo. Es interesante saber, por ejemplo, qué nodos están conectados a uno dado. Si bien ya lo podemos saber a través de la
relación de adyacencia χ , será más útil dar una expresión explícita. Así definimos el
conjunto N ( x ) = { xi ∈ V tal que χ ( x , xi ) = 1} como el conjunto de los nodos conectados o adyacentes a x. Y siguiendo la notación anterior, llamaremos grado del nodo x a
N ( x ) . Esto es muy importante porque como se verá más adelante el grado del nodo
condiciona la dimensión del tensor métrico y de los campos vectoriales asociados.
2. FUNDAMENTOS TEÓRICOS
23
De hecho, podemos avanzar un poco más generalizando la definición de conjunto de
nodos adyacentes. Así, el conjunto de nodos adyacentes a un conjunto de nodos Α dado
vendrá determinado por la expresión:
N (A ) =
U N ( y ) = U { x ∈V tal que χ ( y , x ) = 1}
i
y∈A
(2.1)
i
y∈A
En estas condiciones también nos podremos plantear la composición de conjuntos adyacentes, es decir, el conjunto determinado por N ( N ( x ) ) , que representa los nodos adyacentes a los adyacentes de x. Diremos entonces que estos nodos están en relación de
adyacencia 2 respecto de x. Y por tanto, una notación lógica y natural será denotar
N ( N ( x ) ) por N 2 ( x ) . Extrapolando esta notación, se obtiene que el conjunto de nodos en orden de adyacencia j se expresará como N
nos será muy útil al hablar de las plantillas de cálculo.
j
( x) = N
( N (K N ( x ) ) ) , lo cual
j veces
Así pues, y como resumen, tenemos que una retícula genérica, irregular o no, viene determinada por un conjunto arbitrario V de puntos contenido en n , con una relación de
adyacencia arbitraria χ , que define un grado arbitrario para cada nodo.
Está claro que una retícula de este estilo será un tipo muy general y nos interesará más
trabajar con retículas que guarden ciertas propiedades útiles para simplificar el esfuerzo
de cálculo. Aparece, pues, la idea de retícula uniforme.
Definición 2.2: Retícula uniforme no acotada.
Sea W = {w 1 , w 2 ,K , w n } una base de
n
. W quedará caracterizada por la matriz de
{e1 , e2 ,K , en } . Sea M la matriz de cambio.
cambio respecto la base cartesiana
↑
↑ 


M =  w1 L w n 
↓
↓ 

Entonces, diremos que la retícula Γ h = (Vh , χ ) es homogénea no acotada de paso h>0 si
cumple:
i)
ii)
Vh = { x = h ⋅ λM : λ ∈
∀x, y ∈ Vh , ∀λ ∈
n
}
conjunto infinito numerable
n
χ ( x , y ) = χ ( x + h ⋅ λM, y + h ⋅ λM )
Que podríamos denominar como regla del paralelogramo,
tal y como puede observarse en la figura adjunta.
2. FUNDAMENTOS TEÓRICOS
λM
y
x
λM
24
y+hλM
x+hλM
Figura 2.2. Regla del paralelogramo en las retículas uniformes.
Fuente: elaboración propia.
En estas condiciones diremos que W es la base de la retícula, donde las direcciones de
los vectores base {w i } , i = 1,K , n definen las líneas de red. Diremos también que Vh es
el conjunto homogéneo no acotado de nodos o puntos de red, y que χ es la relación
homogénea de adyacencia (o de conexión).
Asimismo diremos que la retícula uniforme es ortogonal (respectivamente ortonormal)
si W es una base ortogonal (respectivamente ortonormal).
Sería interesante dar algunos ejemplos de las retículas uniformes que nos podemos encontrar. Por lo general, las más útiles van a ser las retículas rectangulares y las triangulares, por ser las más simples y las que rellenan mejor el espacio. En las figuras adjuntas
pueden observarse los grafos (locales) de estas retículas tanto en 2D como en 3D. Es
importante remarcar que la teoría que se desarrolla en esta tesina es válida para espacios
de dimensión n, lo que implica que las formulaciones para dos y tres dimensiones no
son más que casos particulares.
En primer lugar podemos ver las retículas rectangulares en 2D. En la figura 2.3 se
muestra el caso más sencillo, que corresponde a tener una base de la retícula ortonormal
y tener los pasos iguales tanto en la dirección x como en la dirección y. De esta forma la
retícula se compone de cuadrados.
h
h
w2
w1
Figura 2.3. Retícula rectangular 2D con una base ortonormal.
Fuente: elaboración propia.
2. FUNDAMENTOS TEÓRICOS
25
El caso más general de retícula rectangular en 2D se muestra en la figura 2.4. Consiste
en tener una base de retícula cuyos vectores formen un ángulo θ y donde los pasos sean
distintos según cada línea de red. En este caso se componen paralelogramos.
w2
θ
w1
h2
h1
Figura 2.4. Caso general de la retícula rectangular 2D.
Fuente: elaboración propia.
A las retículas triangulares en 2D se les va a dedicar especial atención en este estudio.
Precisamente el anexo 1 se dedica específicamente a hacer un análisis aplicado a este
tipo de retículas como ejemplo de todo el desarrollo teórico de la tesina. Esta particular
elección se debe a que las retículas triangulares representan un caso más general que las
rectangulares, y contienen las complicaciones de cálculo típicas que nos podremos encontrar. De hecho, se puede decir que, en general, cuanto más conexa sea una retícula
más dificultad habrá en el cálculo. Pues bien, en la figura 2.5 se puede observar una
retícula triangular genérica en la que la base no es ortogonal y los pasos son distintos
según cada línea de red.
h2
h1
w2
θ
w1
Figura 2.5. Retícula triangular en 2D. Fuente: elaboración propia.
2. FUNDAMENTOS TEÓRICOS
26
Veamos ahora algunos ejemplos en 3D. Al aumentar de dimensión van a ocurrir dos
cosas. La primera es que vamos a tener muchos más nodos que antes. La segunda es que
vamos a poder establecer relaciones de conexión mucho más complejas entre los nodos.
Todo esto nos va a llevar a una mayor dificultad para obtener las formulaciones de cálculo. Sin embargo, va a ser interesante comprobar la potencia del método trabajando
con este tipo de retículas.
En la figura 2.6 se observa el grafo de la retícula rectangular en 3D teniendo en cuenta
una base ortonormal, es decir, una retícula compuesta por cubos. Si en el caso 2D teníamos 13 nodos en un ámbito local, ahora podemos contar 25. Posteriormente se matizará qué se entiende por ámbito local cuando se hable de las plantillas de cálculo. Ahora
baste con saber que son aquellos nodos más “cercanos” a uno dado según la relación de
adyacencia o conexión.
w3
w1
w2
Figura 2.6. Retícula rectangular en 3D con una base ortonormal.
Fuente: elaboración propia.
Como último ejemplo se muestra el caso de una retícula triangular en 3D o, mejor dicho, una retícula de tetraedros. Ahora las relaciones de conexión son muy numerosas y
difíciles de dibujar. El entramado de ramas muestra visualmente la dificultad que implica trabajar operativamente en estas condiciones.
2. FUNDAMENTOS TEÓRICOS
27
w3
w1
θ
w2
Figura 2.7. Retícula de tetraedros en 3D con una base no ortogonal.
Fuente: elaboración propia.
De las retículas uniformes podemos obtener dos propiedades importantes que se detallan a continuación.
Proposición 2.1: Simetría central local.
Sea Γ h = (Vh , χ ) una retícula uniforme.
Sean los nodos y = x + h ⋅ λM , ŷ = x − h ⋅ λM , para cada x ∈ Vh , λ ∈
El nodo ŷ representa el opuesto de y respecto el centro x .
y
x
ŷ
Figura 2.8. Un nodo y su opuesto.
Fuente: elaboración propia.
Entonces se cumple:
y ∈ N ( x ) ←
→ yˆ ∈ N ( x )
n
2. FUNDAMENTOS TEÓRICOS
28
Es decir, si un nodo “y” está conectado a “x”, también lo estará su opuesto. Y viceversa,
si un nodo “y” no está conectado a “x”, entonces su opuesto tampoco lo estará.
Demostración:
En virtud de que es una retícula
uniforme
y ∈ N ( x ) ←
→ χ ( x , y ) = 1←
→ χ ( x + h ⋅ λ *M , y + h ⋅ λ *M ) = 1 , ∀λ * ∈
n
Ya que la elección de λ * es arbitraria, vamos a escoger λ * = −λ . Con lo que ahora se
obtiene lo siguiente:
1 = χ ( x − h ⋅ λM , y − h ⋅ λM ) = χ ( yˆ , x ) = χ ( x , yˆ ) ←
→ yˆ ∈ N ( x ) ■
Por la propiedad de simetría
Proposición 2.2: El grado de cualquier nodo de una retícula uniforme no acotada es el
mismo y es un número par. Matemáticamente lo expresamos como:
Grado ( x ) = N ( x ) = 2m , m ∈
∀x ∈ Γ h
Demostración:
De la definición de relación de adyacencia homogénea se desprende que el grado de
cualquier nodo ha de ser el mismo. Por otro lado, el que sea un número par, viene de la
propiedad anterior, puesto que si en una determinada dirección hay un nodo conectado,
también lo estará su opuesto, y por tanto habrá dos nodos conectados por cada una de
las m direcciones de la retícula. ■
Por el momento, todo lo visto hasta ahora hace referencia a retículas no acotadas, es
decir, no hemos impuesto restricciones en cuanto a su extensión. Pero en la práctica nos
vamos a encontrar siempre con retículas que sí que van a estar acotadas, puesto que vamos a resolver problemas de contorno de la física e ingeniería. Esto obliga a profundizar
más en la definición de las retículas, distinguiendo claramente qué entendemos por nodos interiores y nodos frontera.
Si uno observa cómo se actúa en el continuo, parece claro que debemos encontrar una
analogía para la bola de radio epsilon entorno a un punto, tan útil para definir nodos
interiores, nodos frontera, puntos de adherencia, de acumulación, etc. Este concepto lo
vamos a encontrar en la plantilla de cálculo o también llamada molécula computacional,
que se define a continuación.
Definición 2.3: Plantilla de cálculo de orden α (o molécula computacional de orden
α ).
Sea Γ h = (Vh , χ ) una retícula homogénea no acotada.
2. FUNDAMENTOS TEÓRICOS
29
Llamaremos plantilla de cálculo de orden α , α ∈ , al conjunto de nodos definido como:
S α ( x ) = { x} U N ( x ) U N 2 ( x )UK U N α ( x )
(2.2)
o equivalentemente:
S
α
α
( x ) = { x} U N j ( x )
(2.3)
j =1
Esto es lo que entendemos como ámbito local, término al que se ha hecho alusión en los
párrafos precedentes. La idea consiste en caracterizar al nodo x y a todos sus “vecinos”,
entendiendo por vecinos aquellos que están conectados según χ .
El tamaño de la plantilla de orden α se define como S α ( x ) y es interesante conocerlo
para determinar la cantidad de nodos que intervienen al ensamblar la matriz del sistema
lineal asociado, el cual nos proporcionará la solución aproximada a nuestra edp.
Una propiedad importante que se tratará más adelante es que el tamaño de la plantilla de
orden α se corresponde con el número de nodos del esquema en diferencias para el
operador diferencial de orden α . Si bien esto se enunciará para los operadores diferenciales de orden 1 y 2, es fácil inferir el resultado para cualquier α ∈ .
Pues bien, ahora ya estamos en condiciones de presentar las retículas uniformes acotadas con sus nodos interiores y nodos frontera.
Definición 2.4: Retícula uniforme acotada.
La retícula Γ h = (Vh , χ ) de paso h>0 se dice que es uniforme acotada si se cumple:
i) Vh = { x = h ⋅ λM : λ ∈ Η ⊂
n
} donde H es un conjunto finito.
ii) χ es la relación de adyacencia homogénea definida para Γ h = (Vh , χ ) .
Está claro que si Γ h = (Vh , χ ) es una retícula uniforme no acotada y Γ h = (Vh , χ ) es
una retícula uniforme acotada, entonces se cumple que Γ h ⊂ Γ h , es decir, la primera es
una parte de la segunda.
Definición 2.5: Los nodos interiores y los nodos frontera.
Sea Γ h = (Vh , χ ) una retícula uniforme acotada. Entonces diremos que
→ {∃α ∈
x es nodo interior de Vh ←
tal que S α ( x ) ⊂ Vh }
o
Al conjunto de nodos interiores de Vh se le denotará por Vh , mientras que al conjunto
o
formado por los nodos frontera se le denotará por ∂Vh y se definirá como ∂Vh = Vh / Vh .
2. FUNDAMENTOS TEÓRICOS
30
Llegados a este punto puede observarse cómo la analogía es evidente entre la plantilla
de cálculo y la bola del espacio euclideano continuo. Así, el orden α hace las veces de
radio epsilon, teniendo en cuenta para un caso la función de conexión y para el otro la
función de distancia. El concepto resulta más claro al observar las figuras adjuntas. La
primera se corresponde con la red rectangular 3D y la segunda con la red de tetraedros,
ambas para un orden de adyacencia 1. Se incluyen sólo los casos de tres dimensiones
por ser los más explicativos gráficamente.
ε
ε
Figura 2.9. Plantilla para la retícula rectangular 3D. Figura 2.10. Plantilla para la retícula de tetraedros.
Fuente: elaboración propia.
Después de lo visto, quisiera hacer una observación al respecto. Con esta definición
dada para la retícula uniforme acotada resulta que no hay adaptación perfecta a una
geometría curva de la frontera. Esto se debe a que en todo momento se exige la regularidad, y el contorno no tiene porqué serlo. De todas formas, la mejor o peor adaptación
al contorno vendrá determinada por el paso h que se escoja para la retícula y por la relación de adyacencia de los nodos. Así, las retículas triangulares por lo general se adaptarán mejor a los contornos irregulares. No obstante, otra posibilidad hubiera sido hacer
los pasos variables. Estaríamos entonces ante una retícula localmente irregular, o como
también suelo llamar, una retícula casi-uniforme. Pero éste ya no es el ámbito de estudio
de esta tesina.
Nodos
frontera
Plantilla de
orden 1
Nodos
exteriores
Nodos
interiores
Figura 2.10. Adaptación a la frontera. Fuente: elaboración propia.
∂Ω
2. FUNDAMENTOS TEÓRICOS
31
De hecho, a efectos de esta tesina, el desarrollo teórico del cálculo vectorial discreto se
va a realizar en los nodos interiores de la retícula, evitando deliberadamente el análisis
en el contorno. Esto es porque la herramienta de cálculo vectorial posibilita un desarrollo teórico extenso en la frontera, tal como por ejemplo los teoremas integrales de la
divergencia y de Stokes. Por ello, para los objetivos de esta tesina, resultaría excesivo
abordarlo.
Además sólo voy a trabajar con retículas uniformes, por ser las que ayudan a comprender mejor el desarrollo teórico, pero también porque hay dos motivos por los que
resulta más conveniente trabajar con ellas. Así, en primer lugar, el trabajar con retículas
irregulares nos haría variar las expresiones de los operadores nodo a nodo, con las dificultades que implicaría de cara a su implementación informática. En segundo lugar,
para mayor eficiencia computacional de cálculo, los ordenadores requieren retículas con
una estructura simple para poder realizar fácilmente las operaciones. Por tanto, las retículas uniformes constituyen el ámbito de estudio más idóneo para su aplicación computacional.
Por último, es importante observar que, lo que puede condicionar o no la regularidad de
nuestra retícula, es la geometría de la frontera. En cambio, en el interior del dominio, no
hay ningún impedimento para usar retículas uniformes, y es precisamente allí donde se
va realizar el estudio. Así que está más que justificado el estudio de las retículas uniformes en esta tesina.
2.2. LA NOTACIÓN DE RETÍCULA.
El propósito de este apartado es establecer una notación práctica que nos permita definir
fácilmente las expresiones de los operadores diferenciales en la retícula.
Como cualquier otra estrategia numérica, la discretización mimética del continuo va a
generar un sistema de ecuaciones algebraicas que deberemos resolver para obtener la
solución aproximada de nuestro problema elíptico. Esto lo podemos representar matricialmente con la expresión:
R ⋅u = f
(2.4)
donde R es la matriz del sistema lineal, que llamaremos matriz de rigidez
u es el vector cuyas componentes son los valores que toma la función incógnita
en los nodos de la retícula.
f es el vector cuyas componentes son los valores que toma la función dato en
los nodos de la retícula.
Lo lógico será numerar los nodos de forma que nos permita ensamblar fácilmente el
sistema lineal de ecuaciones. Por ello definiremos la función de numeración Num ( x )
como aquella aplicación que asigna a cada nodo de la retícula un número natural J que,
por cuestiones de notación, irá en mayúsculas por hacer referencia a la retícula global.
Matemáticamente se expresará como sigue:
2. FUNDAMENTOS TEÓRICOS
32
Num : V 
→
x 
→ Num ( x ) = J
Como es una aplicación biyectiva, también podremos pensar en la función inversa
Num −1 , esto es, dado un número J podremos conocer cuál es el nodo asociado
x = Num −1 ( J ) . Esto nos puede resultar útil cuando trabajemos con algoritmos de cara a
una implementación en programas informáticos.
Pero con esto aún no es suficiente para trabajar de forma operativa. Necesitamos conocer la relación de conexión χ entre los nodos. Para ello haremos uso de la notación
matricial, que caracterizará biunívocamente las retículas uniformes. Así definiremos F
como la matriz de adyacencia o conexión de la retícula, cuya estructura será la siguiente:
↑

F =  v1
↓

↑
↑

v 2 L vm 
↓
↓ 
Donde los vectores v k ∈
n
n filas
serán aquellos que cumplan:
χ ( x , x + h ⋅ v k ) = 1 ∀x ∈ Vh
es decir, x + h ⋅ v k es un nodo adyacente. Estos vectores vendrán referidos generalmente
a la base W de la retícula. Para obtener su expresión en la base cartesiana, será suficiente
con multiplicar la matriz F por la matriz de cambio M.
Hay que observar que por la propiedad de simetría central local de las retículas uniformes, sólo es necesario dar m vectores, puesto que los m restantes se obtienen de:
v m+k = − v k k = 1,K , m . A los nodos v k se les llamará vectores de adyacencia o conexión. A modo de ejemplo, en la figura 2.12 se muestran los vectores de adyacencia
para la retícula triangular 2D y en la figura 2.13 para la retícula de tetraedros en 3D.
v3
v2
v1
w2
w1
Figura 2.12. Vectores de adyacencia para la retícula triangular 2D.
Fuente: elaboración propia.
2. FUNDAMENTOS TEÓRICOS
33
Expresión de los vectores: m = 3 n = 2
v1 = w1 = (1, 0 )W
v 2 = w1 + w 2 = (1,1)W
v 3 = w 2 = ( 0,1)W
v5
v4
v6
v7
v1
v2
v3
w3
w1
w2
Figura 2.13. Vectores de adyacencia para la retícula de tetraedros.
Fuente: elaboración propia.
Expresión de los vectores: m = 7 n = 3
v1 = w1 = (1, 0, 0 )W
v 5 = w 3 − w 2 = ( 0, −1,1)W
v 2 = w1 + w 2 = (1,1, 0 )W
v 6 = w 3 − w 2 − w 1 = ( −1, −1,1)W
v 3 = w 2 = ( 0,1, 0 )W
v 7 = w 3 − w1 = ( −1, 0,1)W
v 4 = w 3 = ( 0, 0,1)W
En ambos casos m es mayor que n, es decir, hay más vectores de adyacencia que dimensión en el espacio. Esto quiere decir que se trata de retículas con gran conectividad, y
hace que aparezcan relaciones dependientes entre los vectores adyacentes. Este hecho
será de gran importancia en los capítulos 3 y 5.
Volviendo con la matriz de adyacencia F, es importante comentar que es invariante en
cada uno de los nodos interiores de la retícula si ésta es uniforme. Es por eso que nos va
a ser muy útil para obtener las expresiones de los operadores discretos.
Para terminar con esta sección dedicada a la notación de las retículas deberíamos hablar
de las plantillas de cálculo. Pero esto se hará más adelante, en la sección 3.2. El hecho
2. FUNDAMENTOS TEÓRICOS
34
de que no se incluya aquí se debe a que faltan aún por definir algunos conceptos que
ayudarán a explicar mejor el porqué de su elección.
2.3. EL ESPACIO TANGENTE.
Uno de los conceptos clave para desarrollar todo el cálculo vectorial que aquí se presenta es el de espacio tangente, pues nos va a proporcionar el rigor matemático para cimentar todo el aparato vectorial. Primero veamos cuál es su definición formal y después
trataremos de entenderlo mejor.
Definición 2.6: El espacio tangente.
Para cada x ∈ Vh , definimos el espacio tangente en x como el espacio vectorial Tx ( Γ h )
de las combinaciones lineales formales de los segmentos incidentes con x, que denominamos sxy , donde y ∈ N ( x ) . Estos segmentos representan las ramas de conexión en el
nodo x. De esta forma, el sistema {sxy }
Tx ( Γ h ) son de la forma v =
∑vs
y xy
y∈N ( x )
es una base de Tx ( Γ h ) y los elementos de
donde v y ∈
para cada y ∈ N ( x ) . Notar que
y∈N ( x )
con esta definición dim Tx ( Γ h ) = N ( x ) para cada x ∈ Vh .
Se trata, pues, de un espacio vectorial con la particularidad que está definido sobre una
variedad discreta. Para mayor información sobre las variedades discretas consultar [8].
Lo que nos va a ayudar más para entender el concepto de espacio tangente es la geometría diferencial. Así, por ejemplo, el plano tangente a una superficie representa un subespacio vectorial del espacio tangente y gráficamente es como mejor podremos representarlo. En estas condiciones, el plano tangente está generado formalmente por las derivadas parciales de las variables en las que se parametriza la superficie (ver figura).
∂
∂θ
∂
∂ϕ
Figura 2.14. Representación del plano tangente a una superficie. Fuente: elaboración propia.
2. FUNDAMENTOS TEÓRICOS
35
Dado que en nuestro caso, el medio discreto, no podemos hablar de diferenciabilidad,
deberemos escoger unos elementos que hagan un papel análogo. Estos elementos son
los segmentos sxy , ya que informan de una determinada direccionalidad y sentido. Así
que podemos definirlos como base de nuestro espacio tangente y por tanto la dimensión
del mismo vendrá determinada por el número de segmentos que tengamos, es decir, el
grado del nodo.
Una forma simplificada de entenderlo sería considerar un plano tangente en cada uno de
los nodos de la retícula, pero teniendo en cuenta que en realidad no es un plano, sino
todo un espacio de dimensión N ( x ) , que en las retículas uniformes equivale a 2m.
Por último quisiera aclarar algunos conceptos que pueden dar lugar a confusión. No se
deben confundir los elementos del espacio tangente con los vectores de n . Me explico.
Los vectores de adyacencia v k k = 1,K , m se obtienen como combinación lineal de los
vectores de la base vectorial W , es decir, pertenecen a n . Sin embargo no pertenecen
al espacio tangente, pues no se pueden obtener como combinación lineal de los segmentos sxy . Recordemos que la dimensión de los vectores v k es n y la dimensión del espacio tangente es 2m.
2.4. DESARROLLO DEL CÁLCULO VECTORIAL DISCRETO.
En esta sección se van a definir los operadores que constituirán la herramienta de cálculo para obtener la solución aproximada de las edp’s elípticas. Pero antes será necesario
introducir unos conceptos previos de carácter formal.
Siguiendo el espíritu de la tesina, se presentará la teoría como la analogía al caso continuo. Por ello, se considerarán las retículas Γ h de paso h>0 como variedades discretas.
Así, se hará uso del espacio tangente, ya visto en la sección anterior. Esto nos permitirá
definir los conceptos de campo vectorial, campo de matrices y los productos internos en
el espacio de funciones de retícula y en el espacio de los campos vectoriales.
Empezaremos denotando por C (Vh ) al conjunto de las funciones reales definidas en Vh ,
esto es, en el conjunto de nodos de la retícula. De esta forma las funciones que pertenezcan a C (Vh ) tendrán un soporte definido por los nodos cuyo valor de la función sea
no nulo, es decir, si u ∈ C (Vh ) , entonces Soporte ( u ) = { x ∈ Vh u ( x) ≠ 0} .
Está claro que las funciones pertenecientes al conjunto C (Vh ) tendrán un soporte discreto, y por tanto, no necesitaremos hablar de propiedades de diferenciabilidad.
El siguiente paso es definir qué se entiende por campo vectorial. Lo podemos ver como
una aplicación que asigna a cada nodo un vector de su espacio tangente. Así, si f es un
campo vectorial entonces para cada x ∈ Vh , f ( x ) = ∑ f ( x, y ) s xy , donde recordemos
y∈N ( x )
36
2. FUNDAMENTOS TEÓRICOS
que {s xy }
y∈N ( x )
representa la base del espacio tangente en x. De forma equivalente, un
campo vectorial puede representarse a través de su función componente, es decir, la función f : Vh × Vh 
→ tal que para cada x ∈ Vh , f ( x , y ) = 0 si y ∉ N ( x ) .
El conjunto de campos vectoriales definidos en Γ h lo denotaremos por X (Vh ) . Al igual
que antes, el soporte de un campo vectorial f perteneciente a X (Vh ) será el conjunto
finito de nodos donde el valor de f es diferente de 0.
Profundicemos un poco más en los campos vectoriales. Si f es un campo vectorial diremos que es no negativo si su función componente es no negativa, y lo denotaremos
por f ≥ 0 . Además, diremos que f representa un campo vectorial no negativo, cuya
función componente vendrá dada por f . Por extensión, si f, g ∈ X (Vh ) , la inecuación
f ≥ g representará que el campo vectorial dado por la diferencia es no negativo, es decir, f − g ≥ 0 .
En el contexto de las retículas uniformes, llamaremos flujo a un campo vectorial f cuya
función componente satisfaga que f ( x , y ) = − f ( y , x ) para todo x, y ∈ Vh . El nombre
de flujo se debe a la mecánica de medios continuos, donde el flujo saliente de una propiedad a través de una superficie es igual al flujo entrante cambiado de signo, puesto
que la normal a la superficie tiene el signo opuesto.
Continuando con más definiciones para las retículas uniformes, diremos que el campo
vectorial f es homogéneo si existe un vector b = ( b j ) ∈ 2m , donde 2m es el grado de
los nodos de la retícula, tal que f ( x , x j ) = b j para cada x ∈ Vh y cada x j ∈ N ( x ) ,
j = 1,K , 2m . En este caso diremos que el campo homogéneo f está determinado por el
vector b , y será un flujo si se cumple que b j = −bm + j , j = 1,K , m .
Bien, hemos definido el espacio de funciones de retícula y el espacio de campos vectoriales sobre la retícula. Para seguir con nuestro objetivo debemos ahora definir los productos internos que nos permitirán mediante técnicas de composición y dualidad obtener los operadores en diferencias. Es importante remarcar que estos productos internos
se definen para soportes finitos de las funciones u , v ∈ C (Vh ) y de los campos
f, g ∈ X (Vh ) .
De esta forma, si f , g son las funciones componente de los campos f, g ∈ X (Vh ) , la
expresión ( f , g ) denotará la función perteneciente a C (Vh ) definida por:
( f , g )( x ) = ∑
y∈N ( x )
f ( x , y ) ⋅ g ( x , y ) , x ∈ Vh
Así, definiremos el producto interno en el espacio X (Vh ) como:
(2.5)
37
2. FUNDAMENTOS TEÓRICOS
f ,g
X
=
1
1
( f , g )( x ) = ∑ ∑ f ( x , y ) ⋅ g ( x , y )
∑
2 x∈Vh
2 x∈Vh y∈N ( x )
(2.6)
donde el factor ½ es debido a que cada segmento adyacente se tiene en cuenta dos veces. Recordemos que la relación de adyacencia cumple la propiedad de simetría. También hay que observar que esta definición recuerda en gran medida al producto escalar
de los campos vectoriales continuos, sin más que cambiar el sumatorio por la integral.
En la misma línea definiremos el producto estándar en C (Vh ) , es decir, con el caso continuo como modelo a seguir, pero teniendo en cuenta un soporte discreto, lo cual nos
llevará a la siguiente expresión:
u,v
C
=
∑ u ( x) ⋅ v ( x)
x∈Vh
u, v ∈ C (Vh )
(2.7)
Sin embargo, con todo esto, aún no es suficiente para determinar los operadores en diferencias. Falta por introducir el concepto de campo de matrices en la retícula.
Así, al igual que los campos vectoriales, un campo de matrices sobre Γ h se define como
la aplicación A que asigna a cada nodo x ∈ Vh una matriz A ( x ) cuadrada de orden
dim Tx (Vh ) , que equivale justamente al grado del nodo x. Veamos qué tipos de campos
de matrices nos podemos encontrar y cuales van a ser los más idóneos para nuestros
propósitos.
Diremos que A es un campo de matrices no singulares si para cada x ∈ Vh la matriz
A ( x ) es no-singular. En este caso denotaremos por A −1 al campo de matrices inversas.
También diremos que el campo de matrices A es diagonal, simétrico o definido positivo si para cada x ∈ Vh la matriz A ( x ) es diagonal, simétrica o definida positiva, respectivamente. De esta forma diremos que A es un tensor métrico si es un campo simétrico
y definido positivo. Además, un tensor métrico A es diagonal si y sólo si la base del
espacio tangente es ortogonal. Observar que esto no son más que propiedades que se
infieren del caso continuo.
En el contexto de las retículas uniformes, la matriz A ( x ) será siempre de orden 2m y
tendrá la siguiente estructura:
 A ( x ) A2 ( x ) 
A ( x) =  1

 A3 ( x ) A4 ( x ) 
donde A1 ( x ) , A 2 ( x ) , A 3 ( x ) , A 4 ( x ) son matrices cuadradas de orden m.
En estas condiciones, diremos que el campo de matrices A es isotrópico por bloques si
para cada x ∈ Vh existen A1 ( x ) , A 2 ( x ) matrices simétricas de orden m tales que:
38
2. FUNDAMENTOS TEÓRICOS
 A ( x ) A2 ( x ) 
A ( x) =  1

 A 2 ( x ) A1 ( x ) 
Además diremos que A es isotrópico si es isotrópico por bloques y existen funciones
a1 , a2 , a3 ∈ C (Vh ) tales que para cada x ∈ Vh , A1 ( x ) = a1 ( x ) I + a2 ( x ) (J − I ) y
A 2 ( x ) = a3 ( x ) I + a2 ( x ) (J − I ) , donde I es la matriz identidad y J es la matriz cuyos
coeficientes son iguales a 1.
Siguiendo la notación que se propone para las retículas, el coeficiente de la matriz
A ( x ) que corresponde a la fila j-ésima y a la columna k-ésima se expresará como
a jk ( x ) , donde los índices j,k hacen referencia a la numeración local de los nodos adya-
centes a x, es decir, si x j = x + h ⋅ v j y xk = x + h ⋅ v k , entonces:
a jk ( x ) ≡ ax ( x j , xk ) = ax ( x + h ⋅ v j , x + h ⋅ v k )
j , k = 1,K , 2m
Sobre el tema de la numeración local de la plantilla de cálculo se dará más información
en la sección 3.2. Ahora sigamos con más definiciones de los campos de matrices.
En el contexto de las retículas uniformes, al igual que con los campos vectoriales, podemos definir campos de matrices homogéneos de la siguiente forma; si existe una matriz A = ( aij ) de orden 2m tal que aij ( x ) = aij , i, j = 1,K , 2m para cada x ∈ Vh , entonces diremos que el campo de matrices es homogéneo y está determinado por la matriz
A = ( aij ) .
Por supuesto todas las propiedades hasta ahora expuestas son aplicables al campo de
matrices inverso siempre y cuando el campo sea de matrices no-singulares. En particular, si A es un tensor métrico, entonces es no-singular y A −1 también es un tensor métrico.
Ahora ya estamos en condiciones de definir el cálculo operativo sobre las retículas. Empezaremos por definir el gradiente como nuestro operador básico a partir del cual obtendremos todos los demás, tal y como es usual en el caso continuo.
Definición 2.7: El operador gradiente.
El operador gradiente sobre Γ h asigna a cada u ∈ C (Vh ) el campo vectorial ∇u dado
por la siguiente expresión:
2m
∇u ( x ) =
∑
j =1
( u ( x j ) − u ( x) ) ⋅ s
xj − x
xx j
x ∈ Vh , x j ∈ N ( x ) , s xx j ∈ Tx (Vh )
que en las retículas uniformes podemos definir de forma equivalente como
(2.8)
39
2. FUNDAMENTOS TEÓRICOS
2m
∇u ( x ) =
∑
( u ( x + h ⋅ v j ) − u ( x) ) ⋅ s
h ⋅ vj
j =1
donde
xx j
x ∈ Vh , s xx j ∈ Tx (Vh )
hace referencia a la norma euclideana en
n
(2.9)
y h es la frecuencia de paso de
la retícula.
Está claro que ∇u ∈X (Vh ) cuando u ∈ C (Vh ) y que ∇u es siempre un flujo. Además,
∇u = 0 si y sólo si u es una función constante, tal y como era de esperar. Notar, además,
que haciendo el límite cuando h tiende a 0, la definición propuesta se corresponde con la
derivada direccional del cálculo diferencial clásico.
Definamos ahora el operador divergencia como aquel que asigna a cada f ∈X (Vh ) la
función div f ∈ C (Vh ) determinada por la relación:
u , div f
C
= − f , ∇u
para cada u ∈ C (Vh ) .
X
(2.10)
Que se puede expresar de forma equivalente como:
∑ u ( x ) ⋅ div f ( x ) = − 2 ∑ (f , ∇u )( x )
1
x∈Vh
(2.11)
x∈Vh
Esta relación se obtiene de la fórmula de Green, que puede utilizarse también como caracterización de la divergencia. De esta forma puede definirse formalmente la divergencia como el opuesto del operador gradiente adjunto respecto a los productos internos
definidos en C (Vh ) y X (Vh ) . Esto es, matemáticamente, div = −∇∗ en C (Vh ) .
Así pues, si f es la función componente del campo vectorial f , entonces:
1
div f ( x ) =
2
2m
∑
f ( x , xj ) − f ( xj , x)
xj − x
j =1
, x ∈ Vh , x j ∈ N ( x )
(2.12)
o de forma equivalente para las retículas uniformes:
1
div f ( x ) =
2h
2m
∑
f ( x , xj ) − f ( xj , x)
vj
j =1
, x ∈ Vh , x j ∈ N ( x )
(2.13)
Vamos a definir también dos operadores más, denotados por A∇ y ( b ,∇ ) y dados por
las siguientes expresiones:
A∇u ( x ) =
2m



2m
∑∑
i =1
j =1
aij ⋅
u ( x j ) − u ( x) 
xj − x
 ⋅ s xxi , x ∈ Vh

(2.14)
2. FUNDAMENTOS TEÓRICOS
2m
( b , ∇u )( x ) = ∑ b j ⋅
j =1
u ( x j ) − u ( x)
xj − x
40
, x ∈ Vh
(2.15)
En notación matricial, según las definiciones dadas, tenemos que:

A∇u ( x ) = 

A
 

 ⋅ ∇u ( x )  y
 

 

( b , ∇u )( x ) = [
b



] ⋅ ∇u ( x) 


Por supuesto, A∇u ∈ X (Vh ) y ( b , ∇u ) ∈ C (Vh ) . Notar además que se toman homogéneos el campo de matrices A y el campo de vectores b, puesto que sólo vamos a trabajar
con operadores diferenciales elípticos de coeficientes constantes.
Por último, para terminar con esta sección, comentar que todos los operadores definidos
hasta ahora, a saber, ∇ , div, A∇ y ( b ,∇ ) , sólo toman los valores de la función incógnita u en x y en los nodos de orden de adyacencia 1, es decir, en la plantilla de orden 1.
Esto sugiere que una forma de indicar el orden de estos operadores discretos sea la plantilla de cálculo en la que se expresan. Así, un operador que tome valores en la plantilla
de orden 1 será de orden 1, y se deberá cumplir que un operador que tome valores en la
plantilla 2 es de orden 2. Pero esto lo veremos en el siguiente capítulo.