sábado, 29 de octubre de 2016

Los cuatros espacios fundamentales de una matrix

Cada que se considera una matriz $A\in \mathbb{F}^{n,m}$, donde $\mathbb{F}$ puede ser $\mathbb{R}$ o $\mathbb{C}$, esta inducen naturalmente dos aplicaciones lineales $A:\mathbb{F}^{m}\to \mathbb{F}^{n}$ y $A^{\top}:\mathbb{F}^{n}\to \mathbb{F}^{m}$, las que cuales permiten determinar cuatro subespacios vectoriales, dos subespacios vectoriales en $\mathbb{F}^{n}$ y dos subespacios vectoriales en $\mathbb{F}^{m}$, conocidos usualmente como espacio nulo, espacio fila, espacio nulo izquierdo y espacio columna de $A$ respectivamente. Por lo tanto el objetivo en esta ocasión es estudiar estos conjuntos y las relaciones entre ellos.

Para comenzar veamos la definición de nulo y rango de una matriz $A$. El espacio nulo de la aplicación $A:\mathbb{F}^m\to \mathbb{F}^n$ es el conjunto $N(A)\subseteq \mathbb{F}^{m}$ dado por \begin{equation} N(A)=\{x\in \mathbb{F}^{n}:Ax=0\} \end{equation} y el espacio column o rango es el conjunto $R(A)\subseteq \mathbb{F}^{n}$ dado por \begin{equation} R(A)=\{y\in \mathbb{F}^{n}:y=Ax, \forall\,x\in \mathbb{F}^{m}\}. \end{equation}

En la Figura 1 se muestra el espacio nulo y el espacio rango asociados con la matriz $A$. En el gráfico se puede apreciar que para la función $A:\mathbb{F}^{m}\to \mathbb{F}^{n}$ el espacio nulo es un subconjunto de $\mathbb{F}^{m}$, mientras que el espacio columna es un subconjunto de $\mathbb{F}^{n}$. Recuerda que la ecuación homogénea $Ax=0$ siempre tiene como solución trivial $x=0$. Esta propiedad implica que $O_{_{\mathbb{F}^{m}}}\in \mathbb{F}^{m}$ también pertenece a $N(A)$ y $O_{_{\mathbb{F}^{n}}}\in \mathbb{F}^{n}$ también pertenece a $R(A)$.

Figura 1. El nulo y el rango de la aplicación $A:\mathbb{F}^m\to \mathbb{F}^n$.

El lector puede observar que realmente $A$ tiene «cuatro subespacios vectoriales» asociados que son $N(A)$, $R(A)$, $N(A^{\top})$ y $R(A^\top)$. Observe que el espacio nulo y el espacio columna son subespacios de espacios diferentes. Más precisamente una matriz $A\in \mathbb{F}^{m,n}$ define las aplicaciones $A:\mathbb{F}^{m}\to \mathbb{F}^{n}$, $A^{\top}:\mathbb{F}^{n}\to \mathbb{F}^{m}$ y los conjuntos $N(A)$, $R(A^\top)\subseteq \mathbb{F}^{m}$ y $N(A^\top), R(A)\subseteq \mathbb{F}^{n}$.

En efecto para ver que $N(A)$, $R(A)$, $N(A^{\top})$ y $R(A^\top)$ son subespacios vectoriales, recordemos que un subconjunto $U\subset \mathbb{F}^{m}$ es cerrado bajo combinaciones lineales si para cada par de elementos $x,y\in U$ y todos los escalares $a, b\in \mathbb{F}$ se tiene que $ax+by\in U$.

Los conjunto $N(A)$ y $R(A)$ son cerrados bajo combinaciones lineales dado la función inducida por $A$ es una transformación lineal. En efecto, considere dos elementos arbitrarios $x_{_1}, x_{_2}\in N(A)$, entonces $Ax_{_1}=0$ y $Ax_{_2}=0$. Entonces, para cualquier par de elementos $a,b\in \mathbb{F}$ se tiene \begin{equation} A(ax_{_1}+bx_{_2})=aAx_{_1}+bAx_{_2}=0 \end{equation} por lo tanto $(ax_{_1}+bx_{_2})\in N(A)$. Por lo tanto, $N(A)\subseteq \mathbb{F}^m$ es cerrado bajo combinaciones lineales. Análogamente, considere dos elementos $y_{_1}, y_{_2}\in R(A)$, esto es, existe $x_{_1}, x_{_2}\in \mathbb{F}^{m}$ tal que $y_{_1}=Ax_{_1}$ y $y_{_2}=Ax_{_2}$. Entonces, para cualquier $a, b\in \mathbb{F}$ se tiene \begin{equation} (ay_{_1}+by_{_2})=aAx_{_1}+bAx_{_2}=A(ax_{_1}+bx_{_2}) \Rightarrow (ay_{_1}+by_{_2})\in R(A). \end{equation} Por lo tanto, $(A)\subset \mathbb{K}^{n}$ es cerrado bajo combinaciones lineales.

Por otro lado si $A=[A_{_{:1}},\dots, A_{_{:n}}]$, entonces cualquier elemento $y\in R(A)$ puede ser expresado como $y=Ax$ para algún $x\in \mathbb{F}^{n}$, esto es, \begin{equation} y=Ax =[A_{_{:1}},\dots, A_{_{:n}}]\left[\begin{array}{c} x_{_1}\\ \vdots \\ x_{_n} \end{array}\right]=A_{_{:1}}x_{_1}+\cdots + A_{_{:n}}x_{_n}\in gen(\left\{A_{_{:1}}x_{_1},\dots, A_{_{:n}}x_{_n}\right\}). \end{equation} Por está razón es que el espacio columna y el rango son el mismo espacio. El lector puede verificar que cualquier elemento de la forma $y=A_{_{:1}}x_{_1}+\cdots + A_{_{:n}}x_{_n}$ es también un elemento de $R(A)$, por lo tanto $y=Ax$.

El nulo y el rango de una matriz cuadrada caracterizan si la matriz es invertible o no. Esto se establece en el siguiente resultado:

Teorema 1. Dada una matriz $A\in \mathbb{F}^{n,n}$, entonces las siguientes afirmaciones son equivalentes:
  • La matriz $A^{-1}$ existe;
  • $N(A)=\{0\}$;
  • $R(A)=\mathbb{F}^{n}$.

Veamos ahora algunas de las relaciones que existen entre los cuatro espacios asociados a un par de matrices $A$ y $B$, tales que la matriz $B$ se pueda obtener mediante operaciones de fila sobre la matriz $A$. Asumiendo la siguiente notación: $A\stackrel{fila}{\longleftrightarrow}B$ para indicar que la matriz $A$ se puede transformar en la matriz $B$ mediante operaciones de Gauss sobre las filas de $A$. Entonces se tiene el siguiente teorema:

Teorema 2. Sea la matrices $A, B \in \mathbb{R}^{m,n}$, entonces:
  • $A\stackrel{fila}{\longleftrightarrow}B\Longleftrightarrow N(A)=N(B)$;
  • $A\stackrel{fila}{\longleftrightarrow}B\Longleftrightarrow R(A^\top)=R(B^\top)$.

Es fácil de ver que este último resultado es verdadero, dado que las operaciones de Gauss no cambian las soluciones de sistemas lineales, por lo que el sistema lineal $Ax = 0$ tiene exactamente las mismas soluciones de $x$ como el sistema lineal $Bx = 0$, es decir, $N (A) = N (B)$. La segunda propiedad también es cierta, ya que las operaciones de Gauss en las filas de $A$ son equivalentes a las operaciones de Gauss sobre las columnas de $A^{\top}$. Ahora bien, es fácil de ver que cada una de las operaciones de Gauss sobre las columnas de $A^{\top}$ no cambian el $R(A^{\top})$, por lo tanto, $R(A^{\top})=R(B^{\top})$.

Sin embargo no podemos pensar que el párrafo anterior es suficiente para una prueba del teorema. Veamos una presentación más detallada de las ideas anteriores. Una forma es utilizar la multiplicación de matrices para expresar la propiedad de que las operaciones de Gauss no cambian las soluciones de sistemas lineales. Uno puede probar lo siguiente: Si se tiene las matrices $A$ y $B$ de orden $n\times m$ relacionadas por las operaciones de Gauss en sus filas, entonces existe una matriz $G$ de orden $m × m$ invertible $GA = B$. La prueba de esta propiedad es simple. Cada una de las operaciones de Gauss se asocia con una matriz invertible, $E$, llamada una matriz de Gauss elemental. Cada matriz de Gauss elemental es invertible, ya que cada operación de Gauss siempre se puede revertir. El resultado de varias operaciones de Gauss sobre una matriz $A$ es el producto de las matrices de Gauss elementales apropiados en el mismo orden que se realizan las operaciones de Gauss. Si se obtiene la matriz $B$ de la matriz $A$, haciendo operaciones de Gauss dadas por matrices $E_{_i}$, para $i = 1,\dots k$, en ese orden, podemos expresar el resultado del método de Gauss de la siguiente manera: \begin{equation} E_{_k}\cdots E_{_1}A=B\hspace{0.5cm} G=E_{_k}\cdots E_{_1} \Longrightarrow GA =B. \end{equation} Donde cada matriz elemental de Gauss es una matriz invertible, y por lo tanto $G$ también es invertible.

Considere dos matrices $A$ y $B$ de orden $m\times n$ estan relacionadas por operaciones de Gauss en sus filas, entonces existe una matriz de orden $G$ de orden $m\times m$ tal que $GA = B$. Esta observación es la clave para mostrar que $N(A) = N(B)$, ya que dado cualquier elemento $x \in N(A)$ \begin{equation} Ax=0\Longleftrightarrow GAx=0, \end{equation} donde la equivalencia se sigue del hecho de que $G$ es invertible. Entonces es simple (¡Lo simple es una provocación para que tu lo verifiques!) ver que \begin{equation} 0=GAx=Bx \Longleftrightarrow x\in N(B). \end{equation} Por lo tanto, se tiene que $N(A)=N(B)$. Ahora veamos la afirmación apuesta. Si $N(A)=N(B)$ esto significa que sus formas escalonada reducidas $E_{_A}, E_{_B}$ son los mismas, es decir, $E_{_A}=E_{_B}$. Esto significa que existen operaciones de Gauss en las filas $A$ que la transforman en la matrix $B$

Ahora mostraremos que $R(A^{\top})= R(B^{\top})$. Considere un elemento $x\in R(A^{\top})$, entonces existe un elemento $y\in \mathbb{F}^{m}$ tal que \begin{equation} x=A^{\top}y = A^{\top}G^{\top}(G^{\top})^{-1}y=(GA)^{\top}\bar{y}=B^{\top}\bar{y},\hspace{0.5cm} \bar{y}=(G^{\top})^{-1}\bar{y} \end{equation} Veamos que dado un $x\in R(A^{\top})$, entonces $x\in R(B^{\top})$, esto es, $R(A^{\top})\subset R(B^{\top})$. La implicación opuesta se prueba de igual forma: Considere $x\in R(B^{\top})$ entonces existe $\bar{y}\in \mathbb{F}^{m}$ tal que \begin{equation} x=B^{\top}\bar{y}=B^{\top}(G^{\top})^{-1}G^{\top}\bar{y}=(G^{-1}B)^{\top}y=A^{\top}y,\hspace{0.5cm} y=G^{\top}\bar{y}. \end{equation} Por lo tanto se ha mostrado que para cualquier $x\in R(B^{\top})$, entonces $x\in R(A^{\top})$, es decir, $R(B^{\top})\subset R(A^{\top})$. Por lo tanto $R(A^{\top})=R(B^{\top})$. Ahora si se asume que $R(A^{\top})=R(B^{\top})$. Esto quiere decir que cada fila de $A$ es una combinación lineal de las filas de $B$. Es to quiere decir, quiere decir que existe operaciones de Gauss sobre las filas de $A$ que transforma a $A$ en $B$. Y por lo tanto con eso se establece el resultado final del teorema 2.

Un argumento similar también dice que $E_{_A^{\top}}=E_{_{B^{\top}}}$ si y solo si $A^{\top}\stackrel{fila}{\longleftrightarrow} B^{\top}$. También se puede concluir que $E_{_{A^{\top}}}=E_{_{B^{\top}}}$ es equivalente a $R(A)=R(B)$ y esto también es equivalente a $N(A^{\top})=N(B^{\top})$.

Otro resultado con respecto a la transpuesta de $A$ dice:

Teorema 3.Para cada matriz $A\in \mathbb{F}^{n,m}$ se cumple que $\dim R(A)=\dim R(A^{\top})$.
Sabiendo que una matriz se dice que es de rango completo si y solo si $\dim R(A)=\min(m,n)$. Entonces podemos enunciar el siguiente teorema donde se establecen algunas de las relaciones principales del los cuatros espacios fundamentales.
Teorema 3. Si una matriz $A\in \mathbb{F}^{n,m}$ es de rango completo, entonces:
  • Si $m=n$, entonces: \[\dim R(A)=\dim R(A^{\top})=n=m\Leftrightarrow \{0_{_{\mathbb{F}^{m}}}\}=N(A)=N(A^{\top})\subset \mathbb{F}^{m};\]
  • Si $n < m$, entonces: \[\dim R(A)=\dim R(A^{\top})=n < m \Leftrightarrow \left\{\begin{array}{l} \{0_{_{\mathbb{F}^{m}}}\}\varsubsetneq N(A)\subset \mathbb{F}^{m},\\ \{0_{_{\mathbb{F}^{n}}}\}=N(A^{\top})\subset \mathbb{F}^{n};\end{array}\right.\]
  • Si $m>n$, entonces: \[ \dim R(A)=\dim R(A^{\top})= m < n \Leftrightarrow \left\{\begin{array}{l} \{0_{_{\mathbb{F}^{m}}}\}= N(A)\subset \mathbb{F}^{m},\\ \{0_{_{\mathbb{F}^{n}}}\}\varsubsetneq N(A^{\top})\subset \mathbb{F}^{n}; \end{array}\right.\]

Recordemos que el rango de una matriz $A$ es el número de pivotes columna de la matriz escalonada reducida $E_{_A}$ y que coincide con la dimensión de $R(A)$. Si una matriz $A$ de orden $m\times n$ tiene $rang(A)=n$, esto quiere decir dos cosas: Primero, $n\leq m$; y segundo, que cada columna de $E_{_A}$ tienes un pivote. Es decir, que no hay variables libres en la solución de la ecuación $Ax=0$, y así $x=0$ es la única solución. Por lo tanto $N(A)=0$. Por otro lado al estudiar $N(A^{\top})$ es necesario considerar dos casos $n=m$ o $n

Figura 2. Relaciones entre los cuatro espacios fundamentales de $A:\mathbb{F}^m\to\mathbb{F}^n$.

Para terminar enunciamos el siguiente resultado que relaciona las dimensiones de los espacios nulos y el rango de una transformación lineal en espacios vectoriales de dimensión finita. Este resultado se suele llamar Teorema de la Nulidad y el Rango, donde la nulidad de una transformación lineal es la dimensión de su espacio nulo, y el rango es la dimensión de su espacio de columna. Este resultado también se le conoce como el Teorema de la Dimensión.

Teorema 4. Para cada transformación lineal $T:V\to W$ entre espacios vectoriales $V$ y $W$ se cumple que \begin{equation} \dim N(T)+\dim R(T) = \dim V. \end{equation}

Si consideramos el producto punto de dos vectores $u, v\in \mathbb{F}^{m}$ como el producto matricial definido por $u\cdot v =u^{\top}v$. Podemos observar lo siguiente: \begin{equation} Ax=\left[\begin{array}{c} A_{_{1:}}x\\ \vdots \\ A_{_{m:}}x \end{array}\right] \end{equation} Note que las filas $A$ son las columnas de $A^{\top}$ por lo tanto $A^{\top}_{_{i:}}\cdot x= A_{_{i*}}x$ Si $x\in N(A)$, entonces se puede concluir que $N(A)$ es el complemento ortogonal de $R(A^{\top})$, por lo tanto $N(A)\cap R(A^{\top})=\emptyset$. Por otro lado si la dimensión de $\dim R(A)=r$ entonces la dimensión del nulo $\dim N(A) = m-r$. Y $\dim R(A^{\top})=r$ luego $\dim N(A^{\top})=n-r$. Note que $\dim N(A)\oplus R(A^{\top}) = m$ y $\dim N(A^{\top})\oplus R(A) = n$ es decir que $ N(A)\oplus R(A^{\top})\cong \mathbb{F}^{m}$ y $N(A^{\top})\oplus R(A) \cong \mathbb{F}^{n}$. Estas relaciones se resumen en la Figura 2.

Bueno, ya nos hemos extendido los suficiente por esta ocasión. Esperamos que aprovechen mucho este post.

Referencias

viernes, 7 de octubre de 2016

Diferencias finitas y el problema de valor de frontera de segundo orden

Se sabe que las ecuaciones diferenciales ordinarias permiten describir una gran variedad de problemas y fenómenos que surgen en las matemáticas, la física, la química, la economía, etc. Por lo cual hay un gran interés por encontrar soluciones a este tipo de problemas. Por lo tanto el objetivo de esta entrada es mostrar como encontrar soluciones aproximadas al problema de valor de frontera de segundo orden (PVF2) mediante diferencias finitas.

El problema que vamos a abordar en esta ocasión es el PVF2, definido como: \begin{equation}\left\{\begin{array}{l} y''=f(x,y,y'),\\ y(a)=\alpha,\\ y(b)=\beta. \end{array}\right. \end{equation} Lo anterior también se puede expresar de la siguiente forma \begin{equation}\label{eq02} y''(x)=p(x)y'(x)+q(x)y(x)+r(x),\,a\leq x\leq b,\, y(a)=\alpha,\, y(b)=\beta, \end{equation} en este caso se dice que el PVF2 es de tipo lineal, en caso contrario, se dice que el PVF2 es del tipo no lineal. Este tipo de problemas son muy comunes en la ciencias básicas o en ingeniería. Sin embargo, su principal dificultad radica en que la solución a dichos problemas no es sencilla. Debido a esto, resulta importante definir métodos numéricos para dar estimaciones numéricas de las soluciones. Existe una gran diversidad de métodos para hacer esto: el método del disparo, el método de elementos finitos, el método de Galerkin discontinuo, el método de diferencias finitas, etc. Aquí solo nos vamos a centrar en el método de diferencias finitas.

El método de diferencias finitas se basa en la idea de aproximar las derivadas que aparecen en el PVF2 mediante una partición del intervalo $[a, b]$. A menudo a estas particiones se le denominan mallas o grillas. Una vez que se hace estas aproximaciones, se sustituyen en la ecuación diferencial del PVF2 y posteriormente se resuelve un sistema de ecuaciones lineales o no lineales, esto depende de si el PVF2 es lineal o no.

El caso lineal

En el caso del problema $\eqref{eq02}$, se procede a discretizar el intervalo $[a,b]$ en una malla o grilla de puntos $a=x_{_{0}} < x_{_{1}} < \cdots < x_{_{n+1}}=b$, en la cual $x_{_{i+1}}-x_{_i}=h$, para todo $i=0,\dots, n$. En general los puntos no necesitan estar distribuidos de esta manera.

En los puntos de la forma $x_{i}=a+ih$ con $i=1,2,\dots,n$. Así la ecuación diferencial a aproximar es \begin{equation}\label{eq:03} y''(x_{_i})= p(x_{_{i}})y'(x_{_i})+q(x_{_i})y(x_{_i})+r(x_{_i}). \end{equation} Ahora empleando el polinomio de Taylor de orden 3 para $y(x)$, alrededor de $x_{_i}$, y se evalúa en $x_{_{i+1}}$ y $x_{_{i-1}}$, entonces, suponiendo que $y\in C^{4}[x_{_{i-1}},x_{_{i+1}}]$, se tiene \begin{equation}\label{eq:04}y(x_{_{i+1}})=y(x_{_i}+h)=y(x_{_i})+hy^{(1)}(x_{_i})+\frac{h^2}{2}y^{(2)}(x_{_i})+\frac{h^3}{6}y_{_i}^{(3)}(x_{_i})+\frac{h^4}{24}y^{(4)}(\xi_{_i}^{+}) \end{equation} para algún $\xi_{_i}^{+}$ en $(x_{_i},x_{_{i+1}})$, y \begin{equation}\label{eq:05}y(x_{_{i-1}})=y(x_{_i}-h)=y(x_{_i})-hy^{(1)}(x_{_i})+\frac{h^2}{2}y^{(2)}(x_{_i})-\frac{h^3}{6}y_{_i}^{(3)}(x_{_i})+\frac{h^4}{24}y^{(4)}(\xi_{_i}^{^-}) \end{equation} para algún $\xi^{^-}_{_i}$ en $(x_{_{i-1}},x_{_i})$. Si se suman adecuadamente $\eqref{eq:04}$ y $\eqref{eq:05}$ se obtiene \[y(x_{_{i+1}})+y(x_{_{i-1}})=2y(x_{_i})+h^2y^{(2)}(x_{_i})+\frac{h^4}{24}\left(y^{(4)}(\xi_{_i}^{+})+y^{(4)}(\xi_{_i}^{^-})\right),\] y al despejar $y^{(2)}(x_{_i}),$ \[y^{(2)}(x_{_i})=\frac{1}{h^2}\left(y(x_{_{i+1}})-2y(x_{_i})+y(x_{_{i-1}})\right)-\frac{h^4}{24}\left(y^{(4)}(\xi_{_i}^{+})+y^{(4)}(\xi_{_i}^{^-})\right).\] Ahora, aplicando el teorema de valor intermedio, se tiene que \begin{equation}\label{eq:06} y^{(2)}(x_{_i})=\frac{1}{h^2}\left(y(x_{_{i+1}})-2y(x_{_i})+y(x_{_{i-1}})\right)-\frac{h^4}{12}y^{(4)}(\xi_{_i}), \end{equation} para algún $\xi_{_i}\in (x_{_{i-1}},x_{_{i+1}})$. A este esquema se le llama fórmula de diferencias centradas para $y^{(2)}(x_{_i})$. Ahora si a $\eqref{eq:04}$ se le resta $\eqref{eq:05}$. De manera semejante se obtiene una fórmula de este tipo para $y^{(1)}(x_{_i})$, la cual viene dada por \begin{equation}\label{eq:07} y^{(1)}(x_{_i})=\frac{1}{2h}\left(y(x_{_{i+1}})-y(x_{_{i-1}})\right)-\frac{h^3}{3}y^{(3)}(\eta_{_i}) \end{equation} para algún $\eta_{_i}$ en $(x_{_{i-1}},x_{_{i+1}})$.

Ahora si se sustituyen $\eqref{eq:06}$ y $\eqref{eq:07}$ en $\eqref{eq:03}$, entonces se obtiene la igualdad \[\frac{1}{h^2}\left(y(x_{_{i+1}})-2y(x_{_i})+y(x_{_{i-1}})\right)=\frac{p(x_{_i})}{2h}\left(y(x_{_{i+1}})-y(x_{_{i-1}})\right)-\frac{p(x{_i})h^3}{3}y^{(3)}(\eta_{_i})+\frac{h^4}{24}\left(y^{(4)}(\xi_{_i}^{+})+y^{(4)}(\xi_{_i}^{^-})\right)\] Por lo tanto hay un error del orden $\mathcal{O}(h^3)$, con $y_{_i}=y(x_{_i})$, $y'_{_i}=y'(x_{_i})$, $y''_{_i}=y'(x_{_i})$, $p_{_i}=p(x_{_i})$, $q_{_i}=q(x_{_i})$ y $r_{_i}=r(x_{_i})$, podemos escribir simplemente \[\frac{1}{h^2}\left(y_{_{i+1}}-2y_{_i}+y_{_{i-1}}\right)=\frac{p_{_i}}{2h}\left(y_{_{i+1}}-y_{_{i-1}}\right)+q_{_i}y_{_i}+r_{_i}+\mathcal{O}(h^{3})\] y por lo tanto de esta manera obtenemos el conjunto de ecuaciones \[\frac{1}{h^2}\left(y_{_{i+1}}-2y_{_i}+y_{_{i-1}}\right)-\frac{p_{_i}}{2h}\left(y_{_{i+1}}+y_{_{i-1}}\right)-q_{_i}y_{_i}=r_{_i},\] para cada $i=1,\dots, n$, junto con $y_{_0}=\alpha$ y $y_{_{n+1}}=\beta$. Reescribiendo la ecuación anterior, se tiene que \[-\left(1+\frac{h}{2}p_{_i}\right)y_{_{i-1}}+\left(2+h^2q_{_i}\right)y_{_i}-\left(1-\frac{h}{2}p_{_i}\right)y_{_{i+1}}=-h^2r_{_i}\] El sistema de ecuaciones resultante se puede expresar en la forma $Ay=b$, donde, \[A=\left[\begin{array}{cccccc} \gamma_{_1} & \delta_{_1} & 0 & \cdots & \cdots & 0 \\ \epsilon_{_2} & \gamma_{_2} & \delta_{_2} & \ddots & & \vdots \\ 0 & \epsilon_{_3} & \gamma_{_3} & \delta_{_3} & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & 0 \\ \vdots & & \ddots & \ddots & \ddots & \delta_{_{n-1}} \\ 0 & \dots & \dots & 0 & \epsilon_{_n} & \gamma_{_n} \\ \end{array}\right],\,\,\begin{array}{ll} \gamma_{_i}=2+h^2q_{_i}, & 1\leq i\leq n,\\ \delta_{_i}=-1+\frac{h}{2}p_{_i}, & 1\leq i\leq n-1,\\ \epsilon_{_i}=-1-\frac{h}{2}p_{_i}, & 2\leq i\leq n,\end{array} \] \[\hat{y}=\left[\begin{array}{c} y_{_1} \\ y_{_2} \\ \vdots \\ y_{_n} \end{array}\right] \mbox{ e } b=\left[\begin{array}{c} -h^2r_{_i}+(1+\frac{h}{2}p_{_i})y_{_0} \\ -h^2r_{_{2}} \\ -h^2r_{_{3}} \\ \vdots \\ -h^2r_{_{n-1}}\\ -h^2r_{_i}+(1+\frac{h}{2}p_{_i})y_{_{n+1}} \end{array}\right].\]

Para revolver este sistema lineal de ecuaciones recomendamos leer factorización LU y para la reconstrucción aproximada de la curva solución recomendamos leer sobre interpolación polinómica de Lagrange.

¿Bajo qué condiciones el sistema lineal $Ay=b$ tiene única solución? La respuesta nos la da el siguiente teorema.

Teorema. Suponga que $p(x)$, $q(x)$ y $r(x)$ son continuas en $[a, b]$. Si $q(x)\geq 0$ en $[a,b]$, entonces el sistema triagonal $Ay=b$ tiene solución única siempre y cuando $h\leq \frac{2}{L}$, donde $L=\max_{a\leq x\leq b}|p(x)|$.

En cuando a los errores el siguiente teorema nos ayuda a deducir las cotas para el error de la aproximación de las ecuaciones diferenciales mediante este método de diferencias finitas de segundo orden.

Teorema Sea $Q^*$ tal que $q(x)>Q^*>0$, para $a\leq x\leq b$, $L=\max_{a\leq x\leq b}|p(x)|$, $M_{_3}=\max_{a\leq x\leq b} |y^{(3)}(x)|$ y $M_{_4}=\max_{a\leq x\leq b}|y^{(4)}(x)|$. Entonces, para $i=0,1,\dots, n+1$, se cumple: \[|\hat{y}_{_i}-y(x_{_i})|\leq h^{3}\left(\frac{M_{_4}+2LM_{_3}}{12Q^*}\right).\]

Este teorema nos dice que la solución de diferencias finitas converge a la solución exacta cuando $h\to 0$ y el error es $\mathcal{O}(h^{3})$. Cuando en la ecuación diferencial lineal, $p(x)=0$, entonces se puede demostrar que el error cometido es del orden de $\mathcal{O}(h^4)$.

El caso no lineal

Para el problema no lineal general con valor de frontera: \[y''=f(x, y,y'),\,a\leq x\leq b,\, y(a)=\alpha,\, y(b)=\beta,\] el método de diferencias finitas es similar al caso lineal. Sin embargo, el sistema de ecuaciones no será lineal y, por lo tanto, se requiere un proceso iterativo para resolverlo. Nuevamente, se procede a discretizar el intervalo $[a,b]$ en una partición de puntos $a=x_{_0}\leq x_{_1}\leq\cdots\leq x_{_{n+1}}=b$, tal que $x_{_{i+1}}-x_{_i}=h$.

En los puntos interiores de la partición $x_{_i}=a+ih$, $i=1,2,\dots,n$, la ecuación diferencial a aproximar es: \begin{equation}\label{eq:08} y''(x_{_i})=f(x_{_i},y(x_{_i}),y'(x_{_i})). \end{equation} Luego sustituyendo $\eqref{eq:04}$ y $\eqref{eq:05}$ en $\eqref{eq:08}$ se obtiene: \[\frac{1}{h^2}\left(y(x_{_{i+1}})-2y(x_{_i})+y(x_{_{i-1}})\right)=f\left(x,y(x_{_i}),\frac{1}{2h}\left(y(x_{_{i+1}})-y(x_{_{i-1}})\right)-\frac{h^3}{3}y^{(3)}(\eta_{_i})\right)+\frac{h^4}{12}y^{(4)}(\xi_{_i})\] para algún $\xi_{_i}$ y $\eta_{_i}$ en el intervalo $(x_{_{i-1}},x_{_{i+1}})$, con $i=1,2,\dots,n$. El método de diferencias finitas se obtiene empleando esta ecuación junto con las condiciones finitas de frontera $y(a)=\alpha$ y $y(b)=\beta$. Por lo tanto si nuevamente se hace $y_{_i}=y(x_{_i})$, $y'_{_i}=y'(x_{_i})$, $y''_{_i}=y'(x_{_i})$, $p_{_i}=p(x_{_i})$, $q_{_i}=q(x_{_i})$ y $r_{_i}=r(x_{_i})$ se obtiene el sistema de ecuaciones \[\frac{y_{_{i+1}}-2y_{_i}+y_{_{i-1}}}{h^2}-f\left(x_{_i},y_{_i},\frac{y_{_{i+1}}-y_{_{i-1}}}{2h}\right)=0,\] para cada $i=1,2,\dots,n$ junto con $y_{_0}=\alpha$ y $y_{_{n+1}}=\beta$. El sistema no lineal de $n\times n$ obtenido con este método $F_{_i}(y_{_1},\dots,y_{_n})=0$ con \[F_{_i}(y_{_1},\dots,y_{_n})=y_{_{i+1}}-2y_{_i}+y_{_{i-1}}-h^2f\left(x_{_i},y_{_i},\frac{y_{_{i+1}}-y_{_{i-1}}}{2h}\right)\] tiene solución única si $h<\frac{2}{L}$, donde $L$ es tal que $\left|\frac{\partial f}{\partial y'}(x,y,y')\right|\leq L$.

Para resolver este sistema de ecuaciones usaremos el método de Newton. Es consiste es suponer una solución $y^{(0)}=(y^{(0)}_{_1},\dots, y^{(0)}_{_{n}})^{\top}$, la iteración con este método es: \[y^{(j+1)}=y^{(j)}-\pmb{J}^{-1}(y^{(j)})F(y^{(j)})\] donde $\pmb{J}(y^{(j)})$ es la matriz jacobiana. En la práctica es más conveniente a efectos computacionales iterar resolviendo las ecuaciones: \[\pmb{J}(y^{(j)})(y^{(j+1)}-y^{(j)})=-F(y^{(j)}).\] En cuanto al jacobiano $\pmb{J}_{rs}(y^{(j)})=\frac{\partial F_{_r}}{\partial y_{_s}}$, se tiene: \[J_{_{kk}}=-2-h^{2}\frac{\partial f}{\partial y}\left(x_{_k},y_{_k},\frac{y_{_{k+1}}-y_{_{k-1}}}{2h}\right),\] \[J_{_{k(k+1)}}=1-\frac{h}{2}\frac{\partial f}{\partial y'}\left(x_{_k},y_{_k},\frac{y_{_{k+1}}-y_{_{k-1}}}{2h}\right),\] \[J_{_{(k+1)k}}=1+\frac{h}{2}\frac{\partial f}{\partial y'}\left(x_{_k},y_{_k},\frac{y_{_{k+1}}-y_{_{k-1}}}{2h}\right),\] con $k=1,\dots, n$. El resto de las entradas de la matris $\pmb{J}$ son nulos. En particular, como todo los elementos $|r-s|>1$ son nunos, entonces $J$ es tridiagonal. En cada iteración del sistema lineal: \[\pmb{J}(y^{(j)})v=-F(y^{(j)}),\] con $v=(v_{_1},\dots,v_{_n})^{\top}$, se requiere obtener $v$, porque \[y^{(j+1)}=y^{(j)}+v.\] Puesto que $\pmb{J}$ es triagonal, se puede utilizar métodos para resolver ecuaciones lineales. Nosotros, para efecto de implementación, en un próximo post usaremos la factorización LU para resolver este tipo de problemas. Para finalizar cabe mencionar que nuevamente la solución que se encuentra con este método se aproxima a la solución exacta con un orden de aproximación del orden $\mathcal{O}(h^{3})$.

Referencias

viernes, 30 de septiembre de 2016

Factorización LU

En esta ocasión estudiamos la factorización $LU$ y la eficiencia que promueve en la solución de sistemas de ecuaciones con muchas incógnitas. Como caso particular, se aplica en el método de diferencias finitas que tiene lugar cuando se quiere resolver algunas ecuaciones diferenciales lineales de orden 2 con valores en la frontera. Esperamos que sea de su agrado.

Una de las factorizaciones más elegante e importante del Álgebra Lineal es la factorización $A = LU$ Cuando se habla de la posibilidad de resolver un sistema de ecuaciones por un método más eficiente que la eliminación Gaussiana cuesta creerlo! Pero existe dicho método y esto es gracias a la factorización de la matriz $A$ como producto de dos matrices $L$ (triangular inferior) y $U$ (triangular superior).

Para entender la factorización $LU$ de una matriz $A$ a la que se puede aplicar eliminación Gaussiana sin encontrar ceros en la diagonal principal y sólo con operaciones elementales de fila que consisten en sustituir una fila, por el resultado de sumarle a ésta un múltiplo escalar de otra, una estrategia es pensar en dicha operación elemental como una premultiplicación por una matriz, en este caso, una matriz elemental triangular inferior.

Suponga que durante el proceso de eliminación Gaussiana para la matriz $$A = \left [ \begin{array}{ccc} 2 & 2& 2\\ 4& 7& 7\\6 &18 &22 \end{array} \right ] $$ A las operaciones elementales de fila:

  1. Sustituir fila $2$ por fila $2$ $-$ $2$ veces fila $1$.
  2. Sustituir fila $3$ por fila $3$ $-$ $3$ veces fila $1$.

Se les interpreta como un premultiplicación por la matriz $T = \left [ \begin{array}{ccc} 1 &0 & 0\\-2 & 1 & 0\\ -3& 0 & 1 \end{array} \right ] $. El cálculo directo de $T$ por $A$ es: $$\left [ \begin{array}{ccc} 1 &0 & 0\\-2 & 1 & 0\\ -3& 0 & 1 \end{array} \right ] \left [ \begin{array}{ccc} 2 & 2& 2\\ 4& 7& 7\\6 &18 &22 \end{array} \right ] = \left [ \begin{array}{ccc} 2 & 2& 2\\ 0& 3& 3\\0 &12 &16 \end{array} \right ] $$ El resultado obtenido es el mismo que si se aplican las operaciones elementales de fila.
Pero, cómo se construye en general dicha matriz $T$?, además, es $T$ no singular (invertible)?. Esas preguntas se van a responder y se vará la relación que tienen con la factorización $LU$.

Para ver que lo anterior (en el ejemplo ilustrativo inicial) sigue pasando en general, se define una matriz elemental triangular interior de orden $n$ como $$T _{k} = I - c _{k} e _{k} ^{T} $$ donde $c _{k} $ es columna con $k$ ceros al principio, es decir de la forma $c _{k} = \left [ \begin{array}{c} O _{k} \\ -\mu _{k + 1 } \\ \vdots \\ -\mu _{n} \end{array} \right ] $ y $e _{k} = I _{*k} $ (k-ésima columna de la identidad $I _{n \times n}$).
La forma que tiene $T _{k} $ es similar a la identidad $I$, excepto que debajo del 1 en la $k$-ésima columna tiene los valores $-\mu _{k + 1 },..., -\mu _{n} $, de la siguiente forma: $$T _{k} = \left [ \begin{array}{cccccc} 1 _{1} & & & & & \\ & \ddots & & & & \\ & & 1 _{k} & & & \\ & &-\mu _{k+1} &1 _{k+1} & &\\& &\vdots & &\ddots &\\& & -\mu _{n} & & & 1 _{n} \end{array} \right ] $$ Considerando los espacios vacíos como ceros, note que tiene la misma forma que la matriz $T$ en el ejemplo ilustrativo inicial.

Se puede sospechar que $T$ sea invertible, pero lo que ahora es sorprendente es que $T _{k} ^{-1} = I + c _{k} e _{k} ^{T}$, lo que se puede ver por multiplicación directa $$(I + c _{k} e _{k} ^{T} )(I - c _{k} e _{k} ^{T} ) = I - c _{k} e _{k} ^{T} + c _{k} e _{k} ^{T} - c _{k} e _{k} ^{T} c _{k} e _{k} ^{T} $$ No queda claro aún es por qué el término $c _{k} e _{k} ^{T} c _{k} e _{k} ^{T} $ es la matriz nula. Para ello, note que $$c _{k} e _{k} ^{T} c _{k} e _{k} ^{T} = c _{k} ( e _{k} ^{T} c _{k} )e _{k} ^{T} = c _{k} 0 e _{k} ^{T} $$ Se ha sustituido $e _{k} ^{T} c _{k} $ por el n\'umero $0$ puesto que recuerde que $e _{k} ^{T} $ solo tiene un 1 en la posición $k$, mientras que $c _{k}$ tiene un $0$ en ese mimo lugar.
Se puede decir que $T _{k} ^{- 1} $ también es una matriz elemental triangular inferior. $T _{k} ^{-1} $ es muy similar a $T _{k} $, los números fuera de la diagonal de 1's cambian de signo.

Todo lo anterior da una idea de qué es lo que se hace al poner ceros bajo la diagonal en el proceso de escalonar $A$ para llegar a una matriz $U$ triangular superior. En realidad lo que se hace es premultiplicar por $n - 1$ matrices elementales triangulares inferiores, es decir: $$A \to T _{1} A \to T _{2} T _{1} A \to \cdots \to T_{n - 1 } T _{n - 2} T _{n - 3} \cdots T _{2} T _{1} A = U $$ Y al premultiplicar por las inversas $$T _{n - 2}\cdots T _{2} T _{1} A = T _{n - 1 } ^{-1} U \ \to \ T _{n - 3}\cdots T _{2} T _{1} A = T _{n - 2} ^{-1} T _{n - 1 } ^{-1} U \ \to \ A = T _{1} ^{-1} \cdots T _{n - 1} ^{-1} U $$ Asi que $$A = LU$$ con $L = T _{1} ^{-1} \cdots T _{n - 1} ^{-1} $.

Ahora, $L =T _{1} ^{-1} \cdots T _{n - 1} ^{-1} = (I + c _{1} e _{1} ^{T} )(I + c _{2} e _{2} ^{T} )(I + c _{3} e _{3} ^{T} ) \cdots (I + c _{n-1} e _{n-1} ^{T} ) $. Por multiplicación directa se puede notar que $$L = I + c _{1} e _{1} ^{T} + c _{2} e _{2} ^{T} + c _{3} e _{3} ^{T} + \cdots+ c _{n-1} e _{n-1} ^{T} $$ Puesto que en el producto de las matrices del tipo $I + c _{k} e _{k} ^{T} $ se cancelan muchos térmidos dado que los factores del tipo $ e _{k} ^{T} c _{j} $ son el número cero cuando $k < j$. (Recuerde que $e _{k} ^{T}$ solo tiene un 1 en la posición $k$ mientras que $c _{j} $ tiene posiblemente números diferentes de cero sólo hasta la posición $j + 1$).

Para ver cómo es $L$ se debe percatar de la forma de la suma $$I + c _{1} e _{1} ^{T} + c _{2} e _{2} ^{T}+ c _{3} e _{3} ^{T} + \cdots+ c _{n-1} e _{n-1} ^{T} $$ $$I +\left [ \begin{array}{cccccc} & & & & & \\ l _{21} & & & & & \\ l _{31}& & & & & \\ l _{41} & && & &\\ \vdots& & & & &\\ l _{n1} & & & & & \end{array} \right ]+ \left [ \begin{array}{cccccc} & & & & & \\ & & & & & \\&l _{32} & & & & \\ & l _{42} && & &\\ & \vdots & & & &\\ & l _{n2}& & & & \end{array} \right ] + \cdots + \left [ \begin{array}{cccccc} & & & & & \\ & & & & & \\& & & & & \\ & && & &\\ & & & & &\\ & & & & l _{n,n - 1} & \end{array} \right ] $$ Es decir $$L = \left [ \begin{array}{cccccc} 1 & & & & & \\ l _{21} & 1& & & & \\ l _{31}&l _{32} & 1 & & & \\ l _{41} & l _{42} &l _{43} & \ddots & &\\\vdots & \vdots & \vdots & &1 &\\ l _{n1} & l _{n2}& l _{n3}& & l _{n,n - 1} & 1 \end{array} \right ] $$

A modo de conclusión se tiene que:

  1. $L$ es triangular inferior y $U$ es triangular superior. Además, $l _{ii} = 1 $ y $u _{ii} \neq 0 $ para cada $i = 1,2,..., n$.
  2. Bajo la diagonal de $L$, la entrada $l _{ij} $ es el múltiplo de la fila $j$ que es sustraida de la fila $i$ con el objetivo de anular la posición $ij$ durante la eliminación Gaussiana.
  3. $U$ es el resultado final al escalonar $A$. Además, las matrices $L$ y $U$ son únicas.

En el ejemplo inicial con $$A = \left [ \begin{array}{ccc} 2 & 2& 2\\ 4& 7& 7\\6 &18 &22 \end{array} \right ] $$ se puede terminar la eliminación y notar que $$L= \left [ \begin{array}{ccc} 1 & 0& 0\\ 2& 1& 0\\3 &4 &1 \end{array} \right ] \ \ \ y \ \ \ A = \left [ \begin{array}{ccc} 2 & 2& 2\\ 0& 3& 3\\0 &0 &4 \end{array} \right ] $$ Es importante percatarse de algo genial, y es que las entradas que caracterizan a $L$ (sólo las de la parte inferior izquierda) y las que caracterizan a $U$ (sólo las de la parte superior derecha), suman tantas entradas como las que inicialmente tiene $A$, lo que significa que se necesita la misma capacidad de memoria para almacenar sólo a $A$ que para almacenar a las dos matrices $L$ y $U$. Ahora se da paso al análisis del proceso con el que se resuelve un sistema de ecuaciones $Ax = b$ usando $L$ y $U$ en lugar de $A$. Además, se resalta la importancia de tener las componentes $L$ y $U$ para resolver familias de sistemas $Ax = \overline{b}$ que tengan la misma matriz de coeficientes $A$ y cualquier otro vector de términos independientes $\overline{b}$.
Suponga que ya se tienen $L$ y $U$. El sistema $Ax = b$ se puede ver como $LUx = b$ y si se asume $Ux = y$ entonces el sistema inicial $Ax = b$ se puede interpretar como dos sistemas triangulares $Ly = b$ y $Ux = y $.
Se procede con el sistema $Ly = b$ en donde el vector $y$ de incógnitas se halla por sustitución directa $$ \left [ \begin{array}{cccccc} 1 & & & & & \\ l _{21} & 1& & & & \\ l _{31}&l _{32} & 1 & & & \\ l _{41} & l _{42} &l _{43} & \ddots & &\\\vdots & \vdots & \vdots & &1 &\\ l _{n1} & l _{n2}& l _{n3}& & l _{n,n - 1} & 1 \end{array} \right ] \left [ \begin{array}{c} y _{1} \\ y _{2} \\ y _{3} \\ y _{4}\\ \vdots \\ y _{n} \end{array} \right ] = \left [ \begin{array}{c} b _{1} \\ b _{2} \\ b _{3} \\ b _{4} \\ \vdots \\ b _{n} \end{array} \right ] $$ Es decir, $y _{1} = b _{1}$, $\ \ y _{2} = b _{2} - l _{21} y _{1}$, $\ \ y _{3} = b _{3} - l _{31} y _{1} $, etc.
En otras palabras, $y _{1} = b _{1} \ \ $ y $\ \ y _{i} = b _{i} - \displaystyle \sum_{k = 1 } ^{i - 1 } l _{ik} y _{k} $ para $i = 2,3,...,n$.
Cuando ya se tienen las entradas del vector $y$, se monta el sistema triangular superior $Ux = y$ que se resuelve usando sustitución regresiva comenzando con $x _{n} = \frac{y _{n} }{u _{n} } $ y tomando $$ x _{i} = \frac{1 }{u _{ii} } \left ( y _{i} - \displaystyle \sum_{k = i + 1 } ^{n} u _{ik} x _{k} \right ) \ \ \ para \ \ \ i = n - 1 , n - 2, ..., 1.$$ Cabe notar que con éste método de sustituciones usando $L$ y $U$, sólo se requieren $n ^{2} $ multiplicación y $n ^{2} - n$ adiciones. Por otro lado, para la eliminación en el sistema aumentado $A | b$ se requieren alrededor de $n ^{3} /3$ multiplicaciones, lo que indica que para $n$ grande, usar la factorización $LU$ genera evidente eficiencia. Además, en cualquier otro sistema $Ax = \overline{b}$ con la misma matriz $A$ se pueden aplicar los mismo algoritmos eficientes de sustitución.
Por ejemplo, en el importante problema (atacado con método de diferencias finitas) para la ecuación diferencial $\frac{d ^{2} }{dt ^{2} } x = f(t) ,\ \ \ con \ \ \ 0 \leq t \leq 1, \ \ x(0) = x(1) = 0$, sale a flote la matriz de coeficientes $$A = \left [ \begin{array}{ccccccc} 2& -1&& && & \\ -1& 2&-1& && & \\ &-1 &2& -1&& &\\ & &&\ddots && &-1 \\ & && &&-1 &2 \end{array} \right ] $$ Cuyos factores son $$L = \left [ \begin{array}{ccccccc} 1& && && & \\ \frac{-1}{2} & 1& && & & \\ &\frac{-2}{3} &1& && &\\ & &\frac{-3}{4}&\ddots && & \\ & && && \frac{-(n-1)}{n}&1 \end{array} \right ] \ \ \ y \ \ \ U = \left [ \begin{array}{ccccccc} 2& -1&& && & \\ &\frac{3}{2} &-1& && & \\ & &\frac{4}{3}& -1&& &\\ & &&\ddots && &-1 \\ & && && &\frac{n + 1 }{n} \end{array} \right ] $$ Y cuyas fórmulas de sustitución recursiva para los respectivos sistemas triangulares $Ly = b$ y $Ux = y$ son respectivamente $$y _{1} = b _{1} , y _{k} = b _{k} + \frac{k - 1 }{k} y _{k - 1 }, \ para \ k = 2,...,n. $$ $$x _{n} = \frac{n}{n + 1 }y _{n}, \ x _{k} = \frac{k}{k + 1 } (y _{k} + x _{k+1} ),\ para \ \ k = n - 1, n - 2,...,1. $$

Se recomienda ir a los trabajos de Carl D Meyer y de Gilbert Strang

Referencias

viernes, 23 de septiembre de 2016

Números de coma o punto flotante y sistemas líneales

En esta ocasión estudiaremos algunos aspectos generales de los números de punto flotante y algunas dificultades que se presentan cuando se realizan cálculos elementales en la resolución de sistemas lineales por eliminación Gaussiana y Gauss - Jordan. Esperamos que lo disfruten mucho.

Los números de coma o punto flotante son un conjunto finito de número racionales que se emplean para representar números reales empleando computadoras. Existen tipos diferentes de números de punto flotantes, todos ellos se caracterizan porque tienen un número finito de dígitos cuando se escriben en una base particular. La necesidad de estos números es precisamente porque las computadoras solo pueden representar los números reales con un conjunto finito de dígitos.

Definición. Un número racional $x$ es un número de punto flotante en base $b\in \mathbb{N}$, de precisión $p\in \mathbb{N}$, con rango de exponencial $N\in \mathbb{Z}$ si y sólo si existen enteros $d_{_i}$, para $i=1,\dots, b-1$ y $d_{_1}\neq 0$ tal que $x$ tiene la forma \begin{equation}\label{eq:01} x=\pm 0.d_{_1}\cdots d_{_p}\times b^{n},\hbox{ con } -N\leq n\leq - N. \end{equation} Se denota por $\mathbb{F}_{p, b, N}$ el conjunto de todos los números de punto flotante con precisión $p$, base $b$ y rango exponencial $N$.

Lo primero que se puede observar es que el conjunto $\mathbb{F}_{p, b, N}$ es un conjunto finito y que no hay una distribución homogénea de sus elementos. Por ejemplo si consideramos el conjunto de número flotantes de precisión $1$, base $10$ y rango exponencial $10$, como el lector podrá notar, es fácil listar los elementos positivos de este conjunto. En efecto sus elementos son, \[\{0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09\}=\{0.i\times 10^{-1}\}_{i=1}^{9},\] \[\{0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9\}=\{0.i\times 10^{0}\}_{i=1}^{9},\] \[\{1, 2, 3, 4, 5, 6, 7, 8, 9\}=\{0.i\times 10^{0}\}_{i=1}^{9}.\]

Como se puede observar los elementos de $\mathbb{F}_{1, 10, 1}$ no se encuentran homogéneamente distribuidos sobre el intervalo $[0,10]$. Esta distribución irregular afecta notablemente la posibilidad de hacer cálculos de adicción y multiplicación entre números pequeños y números grandes o entre solamente números grandes. Observe por ejemplo que cuando sumamos $0.01$ con $9$ obtenemos $9.01$, que no pertenece a $\mathbb{F}_{1, 10, 1}$, igualmente se puede verificar que el producto de $8$ por $9$ tampoco pertenece.


Figura 1. Distribución no homogénea de los elementos de $\mathbb{F}_{1, 10, 1}$.

En general los conjuntos de números flotantes $\mathbb{F}_{p, b, N}\subset\mathbb{R}$ no son cerrados bajo la suma o la multiplicación. Por lo tanto al tratar de representar números reales en algún conjunto $\mathbb{F}_{p, b, N}\subset\mathbb{R}$ habrá ciertos cálculos aritméticos que no serán permitidos. Una manera de realizar cálculos con números reales en $\mathbb{F}_{p, b, N}\subset\mathbb{R}$ es primero proyectar lo números reales en los números de punto flotante y que luego se realice el cálculo. Ya que el resultado podría no estar en el conjunto $\mathbb{F}_{_{p, b, N}}\subset\mathbb{R}$, uno debe proyectar de nuevo el resultado en $\mathbb{F}_{p, b, N}\subset\mathbb{R}$. La acción de proyectar un número real en el conjunto de números flotantes es llamada redondear un número real. Existe muchas forma de hacer esto. A continuación presentamos la forma más común de hacerlo.

Definición. Sea $X_{_N}=\max \mathbb{F}_{p, b, N}$, entonces la función de redondeo se define como una aplicación $fl:\mathbb{R}\cap [-X_{_N},X_{_N}]\to \mathbb{F}_{p, b, N}$ definida como: Dado un número real $x\in \mathbb{R}\cap [-X_{_N},X_{_N}]$, con $x=\pm 0.d_{_1}\cdots d_{_p}d_{_{p+1}}\cdots \times b^{n}$ y $-N\leq n\leq N$, se tiene que \begin{equation} fl(x)=\left\{\begin{array}{cl} \pm 0.d_{_1}\cdots d_{_p}\times b^n & \hbox{ si }d_{_{p+1}}<\frac{b}{2}, \\ \pm (0.d_{_1}\cdots d_{_p}+b^{-p})\times b^n & \hbox{ si }d_{_{p+1}}\geq \frac{b}{2}.\end{array}\right. \end{equation}

Por ejemplo si consideramos el conjunto de números flotantes $\mathbb{F}_{3,10,3}$, entonces algunas proyecciones de los números reales serían $fl(0.2103\times 10^{3})=0.210\times 10^3$, $fl(0.21037\times 10^{3})=0.210\times 10^3$, $fl(0.2105\times 10^{3})=0.211\times 10^3$ y $fl(0.2107\times 10^{3})=0.211\times 10^3$. Observe como la función $fl$ no es una función inyectiva, y por lo tanto, diferentes números reales pueden ser representados con el mismo número flotante.

Otra consecuencia que cabe mencionar, es que la aritmética de los número reales cambia notablemente. Es decir, si $x,y\in \mathbb{R}$ entonces la suma de reales sobre un conjunto $\mathbb{F}_{p, b, N}$ queda definida por \begin{equation} x+_{f}y=fl(fl(x)+fl(y)) \end{equation} y la multiplicación por \begin{equation} x\cdot_{f}y=fl(fl(x)\cdot fl(y)) \end{equation} Aquí las operaciones $+_{f}$ y $\cdot_{f}$ son la suma y la multiplicación en el conjunto de coma flotante $\mathbb{F}_{p, b, N}$. Ésta operaciones son diferentes de la suma ($+$) y la multiplicación ($\cdot$) usual de números reales. En efecto, se tiene la siguiente proposición.

Proposición. Sea $\mathbb{F}_{p, b, N}$ entonces siempre existen $x,y\in \mathbb{R}$ tales que \begin{equation} fl(fl(x)+fl(y))\neq fl(x+y),\hspace{0.5cm} fl(fl(x)\cdot fl(y))\neq fl(xy).\end{equation}

Estamos seguro que el lector puede encontrar algunos ejemplos sencillos que verifican la proposición anterior.

Como se ha visto las operaciones aritméticas usuales no son posible en general en los conjuntos $\mathbb{F}_{p, b, N}$, debido a que estos conjuntos no son cerrados bajo estas operaciones. Por lo que ha definido de una manera un poco artificial las operaciones $+_{f}$ y $\cdot_{f}$. A continuación veremos algunos efectos de estas operaciones cuando se trata de resolver sistemas de ecuaciones lineales en algún $\mathbb{F}_{p, b, N}$.

Considera el conjunto de números flotantes $\mathbb{F}_{3,10,3}$ para resolver el siguiente sistema de ecuaciones lineales, \[ \begin{array}{rcc} 5x_{_1}+x_{_2} &=& 6,\\ 9.34x_{_1}+1.57x_{_2} &=& 11. \end{array} \] Observe que al resolver este sistemas $\mathbb{R}$ sin las limitaciones de $\mathbb{F}_{3,10,3}$ encontramos que las soluciones son $x_{_1}=x_{_2}=1$. Veamos ahora que ocurre cuando empleamos $\mathbb{F}_{3,10,3}$ con las operaciones $+_{f}$ y $\cdot_{f}$ al realizar Gauss - Jordan sobre este sistema de ecuaciones. Inicialmente tenemos el sistema matricial, \[\left[\begin{array}{cc|c} 5 & 1 & 6\\ 9.43 & 1.57 & 11\end{array}\right]\] para cambiar la segunda fila podemos hacer la siguiente operación sobre la filas, empleado la suma y la multiplicación de redondeo en $\mathbb{F}_{3,10,3}$. \[\hat{a}_{_{2i}}=fl(a_{_{21}}-fl\left(fl(a_{_{1i}})\cdot fl\left(\frac{9.53}{5}\right)\right)\hbox{ con } i=1,2,3\] donde $\hat{a}_{_{2i}}$ son las entradas de la nueva fila $2$ y $a_{_{1i}}$ son las entradas de la fila $1$ en el sistema inicial. Al hacer estas operaciones se tiene que \[\left[\begin{array}{cc|c} 5 & 1 & 6\\ 9.43 & 1.57 & 11\end{array}\right]\rightarrow \left[\begin{array}{cc|c} 5 & 1 & 6\\ 0.02 & -0.32 & -0.3\end{array}\right].\] En este caso no es posible continuar con el método de Gauss - Jordan al menos que $\hat{a}_{_{21}}=0$. Pero esto no es posible. Si introducimos la modificación de que $\hat{a}_{_{21}}=0$, se tiene para nuestro ejemplo que \[\left[\begin{array}{cc|c} 5 & 1 & 6\\ 9.43 & 1.57 & 11\end{array}\right]\rightarrow \left[\begin{array}{cc|c} 5 & 1 & 6\\ 0 & -0.32 & -0.3\end{array}\right].\] Es importante notar que aquí no se ha redondeado el error, pues el valor $0.02$ pertenece a $\mathbb{F}_{3,10,3}$. Lo que se ha hecho es una modificación al sistema para poder continuar con el proceso de Gauss - Jordan. Si completamos el proceso bajo las operaciones $+_{f}$ y $\cdot_{f}$ de $\mathbb{F}_{3,10,3}$, se encuentra \[\left[\begin{array}{cc|c} 1 & 0 & 1.01\\ 0 & 1 & 0.938\end{array}\right]\rightarrow \begin{array}{l} x_{_1}=0.101\times 10,\\ x_{_2}=0.938.\end{array}\] Note que la solución bajo $\mathbb{F}_{3,10,3}$ difiere de la solución exacta $x_{_1}=x_{_2}=1$. El error se produce por los redondeos y por la modificación hecha para poder completar el procedimiento de Gauss - Jordan.

En general los errores por redondeo son muy importantes cuando se hace sumas entre números pequeños con números grandes o cuando se divide por números muy pequeños. Por ejemplo considere $x=0.100\times 10^4$ y $y=0.400\times 10$, estos números perteneces al conjunto $\mathbb{F}_{3,10,4}$, y observe que $fl(x)=x$ y $fl(y)=y$. Por lo tanto cuando se hace la suma se obtiene: \[fl(x+y)=fl(1000+4)=fl(0.1004\times 10^3)=1\times 10^3=x.\] Es decir, la información de $y$ se pierde completamente.

Hay tres estrategias conocidas como el pivoteo parcial, pivoteo parcial escalado y pivoteo completo que buscan evitar la división por números grandes, para así reducir un poco los errores por redondeo. Considere una matriz $A$ de tamaño $m\times n$ y veamos en que consisten estos dos métodos:

  • Pivoteo Parcial: En cada paso $k$ de la eliminación Gaussiana se elige en calidad de pivote a la entrada con índice $(k,p)$ tal que \begin{equation} |A_{_{kp}}|=\max_{p\leq i \leq m}|A_{_{ip}}|, \end{equation} es decir, la mayor entrada de la $p$ - ésima columna empezando con $(p,p)$ hasta $(m,p)$. Por ejemplo, supongamos que vamos a realizar el paso $p$ de la eliminación Gaussiana, entonces nuestra matriz tendrá una forma parecida a esto: \[\left[\begin{array}{cccccc|c} * & * & * & & * & & * \\ 0 & * & * & & * & & * \\ 0 & 0 & * & & * & & * \\ \vdots & \vdots & \vdots & & \vdots & & \vdots \\ 0 & 0 & 0 & & A_{_{pp}} & & * \\ \vdots & \vdots & \vdots & & \vdots & & \vdots \\ 0 & 0 & 0 & & A_{_{mp}} & & * \\ \end{array}\right]_{m\times n}\] Y suponga que nuestro máximo es $|A_{_{kp}}|=\max_{p\leq i \leq m}|A_{_{ip}}|$, luego la entrada $A_{kp}$ se usará como pivote. Esto significa que se aplicara el intercambio de la filas $R_{_p}\leftrightarrow R_{_k}$ y luego se aplican las operaciones elementales para eliminar las entradas desde $(p+1,p)$ hasta $(m,p)$.
  • Pivoteo parcial escalado: Suponga que estamos en el paso $p$. En cada renglón $i$, con $p\leq i \leq m$, se calcula el valor máximo absoluto en la parte principal de la matriz: \[s_{_i}=\max_{p\leq j\leq m} |A_{_{ij}}|.\] Suponga que $s_{_i}>0$ para todo $i\in \{k,\dots, m\}$. Si se elige el menor entero $q$ con \[\frac{|A_{_{qk}}|}{s_{_q}}=\max_{p\leq i\leq m}\frac{|A_{_{ip}}|}{s_{_i}}.\] En otras palabras, sea $q$ el menor de los índices $i$ en los cuales la expresión $\frac{|A_{_{ip}}|}{s_{_i}}$ alcanza el máximo. Si $q\neq p$, entonces se intercambian los renglones $p$ y $q$ y luego se eliminan las entradas por debajo de $(p,p)$.
  • Pivoteo completo: En el $k$ - ésimo paso de la eliminación Gaussina se buscan los índices $p,q\in \{k,\dots, m\}$ tales que \[|A_{_{rs}}|=\max_{\substack{p\leq i\leq m \\ p\leq j\leq n}}|A_{_{ij}}|.\] En otras palabras, se busca el máximo entre los números $|A_{_{ij}}|$ con $p\leq i \leq m$ y $p\leq j \leq m$. En este caso, antes de efectuar al paso $p$ de la eliminación Gaussiana, la matriz tendrá una forma como esta: \[\left[\begin{array}{ccccccc|c} * & * & * & & * & & * & *\\ 0 & * & * & & * & & * & * \\ 0 & 0 & * & & * & & * & *\\ \vdots & \vdots & \vdots & & \vdots & & \vdots & \vdots\\ 0 & 0 & 0 & & A_{_{pp}} & \cdots & A_{_{pn}} & * \\ \vdots & \vdots & \vdots & & \vdots & & \vdots & *\\ 0 & 0 & 0 & & A_{_{mp}} & \cdots & A_{_{mn}} & *\\ \end{array}\right]_{m\times n}\] Si $(r,s)\neq (p,p)$ entonces se intercambian los renglones y las columnas de tal manera que la entrada $A_{_{rs}}$ se pone en la posición $(p,p)$ y luego se aplican las operaciones elementales para anular las entradas por debajo de $(p,p)$.
Esperamos que esta entrada haya sido de su agrado, nos vemos en otra ocasión.

Referencias