Visualizar las variables activas
--> who
Borrar una cierta variable x
--> clear x
miércoles, 2 de agosto de 2023
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
$$\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.
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.
Funciones en SciLab
Podemos editar una función en el propio panel de comandos, sin necesidad de almacenarlo en un fichero de texto .sci y cargarlo a continuación. Por ejemplo:
-->function f=cuadxy(x,y)
-->f=x**2+y**2;
-->endfunction
Para hacer actuar la función que acabamos de definir, para calcular, pongamos que 2^3, escribiremos:
-->cuadxy(2,3)
ans =
13.
// Comprobación numérica de
// $\lim_{x \rightarrow 0}\,\dfrac{\sin x}{x}=1$
//
//
// Cargar el fichero ( al que se refiere esta imagen )
// en SciLab desde la consola del programa:
// exec('C:\...\mifuncion.sci')
// y, a continuación, para ejecutar la función invocar
// el nombre de la función con los argumentos de entrada
// y de salida
// Ejemplo con un paso incremental de 0.001
// teclear a continuación de "-->":
// [vi,vd]=limites_i_d(0.001)
//
// Joan Aranes Clua
// 2015
// --------------------------------------------------------------------
// valor avanzado de $\lim_{x \rightarrow 0^{+}}\,\dfrac{\sin x}{x}$
function [vi,vd]=limites_i_d(incremento)
disp("límites por la derecha ( vd ) y por la izquierda (vi)")
vd=0;
incremento_neg=-1*incremento
for i=1:incremento_neg:0.0000001
vd=sin(i)/i;
//disp(vd);
end
//
// valor avanzado de $\lim_{x \rightarrow 0^{-}}\,\dfrac{\sin x}{x}$
//vi=0
for j=-1:incremento:-0.0000001
vi=sin(j)/j;
//disp(vi);
end
endfunction
$\diamond$
-->function f=cuadxy(x,y)
-->f=x**2+y**2;
-->endfunction
Para hacer actuar la función que acabamos de definir, para calcular, pongamos que 2^3, escribiremos:
-->cuadxy(2,3)
ans =
13.
Ejemplo
// --------------------------------------------------------------------// Comprobación numérica de
// $\lim_{x \rightarrow 0}\,\dfrac{\sin x}{x}=1$
//
//
// Cargar el fichero ( al que se refiere esta imagen )
// en SciLab desde la consola del programa:
// exec('C:\...\mifuncion.sci')
// y, a continuación, para ejecutar la función invocar
// el nombre de la función con los argumentos de entrada
// y de salida
// Ejemplo con un paso incremental de 0.001
// teclear a continuación de "-->":
// [vi,vd]=limites_i_d(0.001)
//
// Joan Aranes Clua
// 2015
// --------------------------------------------------------------------
// valor avanzado de $\lim_{x \rightarrow 0^{+}}\,\dfrac{\sin x}{x}$
function [vi,vd]=limites_i_d(incremento)
disp("límites por la derecha ( vd ) y por la izquierda (vi)")
vd=0;
incremento_neg=-1*incremento
for i=1:incremento_neg:0.0000001
vd=sin(i)/i;
//disp(vd);
end
//
// valor avanzado de $\lim_{x \rightarrow 0^{-}}\,\dfrac{\sin x}{x}$
//vi=0
for j=-1:incremento:-0.0000001
vi=sin(j)/j;
//disp(vi);
end
endfunction
$\diamond$
martes, 1 de agosto de 2023
Productos en espacios vectoriales
Se distinguen tres tipos de productos en espacios vectoriales de dimensión finita:
- producto escalar
- producto vectorial
- producto tensorial
El producto escalar de dos vectores produce un escalar, el producto vectorial produce un vector y el producto tensorial produce una aplicación lineal. Cuando trabajamos en una base, el producto escalar (es un escalar) se puede calcular como $x^t\,y$ (entendiendo que $x$ e $y$ son vectores columna de $n$ componentes) y el producto tensorial como $x\,y^t$, que es una matriz $n \times n$. El producto vectorial debido a su naturaleza (es un vector) alternante no puede expresarse como un producto matricial.$\diamond$
Matriz de permutación en una factorización LU
Una matriz de permutación $P$ cumple que $PA=LU$ y, además, es simétrica, esto es, $P^tP=I$. Por consiguiente, se cumple que $A=P^t\,L\,U$. $\diamond$
Factorización QR de una matriz
Sea la matriz invertible $A=\begin{pmatrix}2&-1&0\\0&0&-2\\0&2&-1\end{pmatrix}$. Los vectores columna de esta matriz son $a^{1}=\begin{pmatrix}2&0&0\end{pmatrix}$, $a^{2}=\begin{pmatrix}-1&0&2\end{pmatrix}$, $a^{3}=\begin{pmatrix}0&-2&-1\end{pmatrix}$. Los vectores que constituyen una base ortogonal son $p^{i}\;,\,i=1,2,3$ podemos obtenerlos por el método de Gram-Schmidt
$u^{1}=a^1=\begin{pmatrix}2\\0\\0\end{pmatrix}$
$u^{2}=a^2-\langle a^2,u^1 \rangle \,u^1=\begin{pmatrix}0\\0\\2\end{pmatrix}$
$u^{3}=a^3-\langle a^3,u^1\rangle\,u^1 - \langle u^3,u^2\rangle \,u^2=a^{3}-0\,p^{1}+\dfrac{1}{2}\,p^{2}=\begin{pmatrix}0\\-2\\0\end{pmatrix}$
luego la matriz de paso a la nueva base es la siguiente matriz ortogonal $P=\begin{pmatrix}2&0&0\\0&0&-2\\0&2&0\end{pmatrix}$
con lo que una matriz ortonormal se obtiene normalizando los vectores columnas de la anterior (los vectores de la base ortogonal): $q_1=\dfrac{u^1}{\left\|u^1\right\|}$, $q_2=\dfrac{u^2}{\left\|u^2\right\|}$, $q_3=\dfrac{u^3}{\left\|u^3\right\|}$ por tanto
$Q=\begin{pmatrix}1&0&0\\0&0&-1\\0&1&0\end{pmatrix}$, como puede comprobarse por la definición de matriz ortogonal: $Q^t=Q^{-1}$, esto es $QQ^t=I$
Queremos expresar la matriz $A$ de la forma $A=QR$, siendo $Q$ una matriz ortogonal (que ya hemos obtenido a partir de las coordenadas de los vectores de la nueva base $\{p^i\}\,,i=1,2,3$) dispuestos en columnas y normalizados) y $R$ una matriz triangular superior.
Ya tenemos la matriz ortogonal $Q$. Veamos ahora cómo podemos determinar la matriz triangular superior $R$:
De las ecuaciones escritas arriba podemos escribir, en forma matricial, que $A=P-PM$, y por tanto, $A=P(I+M)$, donde $I$ es la matriz identidad y $M$ puede verse fácilmente que es $\begin{pmatrix}0&m_{12}&m_{13}\\0&0&m_{23}\\0&0&0\end{pmatrix}$, siendo sus coeficientes de la forma $m_{ij}=\dfrac{a^j\,p^i}{\left\|p^i\right\|^2}\,\text{para}\,j\gt i$
Entonces, como al normalizar los vectores columna se han dividido sus coordenadas, $p^{i}\,,i=1,2,3$, por la norma respectiva de cada vector, $\left\|p^i\right\|$, tenemos que $A=P(I+M)=QD(I+M)=Q\,(D(I+M))$, y como esto debe ser igual a $QR$, se tiene que $R=D(I+M)$, siendo $D$ una matriz diagonal que, como acabamos de justificar, resulta ser $D=(d_{ii})$, con $d_{ii}=\left\|p^i\right\|$; y, en nuestro caso, resulta ser $D=\begin{pmatrix}2&0&0\\0&2&0\\0&0&2\end{pmatrix}$, por lo que la matriz triangular superior es $R=\begin{pmatrix}2&-1&0\\0&2&-1\\0&0&2\end{pmatrix}$. Así, una descomposción $QR$ de $A$ puede comprobarse que es $A=\begin{pmatrix}1&0&0\\0&0&-1\\0&1&0\end{pmatrix}\begin{pmatrix}2&-1&0\\0&2&-1\\0&0&2\end{pmatrix}$. $\square$
$u^{1}=a^1=\begin{pmatrix}2\\0\\0\end{pmatrix}$
$u^{2}=a^2-\langle a^2,u^1 \rangle \,u^1=\begin{pmatrix}0\\0\\2\end{pmatrix}$
$u^{3}=a^3-\langle a^3,u^1\rangle\,u^1 - \langle u^3,u^2\rangle \,u^2=a^{3}-0\,p^{1}+\dfrac{1}{2}\,p^{2}=\begin{pmatrix}0\\-2\\0\end{pmatrix}$
luego la matriz de paso a la nueva base es la siguiente matriz ortogonal $P=\begin{pmatrix}2&0&0\\0&0&-2\\0&2&0\end{pmatrix}$
con lo que una matriz ortonormal se obtiene normalizando los vectores columnas de la anterior (los vectores de la base ortogonal): $q_1=\dfrac{u^1}{\left\|u^1\right\|}$, $q_2=\dfrac{u^2}{\left\|u^2\right\|}$, $q_3=\dfrac{u^3}{\left\|u^3\right\|}$ por tanto
$Q=\begin{pmatrix}1&0&0\\0&0&-1\\0&1&0\end{pmatrix}$, como puede comprobarse por la definición de matriz ortogonal: $Q^t=Q^{-1}$, esto es $QQ^t=I$
Queremos expresar la matriz $A$ de la forma $A=QR$, siendo $Q$ una matriz ortogonal (que ya hemos obtenido a partir de las coordenadas de los vectores de la nueva base $\{p^i\}\,,i=1,2,3$) dispuestos en columnas y normalizados) y $R$ una matriz triangular superior.
Ya tenemos la matriz ortogonal $Q$. Veamos ahora cómo podemos determinar la matriz triangular superior $R$:
De las ecuaciones escritas arriba podemos escribir, en forma matricial, que $A=P-PM$, y por tanto, $A=P(I+M)$, donde $I$ es la matriz identidad y $M$ puede verse fácilmente que es $\begin{pmatrix}0&m_{12}&m_{13}\\0&0&m_{23}\\0&0&0\end{pmatrix}$, siendo sus coeficientes de la forma $m_{ij}=\dfrac{a^j\,p^i}{\left\|p^i\right\|^2}\,\text{para}\,j\gt i$
Entonces, como al normalizar los vectores columna se han dividido sus coordenadas, $p^{i}\,,i=1,2,3$, por la norma respectiva de cada vector, $\left\|p^i\right\|$, tenemos que $A=P(I+M)=QD(I+M)=Q\,(D(I+M))$, y como esto debe ser igual a $QR$, se tiene que $R=D(I+M)$, siendo $D$ una matriz diagonal que, como acabamos de justificar, resulta ser $D=(d_{ii})$, con $d_{ii}=\left\|p^i\right\|$; y, en nuestro caso, resulta ser $D=\begin{pmatrix}2&0&0\\0&2&0\\0&0&2\end{pmatrix}$, por lo que la matriz triangular superior es $R=\begin{pmatrix}2&-1&0\\0&2&-1\\0&0&2\end{pmatrix}$. Así, una descomposción $QR$ de $A$ puede comprobarse que es $A=\begin{pmatrix}1&0&0\\0&0&-1\\0&1&0\end{pmatrix}\begin{pmatrix}2&-1&0\\0&2&-1\\0&0&2\end{pmatrix}$. $\square$
Suscribirse a:
Entradas (Atom)