Mostrando entradas con la etiqueta análisis numérico. Mostrar todas las entradas
Mostrando entradas con la etiqueta análisis numérico. Mostrar todas las entradas

jueves, 3 de agosto de 2023

Breve nota sobre la interpolación de funciones

Para encontrar el polinomio interpolador de grado $n$, dados $n+1$ puntos, podemos optar por plantear un sistema de $n+1$ ecuaciones con $n+1$ incógnitas. Cuando $n$ es pequeño, este método ( de los coeficientes indeterminados ) es viable, sin embargo es claro que no lo es cuando $n$ ( es grande ). Para encontrar el polinomio interpolador podemos utilizar, entre otros, el método de los polinomios interpoladores de Lagrange.

Aritmética de punto fijo versus aritmética de punto flotante

0. ARITMÉTICA EXACTA
Idealmente, para obtener resultados exactos en las operaciones aritméticas, no debemos prescindir de ninguna cifra decimal al llevar a cabo las sucesivas operaciones; sin embargo, eso no es viable en la mayor parte de los casos, pues la capacidad de representación de las cantidades no es infinita. Por ello, conviene hablar de de dos tipos de aritméticas inexactas ( asumiendo que vamos a redondear los resultados parciales ).

1. ARITMÉTICAS INEXACTAS

1.1 ARITMÉTICA DE PUNTO FIJO
Si en los cálculos fijamos el número de cifras decimales o bien el número de fijas significativas en los cálculos, decimos que operamos en aritmética de punto fijo.

Ejemplo:
Consideremos la operación $9641,45 \cdot 0,012$ -- que en aritmética exacta da $115.6974$ -- que nos proponemos realizar conservando dos cifras decimales en cada factor ( y por tanto, también en el resultado ), entonces debemos hacer: $9641,45 \cdot 0,01$ que es igual a $96.41$ ( con dos cifras decimales ). Como podemos observar, la aritmética de punto fijo puede dar resultados muy imprecisos.

1.2 ARITMÉTICA DE PUNTO FLOTANTE
Esta forma de representar los números viene a rebajar el problema de la imprecisión apuntada en el ejemplo anterior ( al utilizar aritmética de punto fijo ). Si representamos las cantidades de modo que se conserve un número convenido de cifras significativa ( en este caso, son seis, ya es ese el mayor número de c.s. de entre los dos factores ) escribiéndolas en el formato $0,m \times 10^e$, dejando así flotando el punto ( la coma ) decimal:
$9641,45=0,964146 \times 10^4$ y $0,012 = 0,120000 \times 10^{-1}$, al realizar el producto de ambas cantidades tenemos ahora $0,964146 \cdot 0,120000 \cdot 10^{4-1}=0,115697 \cdot 10^2$, es decir, $115,697$, resultado que es mucho más preciso que el obtenido utilizando aritmética de punto fijo.


miércoles, 2 de agosto de 2023

Acerca de la precisión de la máquina (unidad de redondeo)

Se suele utilizar la siguiente notación para realizar las operaciones de truncamiento y redondeo ( si bien, recordemos que lo usual será hacerlas por redondeo ): $$\mathbb{R}\ni \text{float}(x)=\pm\,0.d_{1}d_{1}\ldots \times b^e \approx \left\{\begin{matrix}\mathcal{T}(x;n)=\text{Ent}(b^n \times 0.m)\times b^{e-n} \\ \mathcal{R}(x;n)=\text{Ent}(b^n \times 0.m+\mu)\times b^{e-n}\end{matrix}\right.$$ donde $\mu$ denota media unidad en el sistema de numeracíon de base $b$ que se emplee; así, si $b=10$ es claro que $\mu=0.5$, y, si, por ejemplo, $b=2$, entonces $\mu=0.1_{2}$ ya que $0.1_{2}=2^{-1}=0.5_{10}$


Nota: La función parte entera $\text{Ent}(.)$ suele abreviarse de la forma $[.]$

En el caso de que $x$ sea un número negativo, se define $\mathcal{R}(x)\overset{\text{def}}{=}-\mathcal{R}(-x)$ y se cumple que $|x-\mathcal{R}(x)|\le \dfrac{1}{2}\,b^{e-n}$.

La cantidad $\dfrac{1}{2}\,b^{e-n}$ se conoce como unidad de redondeo o también como precisión de la máquina, y representa una cota razonable del error absoluto.

Así, por ejemplo, es evidente que al aproximar por redondeo el número $\pi$ con tres cifras significativas obtenemos $\mathcal{R}(\pi)=3.14$; de acuerdo con lo dicho arriba, convenimos así que $\mathcal{R}(\pi)=0.314\times 10^1$, con lo cual $e=1$ ( exponente de la potencia de base $10$ ) y $n=3$ ( precisisón de la mantisa - número de dígitos en la mantisa ( con el dígito entero igual a cero ) -, que es $n=3$ ), por lo que el máximo error de redondeo en esta aproximación $|x-\mathcal{R}(x)|\prec \dfrac{1}{2}\,b^{e-n}$ ( precisión de la máquina, o unidad de redondeo, que sabemos que es igual a media unidad del orden de la última cifra considerada ), y, en efecto, como $3\prec \pi \prec 4$, $e:=1$ ( exponente de la potencia de base $10$ ) y al aproximar con $3$ cifras significativas, $n:=3$, luego dicha precisión de la máquina ( o cota razonable de error absoluto ) en dicha aproximación es $\dfrac{1}{2}\,b^{1-3}=\dfrac{1}{2}\,10^{-2}=\dfrac{1}{2}\,10^{-2}=0.005$
$\square$

Acerca de la representación IEEE Standard 754 en punto flotante

Actualmente, la representación que usan la mayoría de los computadores es la llamada IEEE Standard 754 en punto flotante. Está basada en los tres elementos que se ha mencionado antes: signo, mantisa y exponente. Si la base es binaria, $b=2$, el primer dígito (dígito principal) es, siempre, $1$, esto es, $d_1=1$, por lo que no haría falta almacenarlo en memoria (en cuyo caso nos referiríamos a él como dígito principal implícito).

En precisión simple, el estándard utiliza $32$ bits ($4$ bytes), de los cuales $8$ bits ($1$ byte) es para el exponente $E$ (representando números enteros no negativos que varían entre $0$ y $255$, esto es, los que se pueden codificar en base $2$ con $8$ bits); $23$ bits son para la parte fraccionaria, $F$, y $1$ bit para el signo, $S$. Cada bit almacena un $0$ o bien un $1$. Así, el valor asignado a una representación es $$(-1)^{S}\times 2^{E-127} \times \mathbb{1}.F$$

Se utiliza un sesgo de $127$ aplicado al valor del exponente almacenado $E$, dando lugar al exponente verdadero $E-127$, lo cual tiene como finalidad el poder representar tanto exponentes positivos como negativos; de este modo, un valor almacenado $E$ representa el valor verdadero $E-127$.

Como ya hemos visto que $E$ varía entre $0$ y $255$, entonces $E-127$ toma valores en el conjunto de $2^8=256$ elementos: $\{-127,-126,\ldots,0,1,2,\ldots,128\}$, y , por tanto, $2^{E-127}$ varía entre $2^{-127}\approx 10^{-38}$ y $2^{128}\approx 10^{38}$.

El ordinal de cada término de la sucesión $-127,-126,\ldots,0,1,2,\ldots,128$ (empezando a contar desde $0$) se denomina característica. El valor de la característica varía pues de $0$ (para el valor del primer término $-127$) a $255$ (para el valor del último término $128$), y representa el valor del exponente almacenado $E$.

Designando por $n$ al número de bits del campo del exponente (en nuestro caso $n:=8$), esta característica es igual, por tanto, al valor del exponente verdadero $E-127$ más lo que viene a denominarse desplazamiento y que es fácil ver que es igual a $2^{n-1}-1$. Así por ejemplo, si el exponente verdadero es $E-127:=100$, entonces el valor de la característica es $100+2^{8-1}-1=100+128-1=227$; y (otro ejemplo), como seria de esperar si $E-127:=-127$, entonces sus característica es igual a $-127+2^{8-1}-1=-127+128-1=0$; y si $E-127:=-128$, entonces su característica es igual a $128+2^{8-1}-1=128+128-1=255$, que es el último término de la sucesión que hemos comentado.



Por otra parte, El bit del signo es $0$ para positivos y $1$ para negativos.

Ejemplo 1.
La representación $$\underset{\overbrace{1\; \text{bit para el signo}}}{1}\quad \underset{\overbrace{8\; \text{bits para el exponente}}}{01010011} \quad\underset{\overbrace{23\; \text{bits para la mantisa}}} {10011110\; 00101000\; 0101000}$$ es tal que:

i) El signo es negativo por ser el valor del bit $S$ igual a $1$, y, por tanto, $(-1)^S=(-1)^1=-1$

ii) la mantisa $1.F$ es

1+$(1\cdot 2^{-1}+0\cdot 2^{-2}+0\cdot 2^{-3}+1\cdot 2^{-4}+\overset{\underbrace{23}}{\ldots}+0\cdot 2^{-22}+0\cdot 2^{-23})$
donde a partir del segundo sumando se expresa la cantidad $F$, siendo el primer sumando a la izquierda, el primer uno, que viene de la notación $1.F$, luego la mantisa corresponde al número $$\dfrac{2097151}{1048576}$$

iii) El exponente verdadero de la potencia $2^{E-127}$ es igual a 0$1010011_{2}$, o lo que es lo mismo 0b$1010011$ — en base $2$ los números suelen empezar por 0 o bien por 0b— corresponde a $83_{10}-127=-44$
Por tanto, la representación estándard indicada da el valor

$$-( \dfrac{2097151}{1048576} ) \times 2^{83-127}$$ es decir $$-( \dfrac{2097151}{1048576} ) \times 2^{-44}\approx -0.1136868\times 10^{-12}$$ pues con $23$ bits de mantisa el número de dígitos/cifras significativas correctas que podemos asegurar es de $7$


-oOo-



Nota: En doble precisión, se reservan $11$ bits es para representar el exponente $E$, donde el sesgo es ahora de $2^{11-1}-1=2^{10}-1=1023$.
-oOo-
Referencias:
  [1] Moreno, C.: Introducción al Cálculo Numérico, UNED, Madrid, 2011
  [2] García Valle, J.L.: Matemáticas especiales para computación, McGraw-Hill, Madrid, 1988

Uso de arrays con NumPy y SciPy ( Python)

Python: lista de reproducción del curso de Python de la Universitat d'Alacant
-oOo-

NumPy
-oOo-

ScyPy ( ScyPy es un conjunto de paquetes específicos que se basan en los de NumPy )
-oOo-

SymPy ( ScyPy es sistema CAS, que viene como un conjunto de paquetes específicos de cálculo simbólico basados en Python )
-oOo-

Más vídeos similares, aquí

Acerca de los errores de truncamiento empleando representaciones máquina: inestabilidades

Si en un entorno que usa el estándar IEEE Double Precision, se realiza el cálculo $e^{-10}$ mediante el truncamiento de la serie de potencias $$e^x=\displaystyle \sum_{i=0}^{30}\,\dfrac{x^i}{i!} \approx 1+\dfrac{x}{1!}+\dfrac{x^2}{3!}+\dfrac{x^3}{3!}+\ldots+\dfrac{x^{30}}{30!}$$, y el resultado aproximado es $9.703415796025505 \times 10^{-4}$. El error de truncamiento viene dado por $\dfrac{e^\xi}{(n+1)!}$ donde $-10 \prec \xi \prec 0$, y $n=30$, por lo que tomamos el error de dicha aproximación como $E(-10):=\dfrac{1}{31!} \overset{n:=30}{\succ} \dfrac{e^{\xi}}{31!}$.

Ahora bien, el resultado exacto de $e^{-10}$ es igual a $4.539992976248485 \times 10^{-5}$, el error relativo de dicho resultado aproximado es enorme (!): $$\dfrac{|4.539992976248485 \times 10^{-5}-9.703415796025505 \times 10^{-4}|}{4.539992976248485 \times 10^{-5}} \approx 20 = 2\,000\,\% $$ Tal inexactitud se debe a la acumulación de los pequeños errores causados al operar en aritmética finita, los cuales se han magnificado a lo largo del cálculo, por lo que podemos decir que el algoritmo de cálculo que hemos utilizado es muy inestable.

Nota: Tengamos en cuenta que si se calcula $e^{-10}$ de la forma $\dfrac{1}{e^{10}}$, desarrollando en serie de potencias el denominador, se obtiene como resultado $4,539993338712231 \times 10^{-5}$ que es muy próximo al valor exacto $4.539992976248485 \times 10^{-5}$, obteniéndose en ese caso un error relativo $\dfrac{|4.539992976248485 \times 10^{-5}-4,539993338712231 \times 10^{-5} \times 10^{-4}|}{4.539992976248485 \times 10^{-5}}\approx$ $\approx 8\times 10^{-8}=8\times 10^{-6}\,\%$

En conclusión, la causa de que un procedimiento de cálculo produzca resultados tan imprecisos es debida a la inestabilidad del algoritmo de cálculo ( un error en alguna etapa del cálculo se propaga de manera creciente a las siguientes etapas ). Así que, para evitar tal efecto - los errores debidos a los cálculos en artimética finita de un sistema de cálculo automático son inevitable - es necesario emplear algoritmos estables. $\square$

Categorías de los errores en el análisis numérico

CATEGORÍAS DE ERRORES EN LOS CÁLCULOS NUMÉRICOS:

1. Errores humanos (elección del método, mala interpretación del mismo, errores en los algoritmos, etcétera)

2. Errores en los datos -> noción de condicionamiento (si un pequeño error en los datos llevar a errores grandes en los resultados, decimos que el problema está mal condicionado)

3. Errores de truncamiento -> noción de convergencia de los métodos numéricos (se deben al truncamiento de series infinitas, por ejemplo)

4. Errores de redondeo -> noción de estabilidad de los métodos numéricos (se deben a la acumulación de los errores de redondeo) Notas:[1]

La inexactitud y las representaciones máquina finitas

  ¿Qué entendemos por aritmética exacta? Idealmente, para obtener resultados exactos en las operaciones aritméticas, no debemos prescindir de ninguna cifra decimal al llevar a cabo las sucesivas operaciones; sin embargo, eso no es viable en la mayor parte de los casos, pues la capacidad de representación de las cantidades no es infinita. Por ello, conviene hablar de de dos tipos de aritméticas inexactas (asumiendo que vamos a redondear los datos y, por supuesto, los resultados parciales). Ahora bien, la realización de esa aritmética inexacta a la que inexorablemente nos vemos abocamos en la mayor parte de las situaciones de cálculo podemos hacerla de varias manera, según el formato de representación de las cantidades (en punto (decimal) fijo o punto (decimal) flotante). Una y otra presentan ventajas y desventajas, como puede verse en la siguente referencia:
  https://es.wikipedia.org/wiki/Coma_flotante, Wikipedia

Ejemplo de desbordamiento en una representación finita en coma flotante

ENUNCIADO. En un sistema binario ( $b=2$ ), con un número fijo de dígitos, $n=3$, el conjunto de números ( que no son negativos ) que pueden ser representados de este modo, son
Representación: 000 | 001 | 010 | 011 | 100 | 101 | 110 | 111 
Valor real:       0 |   1 |   2 |   3 |   4 |   5 |   6 |   7
Un aspecto relevante de este sistema de representación es que puede producirse un desbordamiento en las operaciones aritméticas debido a que el conjunto es acotado. Por ejemplo, la suma de 3 y 5 produce un número que no puede ser representado en este sistema. $\square$

Aclaración:

$0=0\cdot 2^{0}+0\cdot 2^{1}+0\cdot 2^{2}$
$1=1\cdot 2^{0}+0\cdot 2^{1}+0\cdot 2^{2}$
$2=0\cdot 2^{0}+1\cdot 2^{1}+0\cdot 2^{2}$
$3=1\cdot 1^{0}+1\cdot 2^{1}+0\cdot 2^{2}$
$4=0\cdot 2^{0}+0\cdot 2^{1}+1\cdot 2^{2}$
etcétera

Ejemplo de rango de representación en una aritmética finita de coma/punto flotante

En un sistema binario, $b=2$, si el número de dígitos es $n=2$, el conjunto de números (que no son negativos) que pueden ser representados en punto flotante con el rango $-2\le e\le 2$ es el siguiente:


Así, se obtienen las siguientes cantidades representadas:
$$\mathcal{M}=\{0,\dfrac{1}{16},\dfrac{1}{8},\dfrac{3}{16},\dfrac{1}{4},\dfrac{3}{8},\dfrac{1}{2},\dfrac{3}{4},\dfrac{1}{16},1,\dfrac{3}{2},2,3\}$$

En efecto ( aclaración ):
Para $e=-2$ se obtiene $0.00 \times 2^{-2}=0$
Para $e=-1$ se obtiene $0.00 \times 2^{-1}=0$
Para $e=0$ se obtiene $0.00 \times 2^0=0$
Para $e=1$ se obtiene $0.00 \times 2^1=0$
Para $e=-2$se obtiene $0.00 \times 2^2=0$
---
Para $e=-2$ se obtiene $\dfrac{1}{4} \times 2^{-2}=\dfrac{1}{16}$
Para $e=-1$ se obtiene $\dfrac{1}{4} \times 2^{-1}=\dfrac{1}{8}$
Para $e=0$ se obtiene $\dfrac{1}{4} \times 2^0=\dfrac{1}{4}$
Para $e=1$ se obtiene $\dfrac{1}{4} \times 2^1=\dfrac{1}{2}$
Para $e=-2$ se obtiene $\dfrac{1}{4} \times 2^2=1$
---
Para $e=-2$ se obtiene $\dfrac{1}{2} \times 2^{-2}=\dfrac{1}{8}$
Para $e=-1$ se obtiene $\dfrac{1}{2} \times 2^{-1}=\dfrac{1}{4}$
Para $e=0$ se obtiene $\dfrac{1}{2} \times 2^0=\dfrac{1}{2}$
Para $e=1$ se obtiene $\dfrac{1}{2} \times 2^1=1$
Para $e=-2$ se obtiene $\dfrac{1}{2} \times 2^2=2$
---
Para $e=-2$ se obtiene $\dfrac{3}{4} \times 2^{-2}=\dfrac{3}{16}$
Para $e=-1$ se obtiene $\dfrac{3}{4} \times 2^{-1}=\dfrac{3}{8}$
Para $e=0$ se obtiene $\dfrac{3}{4} \times 2^0=\dfrac{3}{4}$
Para $e=1$ se obtiene $\dfrac{3}{4} \times 2^1=\dfrac{3}{2}$
Para $e=-2$ se obtiene $\dfrac{3}{4} \times 2^2=3$

-oOo-
Comentarios:
En esta representación es de hacer notar que:
1. Los números representados no son equidistantes y se observa mayor densidad cerca del $0$
2. Un mismo número representado puede serlo de varias maneras distintas
3. Como el conjunto es acotado, se puede producir un desbordamiento en las operaciones aritméticas; por ejemplo, la suma de $1$ y $3$, esto es $4$, cae fuera del conjunto de valores que pueden representarse; por ejemplo, $1+3=4 \succ \text{máx}(\mathcal{M})=3$
4. El resultado obtenido de algunas operaciones aritméticas da valores que, aunque no se produzca desbordamiento ( por obtener un número mayor que el máximo de dicha representación ), no se obtiene un número del conjunto; por ejemplo $2+\dfrac{1}{8}=\dfrac{17}{8}$, que si bien es mayor que $0=\text{mín}(\mathcal{M})$ y menor que $\text{máx}(\mathcal{M})=3$, no corresponde a ninguno valor del conjunto de representación, $\mathcal{M}$

-oOo-

2'. Para que un mismo número tenga representaciones distitnas en un mismo sistema, establecemos que el primer dígito de la representación, $d_1$, sea distinto de $0$, por lo que diremos que, siendo así, un sistma de prepresentación en punto flotante está normalizada.

Actualmente, la representación que usan la mayoría de los computadores es la llamada IEEE Standard 754 en punto flotante. Está basada en los tres elementos que se ha mencionado antes: signo, mantisa y exponente. Si la base es binaria, $b=2$, el primer dígito ( dígito principal ) es, siempre, $1$, esto es, $d_1=1$, por lo que no haría falta almacenarlo en memoria ( en cuyo caso nos referiríamos a él como dígito principal implícito ).

En precisión simple, el estándard utiliza $32$ bits ( $4$ bytes ), de los cuales $8$ bits ( $1$ byte ) es para el exponente $E$; $23$ bits son para la parte fraccionaria, $F$, y $1$ bit para el signo, $S$. Cada bit almacena un $0$ o bien un $1$. Así, el valor asignado a una representación es $$(-1)^{S}\times 2^{E-127} \times \mathbb{1}.F$$ donde el sesgo de $-127$ añadido al valor positivo almacenado tiene como finalidad el poder representar tanto exponentes positivos como negativos; de este modo, un valor almacenado $E$ representa un valor real $E-127$. El bit del signo es $0$ para positivos y $1$ para negativos.

Por ejemplo:
La representación $\underset{\overbrace{\text{signo}}}{1}$     $\underset{\overbrace{\text{exponente}}}{01010011}$     $\underset{\overbrace{\text{mantisa}}} {10011110\; 00101000\; 0101000}$ es tal que:
i) El signo es $-$ por ser el bit del signo igual a $1$
ii) la mantisa es
$\mathbb{1}+\left(1\cdot 2^{-1}+1\cdot 2^{-2}+0\cdot 2^{-3}+1\cdot 2^{-4}+\ldots+1\cdot 2^{-20}+0\cdot 2^{-21}+0\cdot 2^{-22}+0\cdot 2^{-23}\right)$
—donde entre paréntesis se expresa la cantidad $F$, siendo el primer sumando a la izquierda, el primer uno, es el uno de la notación $1.F$—, que es igual a $\dfrac{2097151}{1048576}$
iii) y el exponente 0$1010011_{2}$, o lo que es lo mismo 0b$1010011$ ( en base $2$ los números suelen empezar por 0 o bien por 0b ), representa al número $83_{10}$
Por tanto, la representación estándard indicada da el valor
$$-( \dfrac{2097151}{1048576} ) \times 2^{83-127}$$ es decir $$-( \dfrac{2097151}{1048576} ) \times 2^{-44}$$

Nota: En doble precisión, se reservan $11$ bits es para el exponente $E$, donde el sesgo es ahora de $1023$.

Ejemplos de comportamientos en representación finita de coma flotante

ENUNCIADO. Sea $x \succ 0$ un número que se representa exacto en una aritmética finita de punto flotante. Analícese el comportamiento del cociente $$\dfrac{\sin\,(x+h)-\sin\,x}{h}$$ cuando se consideran valores de $h$ próximos a $0$ y los cálculos se realizan en una aritmética de punto flotante.

SOLUCIÓN. Para valores de $h$ suficientemente pequeños tenemos que $\mathcal{R}(x+h)=x$, por lo que el numerador de la expresión se anula, habida cuenta de que $\sin(x+h)\approx \sin\,x$ si $h \ll 1$, mientras que el denominador no lo hace. En consecuencia, el valor del cociente se anula. $\square$

Sobre la transformación de Householder

[ artículo de Wikipedia ]

Matrices de Hessemberg

En este breve artículo de Wikipedia se describen las matrices de Hessenberg y se hacen algunos comentarios sobre su utilidad en los algoritmos de cálculo numérico.

Número condición de una matriz

En este artículo de Wikipedia se explica brevemente qué es el número condición de una matriz

Matriz de Fröbenius

En este artículo de Wikipedia se describe la forma de una matriz de Fröbenius ( asociada a la llamada transformación de Gauss ) y se comentan sus propiedades y utilidades en los algoritmos de cálculo numérico

Cuadratura de Gauss

En este artículo de Wikipedia, que es bastante conciso, se explica el concepto de cuadratura de Gauss, método en el que las abscisas no están equiespaciadas

Visualización de variables activas y borrado de las mismas en SciLab

Visualizar las variables activas
--> who

Borrar una cierta variable x
--> clear x

Factorización LU con SciLab

Quereremos resolver el sistema de ecuaciones
$$\left\{\begin{matrix}x+y+z=3\\2x-y+z=2\\x+2y-2z=1\end{matrix}\right.$$

por lo que podemos expresarlo en forma matricial: $AX=b$ siendo $A$ la matriz de los coeficientes del sistema y $b$ la matriz columna de los términos independientes

Vamos a utilizar el método de factorizació de la matriz $A$ como $LU$, de forma que $AX=b \Leftrightarrow LUX=b$, con lo cual ( una vez obtenida la f. LU ) procederemos en dos etapas:
i) $LY=b \Rightarrow Y=L^{-1}\,B$ ( eliminación progresiva )
ii) $UX=Y \Rightarrow X=U^{-1}\,Y$ ( sustitución retrógrada )

Antes que nada, debemos factorizar
Factorización LU
Nota: L es una matriz triangular inferior y U una matriz triangular superior ( E es la matriz de permutación )
-->A
A =

1. 1. 1.
2. - 1. 1.
1. 2. - 2.

-->[L,U,E]=lu(A)
E =

0. 1. 0.
0. 0. 1.
1. 0. 0.
U =

2. - 1. 1.
0. 2.5 - 2.5
0. 0. 2.
L =

1. 0. 0.
0.5 1. 0.
0.5 0.6 1.

Comprobación LU=A ?
-->L*U
ans =

2. - 1. 1.
1. 2. - 2.
1. 1. 1.
Cuidado, se ha utilizado el método de pivotamiento parcial, por lo que el orden de las filas se ha cambiado, de modo que ahora el sistema a resolver es
$$\left\{\begin{matrix}2x-y+z=2\\x+2y-2z=1\\ x+y+z=3\end{matrix}\right.$$ y, por tanto, $$b=\begin{pmatrix}2\\1\\3\end{pmatrix}$$

Procedamos ahora con los dos pasos anunciados ( eliminación progresiva seguida de sustitución retrógrada ):

-->b=[2;1;3]
b =

2.
1.
3.

-->Y=inv(L)*b
Y =

2.
0.
2.

-->X=inv(U)*Y
X =

1.
1.
1.
que es la solución del sistema de ecuaciones


------------------------
Cuidado:
las funciones tril(A) y triu(A) también dan matrices triangulares inferiores y superiores, respectivamente, pero no son las que cumplen A=LU

Factorización QR en SciLab

Factorización QR
Nota: Q es una matriz ortogonal y R una matriz triangular superior

-->A=[1,2;3,4]
A =

1. 2.
3. 4.
-->[Q,R]=qr(A)
R =

- 3.1622777 - 4.4271887
0. - 0.6324555
Q =

- 0.3162278 - 0.9486833
- 0.9486833 0.3162278

-->Q*R
ans =

1. 2.
3. 4.

jueves, 27 de julio de 2023

Submatrices. Submatrices principales y menores principales

Consideremos una matriz $m \times n$ ( $m$ filas y $n$ columnas ). Sin pérdida de generalidad, pongamos que $m:=3$ y $n:=4$ $$A=\begin{pmatrix}a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34}\end{pmatrix}$$ a la cual le corresponde los siguientes conjuntos de índices de filas y columnas, respectivamente: $I=\{1,2,3,4\}$ y $J=\{1,2,3,4\}$. Pues bien, podemos referirnos a una cierta submatriz $B$ de $A$ seleccionando su subconjunto de índices de fila $I_B$ y un subconjunto de índices de columna $J_B$; pongamos que $I_B=\{2,3\}$ y $J_B=\{3,4,\}$, entendiendo dicha submatriz como una matriz de tamaño $(3-2) \times (4-2)$, que, concretamente es $$B\overset{.}{=}A[I_B;J_B]=A[\{2,3\};\{3,4\}]=\begin{pmatrix}a_{23} & a_{24} \\ a_{33} & a_{34} \end{pmatrix}$$ Diremos que una submatriz $B$ de una matriz $A$ es principal si eliminamos las filas y columnas de igual índice, esto es $I_B=J_B$, y lo abreviaremos con la notación $A[I_B]$. Así, por ejemplo, si eliminamos la primera fila y la primera columna, obtenemos la matriz principal $A[\{1\}]=\begin{pmatrix}a_{22} & a_{23} & a_{24} \\ a_{32} & a_{33} & a_{34} \end{pmatrix}$

Llamamos menor complementario de una matriz $A_{m \times n}$ al determinante de alguna de de sus submatrices cuadradas ( que obtenemos al eliminar una o más de de una de sus filas o columnas ).

-oOo-
Matrices cuadradas de orden $m$

Si en una matriz cuadrada de orden $m$ seleccionamos las filas y columnas $I=J=\{1,\ldots,k\}$ ( con $k\le m$ ), obtenemos una matriz principal superior; y si seleccionamos las filas y columnas con $I=J=\{k,\ldots,m\}$ obtenemos una matriz principal inferior.
Así, para una matriz $B=\begin{pmatrix}b_{11} & b_{12} & b_{13} \\ b_{21} & b_{22} & b_{23} \\ b_{31} & b_{32} & b_{33} \end{pmatrix}$, hay tres submatrices principales superiores: $B[\{1\}]=(b_{11})$, $B[\{1,2\}]=\begin{pmatrix}b_{11} & b_{12} \\ b_{21} & b_{22} \end{pmatrix}$ y $B[\{1,2,3\}]=B$.

Y también tres matrices principales inferiores: $B[\{3\}]=(b_{33})$, $B[\{2,3\}]=\begin{pmatrix}b_{22} & b_{23} \\ b_{32} & b_{33} \end{pmatrix}$ y $B[\{1,2,3\}]=B$

Si se elimina una sola fila y una sola columna, los menores complementarios que así se obtienen se denominan primeros menores de dicha matriz cuadrada, y se necesitan para encontrar la matriz de cofactores, la cual es útil para calcular el determinante y la inversa de matrices cuadradas. Además, en el caso de que la submatriz cuadrada formada sea una submatriz principal, los correspondientes primeros menores se denominan menores principales.

Nota. Los determinantes de las matrices principales superiores/inferiores se denominan menores principales superiores/inferiores.

Notación. Al eleminar la fila $i$ y la columna $j$ de una matriz cuadrada $A_{m \times m}$ obtenemos el menor complementario que suele notarse de la forma $M_{ij}$; o, como también suele decirse, el $(i,j)$-ésimo menor complementario de $A$ ) y que, lógicamente, es de orden $m-1$. Dicho menor puede entenderse como el que resulta de eliminar la fila y la columna a las que pertenece el elemento $a_{ij}$ de la matriz $A$. Recordemos que en el caso de obtener un cierto menor eliminando una única fila y una única columna, hablamos de primeros menores, por lo que si se obtienen eliminando dos filas y dos columnas, los llamaremos segundos menores, etcétera.