INTRODUCCIÓN El objetivo del Análisis Numérico es proporcionar métodos convenientes para obtener soluciones útiles (resultados numéricos) de ciertos problemas matemáticos. En ciertos casos, puede ocurrir que, aunque el problema tenga solución y ésta sea única, no exista un método analítico que nos la permita calcular (por ejemplo, en el cálculo de ciertas integrales). En otros casos, sí existen métodos analíticos que nos conducen a la solución exacta, pero éstos pueden requerir un gran número de operaciones o no ser isibles para su implementación en un ordenador; por ejemplo, según se ha visto en la resolución de sistemas lineales, la regla de Cramer conducía a la solución exacta, pero en la práctica, esto no será cierto para sistemas con una cierta dimensión, por el elevado número de operaciones que requiere. En los casos anteriormente planteados, se buscan soluciones que, sin coincidir plenamente con la solución exacta del problema se aproximen a ella tanto como sea posible. En definitiva, el Análisis Numérico es una parte de las Matemáticas que estudia los métodos constructivos para la resolución de problemas matemáticos de una forma aproximada. El conjunto de reglas que definen el método numérico recibe el nombre de algoritmo. El fin de todo algoritmo es llegar a un programa, que implementado en un ordenador, proporcione la solución del problema planteado. Interesará, obviamente, elegir algoritmos que produzcan resultados lo más precisos posibles.
MATRIZ ORTOGONAL Una matriz ortogonal es una matriz cuadrada cuya matriz inversa coincide con su matriz traspuesta. El conjunto de matrices ortogonales constituyen una representación lineal del grupo ortogonal. Geométricamente las matrices ortogonales representan transformaciones isométricas en espacios vectoriales reales1 (o más exactamente espacios de Hilbert reales) llamadas justamente, transformaciones ortogonales. Estas transformaciones son isomorfismos internos del espacio vectorial en cuestión. En el caso real dichas transformaciones pueden ser rotaciones, reflexiones especulares o inversiones y son usadas extensivamente en computación gráfica. Por sus propiedades, también son usadas para el estudio de ciertos fibrados y en física se las usa en el estudio del movimiento de cuerpos rígidos y en la formulación de ciertas teorías de campos.
TEOREMA A n × n: A es ortogonal ssi AT · A = I. Observe que el teorema anterior se deduce de que para dos vectores x y y en Rn , x • y = x ′ · y:
Con lo anterior se deduce que cuando se hace AT·v se calcula un vector donde cada componente es el producto punto de la columna correspondiente de A con el vector v. Con lo anterior se deduce que cuando se calcula AT · A la matriz resultante tiene en la posición (i, j) justo ai • aj es decir, el producto punto de la columna i de A con la columna j de A. De esta forma: AT · A = I si y solo si se tiene que las columnas de A son ortogonales y que tienen norma 1. Con la observación anterior presente podemos hacer el ejemplo 1 más fácilmente. Ejemplo 1 Indique si el conjunto formado por los siguientes vectores es ortogonal
Y calculamos AT · A:
que sean cero los elementos que están fuera de la diagonal principal indica que el conjunto es ortogonal.
Ejemplo 2 Determina los valores de x, y y z para que el conjunto de vectores
sea ortogonal. Formamos la matriz A cuyas columnas son los vectores:
Y calculamos AT · A:
Por tanto, para que el conjunto se ortogonal debe cumplirse que:
de donde, los ´unicos valores que hacen ortogonal al conjunto son x = −31/5, y = 1/15 y z = −14/5 Ejemplo 3 Determine el vector de coordenadas de v =< 2, 2, −4 > respecto a la base ortonormal
Recordemos que el vector de coordenadas de un vector respecto a una base son los coeficientes de la combinación lineal de la base que da tal vector. Si la base es orto normal entonces los coeficientes de la combinación lineal son los coeficientes de Fourier, es decir los productos punto del vector con cada uno de los elementos de la base. Verifiquemos primero que el conjunto es orto normal. Para ello, formamos la matriz A cuyas columnas son los vectores de B:
y calculamos AT · A:
Por tanto, dando la matriz diagonal el conjunto es ortogonal; dando la identidad el conjunto es orto normal. Para calcular los productos punto del elemento de B con v recurrimos al producto:
Por tanto, c1 = v • u1 = −2/3, c2 = v • u2 = −4/3, y c3 = v • u3 = 14/3 y el vector de coordenadas de v respecto a la base B es < −2/3, −4/3, 14/3 >. Teorema Sea A una matriz n × n, y u y v dos vectores en Rn. Entonces
DEMOSTRACIÓN Teorema Sea A una matriz n × n. Son equivalentes las siguientes afirmaciones: (1) A es ortogonal. (2) A preserva los productos punto:
(3) A preserva norma:
Demostración (1) implica (2) Si A es ortogonal, AT A = I. Así (A u) • (A v) = (A u) T · A v = u TAT · A v = u T · (AT · A)v = u T · I · v = u T · v = u (2) implica (3)
Se tiene ||A v||2 = (A v) • (A v) = v • v = || tomando raíz cuadrada se tiene la igualdad de (3). (3) implica
PROPIEDADES
De la definición, es inmediato que si una matriz es ortogonal, la matriz es no singular o invertible y su transpuesta coincide con su inversa El determinante de una matriz ortogonal A es +1 ó -1. En efecto, de las propiedades del determinante tenemos det (A. At) = det A det At = det A det A = (det A) ² = det I = 1, por tanto, det A = +-1 El conjunto de matrices n x n ortogonales, junto con la operación de producto de matrices es un grupo llamado grupo ortogonal O(n). Supongamos que A y B son matrices ortogonales y sea C igual al producto de A por B. Usando las propiedades del producto de matrices, tenemos
y así, el producto de matrices ortogonales es una matriz ortogonal.
En teoría de grupos, al grupo de matrices ortogonales n por n con coeficientes en el cuerpo K se denomina grupo ortogonal de dimensión n y se representa con O (n, K). En particular el subgrupo formado por las matrices ortogonales de determinante +1, se llama grupo especial ortogonal y se le representa con SO(n,K). Entre las matrices ortogonales se encuentran las matrices de rotación y las de permutación. Cuando el cuerpo es el de los reales K=R entonces se escribe simplemente O(n) := O(n,K) y SO(n) SO(n,K).
SENSIBILIDAD EN SISTEMAS LINEALES Actualmente los modelos matemáticos que se están utilizando para investigar fenómenos físicos son cada vez más realistas. Las nuevas características de los modelos utilizan a menudo parámetros cuyos valores no pueden ser conocidos exactamente. A raíz de esto hay una necesidad de realizar análisis de sensibilidad paramétrica de los modelos representados por sistemas algebraicos diferenciales. Las áreas de uso incluyen optimización, valoración del parámetro, simplificación del modelo, control óptimo, sensibilidad de proceso, análisis de la incertidumbre, y diseño experimental para una amplia gama de los problemas científicos y de ingeniería.
1.CÁLCULO DE LA SENSIBILIDAD PARA SISTEMAS DE ECUACIONES ALGEBRAICOS DIFERENCIALES POR EL MÉTODO ADJUNTO La resolución de problemas modelados con sistemas DAE es un tema de importancia al cual se dedican esfuerzos considerables desde hace 30 años (Brenan 1996 y Kees 1999). Muchos problemas de ingeniería y ciencia se modelan en forma natural con sistemas DAE, donde en general las ecuaciones son del tipo:
Compuesta por una combinación de ecuaciones diferenciales y algebraicas, donde, Nx, x, x, ∈ Rᵑ son las variables de estado, y Np p∈ Rᵑ es un vector de parámetros. El problema del cálculo de la sensibilidad de parámetros para un sistema DAE se plantea de la siguiente forma: dado un sistema DAE del tipo (1) dependiente de parámetros Np p∈ Rᵑ, encontrar dx/dp en el tiempo T para j:1,..,Np . La solución de este problema requiere la solución simultánea del sistema DAE original con los Np sistemas de sensibilidad, obtenidos por diferenciación del DAE original con respecto a cada parámetro (Cao 2003). Para sistemas grandes esto puede parecer mucho trabajo, pero puede ser realizado eficientemente si Np es relativamente pequeño, explotando el hecho de que los sistemas de la sensibilidad son lineales, y todos comparten las mismas matrices Jacobianas del sistema original. No obstante, algunos problemas requieren el estudio de las sensibilidades con respecto a una gran cantidad de parámetros. Para estos problemas, particularmente si el número de las variables del estado Nx es también grande, la estimación previa de la sensibilidad es intratable. En Cao (2003) se demostró que estos problemas se pueden manejar más eficientemente por el método adjunto (Errico 1997). En este sentido el interés es calcular la sensibilidad dG/dp de una función objetivo G (x,p ) definida como:
O en su defecto la sensibilidad dg/dp de una función g (x,T, p) definida solamente en el tiempo T . Por ejemplo, en el modelado matemático de reactores nucleares, un caso común de análisis es el
estudio de sensibilidad de la velocidad del líquido refrigerante Ui, perturbando la potencia térmica demandada Q. Para ello, se puede plantear la siguiente función objetivo: donde u es la media de la velocidad del líquido. A continuación de describe la aplicación del método adjunto en el cálculo de la sensibilidad.
1.1 Sensibilidad para G (x, p) Para hallar una expresión de dG/dp se define una función objetivo aumentado I (x, p)
Donde λ es un multiplicador de Lagrange, * denota la matriz traspuesta, y F (x, x, p, t) = 0 según (1). Derivando ambos con respecto a p, se obtiene que la sensibilidad respecto de G (x, p) es:
Luego se desarrolla el término λ (t)* Fx Xp para obtener una expresión con respecto a Xp:
Reemplazando la ecuación (5) en (4), y reagrupando los términos se obtiene que:
Donde:
Es la ecuación adjunta, y λ es la incógnita o variable adjunta del sistema (7). En vez de trabajar con el sistema (7), en este trabajo se va a utilizar la notación de sistema adjunto aumentado:
Luego, la sensibilidad para G (x,p) queda representada por:
A partir de (9) se observa que para calcular la sensibilidad por el sistema adjunto, el sistema (8) debe ser resuelto hacia atrás hasta el tiempo t = 0. Las condiciones iniciales en el tiempo t T = para calcular dicha sensibilidad surgen de la resolución del sistema:
Para calcular el término integral cuadratura β y la ecuación, adjunto (8):
de la ecuación (9) se utiliza una variable de la cual se incorpora a la resolución del sistema
Obteniendo de esta forma para t = 0:
Con lo cual la expresión de la sensibilidad resulta:
1.2 sensibilidad para g (x, t, p) Para el cálculo de la sensibilidad dg/ dp, se aplica dg /dp= d/dT dG /dp = a la ecuación (9), con lo cual se obtiene que:
Donde la variable adjunta λ (t, T) depende de los parámetros t y T, y donde λ t (t, T) es dλ (t, T) /dT. A partir de (14) se obtiene que el sistema adjunto es:
Y se resuelve de manera similar que el sistema (8). Las condiciones iniciales en el tiempo t= T se obtienen de resolver el sistema:
Nuevamente para calcular la integral cuadratura β y la ecuación
de la ecuación (14), se utiliza una variable de , la cual se incorpora a la resolución del sistema adjunto (15):
Por último la expresión de la ecuación de sensibilidad está dada por:
2) RESULTADOS NUMÉRICOS Para la solución numérica de los sistemas (8) y (15) (problemas de condiciones iniciales) actualmente se dispone de una variedad de métodos numéricos que pueden clasificarse en tres categorías: métodos de un paso, multipaso y de extrapolación. La eficiencia de cada tipo de aproximación depende del problema particular. En este trabajo se resuelven los sistemas (8) y (15) transformándolos en sistemas puramente algebraicos mediante aproximaciones numéricas (Boroni 2005). Para ello se define una transformación reemplazando λ y λ & por funciones que representan estimadores multipaso de las variables adjuntas y sus derivadas, es decir:
Específicamente en este trabajo se emplea un estimador implícito para las variables adjuntas y un estimador multipaso lineal clásico (Jackson 1980) en las derivadas, es decir:
Para analizar la aplicación del método adjunto en el cálculo de la sensibilidad se realizaron una serie de experimentos, donde se examinaron los resultados y errores numéricos. Para dicho estudio se definieron dos sistemas de prueba: un sistema de ecuaciones diferenciales lineales con solución explícita analítica y un sistema oscilatorio conservativo no lineal. Para simplificar la presentación de los resultados se considera únicamente el cálculo de la sensibilidad de la función objetivo ( , ) Gxp .
2.1 Sistema lineal Se plantea el ejemplo de una ecuación diferencial lineal de primer orden dado por:
El cual posee un único parámetro p, y donde la solución explícita es:
El objetivo de trabajar con este sistema lineal es obtener medidas del error asociados a los términos en el cálculo numérico de la sensibilidad dG dp definida en (9). Para el cálculo de la sensibilidad dG/ dp se define la siguiente función objetivo:
Utilizando (9) y (23) se obtiene que la solución de la ecuación de sensibilidad es:
La evolución temporal de la sensibilidad dG /dp puede verse en la figura 1, la cual tiene forma de campana sesgada a izquierda con un valor máximo
En la figura 2 se observa el error absoluto entre el cálculo numérico de y la solución explícita. En la misma se puede ver que el cálculo del término integral es el que genera mayor error, con un máximo aproximado a 1.0 x 10-4, el cual se estabiliza para un tiempo t ≥ 3. En cuanto al error vinculado al término λ (t)t=0, el valor máximo es inferior a 1.5x10−5 y se da cuando t ≈ 0.28; luego de alcanzar este valor el error disminuye de manera asintótica a 0. Por último, los resultados del error son bajos comparados con el paso de tiempo dt =1.0 10−3 utilizado para el cálculo numérico de la sensibilidad. En la figura 3 se muestra el error absoluto entre la solución explícita y numérica de la sensibilidad dG/ dp.