Matemática numérica I

3 Métodos para sistemas lineares 3 Métodos para sistemas lineares 3.2 Norma e número de condicionamento

Ajude a manter o site livre, gratuito e sem propagandas. Colabore!

3.1 Método da decomposição LU

O método da decomposição LU baseia-se na fatoração da matriz de coeficientes de um sistema linear. Dado um sistema $A\bm{x}=\bm{b}$ , buscamos escrever a matriz $A$ como o produto de uma matriz triangular inferior $L$ (do inglês, lower triangular matrix) por uma matriz triangular superior $U$ (do inglês, upper triangular matrix), isto é,

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}A=LU.

(3.9)

Com isso, o sistema pode ser reescrito na forma

	$\displaystyle A\bm{x}=\bm{b}$		(3.10)
	$\displaystyle(LU)\bm{x}=\bm{b}$		(3.11)
	$\displaystyle L(U\bm{x})=\bm{b}.$		(3.12)

Definindo

\bm{y}:=U\bm{x},

(3.13)

obtemos o sistema triangular

L\bm{y}=\bm{b}.

(3.14)

Depois de resolver esse sistema, obtemos a solução do sistema original resolvendo o sistema triangular

U\bm{x}=\bm{y}.

(3.15)

Ou seja, a decomposição LU nos permite resolver um sistema linear por meio da resolução sucessiva de dois sistemas triangulares.

3.1.1 Sistemas triangulares

Antes de estudarmos como calcular a decomposição LU de uma matriz, vamos discutir a resolução de sistemas triangulares.

Sistema triangular inferior

Um sistema linear triangular inferior tem a forma algébrica

\begin{matrix}a_{1,1}x_{1}&&&&&=b_{1}\\ a_{2,1}x_{1}&+a_{2,2}x_{2}&&&&=b_{2}\\ \vdots&&&&&\vdots\\ a_{n-1,1}x_{1}&+a_{n-1,2}x_{2}&+\cdots&+a_{n-1,n-1}x_{n-1}&&=b_{n-1}\\ a_{n,1}x_{1}&+a_{n,2}x_{2}&+\cdots&+a_{n,n-1}x_{n-1}&+a_{n,n}x_{n}&=b_{n}\end{matrix}

(3.16)

Se $a_{i,i}\neq 0$ para todo $i=1,\ldots,n$ , o sistema tem solução única e pode ser resolvido diretamente, de cima para baixo, isto é:

	$\displaystyle x_{1}=\frac{b_{1}}{a_{1,1}}$		(3.17)
	$\displaystyle x_{2}=\frac{b_{2}-a_{2,1}x_{1}}{a_{2,2}}$		(3.18)
	$\displaystyle\text{}\quad\vdots$		(3.19)
	$\displaystyle x_{n}=\frac{b_{n}-a_{n,1}x_{1}-a_{n,2}x_{2}-\ldots-a_{n,n-1}x_{n% -1}}{a_{n,n}}$		(3.20)

Exemplo 3.1.1.

Vamos resolver o sistema triangular inferior

\begin{matrix}x_{1}&&&=2\\ -3x_{1}&+2x_{2}&&=-8\\ -x_{1}&+x_{2}&-x_{3}&=0\end{matrix}

(3.21)

Na primeira equação, temos

x_{1}=2

(3.22)

Então, da segunda equação do sistema

	$\displaystyle-3x_{1}+2x_{2}=-8$		(3.23)
	$\displaystyle x_{2}=\frac{-8+3x_{1}}{2}$		(3.24)
	$\displaystyle x_{2}=\frac{-8+3\cdot 2}{2}$		(3.25)
	$\displaystyle x_{2}=-1$		(3.26)

E, por fim, da última equação

	$\displaystyle-x_{1}+x_{2}-x_{3}=0$		(3.27)
	$\displaystyle x_{3}=\frac{x_{1}-x_{2}}{-1}$		(3.28)
	$\displaystyle x_{3}=\frac{2-(-1)}{-1}$		(3.29)
	$\displaystyle x_{3}=-3$		(3.30)

Concluímos que a solução do sistema é $\bm{x}=(2,-1,-3)$ .

Código 4: solSisTriaInf.py

⬇

1import numpy as np

3def solSisTriaInf(A, b):

4 n = b.size

5 x = np.zeros_like(b)

6 for i in range(n):

7 x[i] = b[i]

8 for j in range(i):

9 x[i] -= A[i,j]*x[j]

10 x[i] /= A[i,i]

11 return x

13# mat coefficientes

14A = np.array([[1., 0., 0.],

15 [-3., 2., 0.],

16 [-1., 1., -1.]])

18# vet termos constantes

19b = np.array([2., -8., 0.])

21# resol sis lin

22x = solSisTriaInf(A, b)

23print(x)

Observação 3.1.1.(Número de operações em ponto flutuante)

A computação da solução de um sistema $n\times n$ triangular inferior requer $O(n^{2})$ operações em ponto flutuante (multiplicações/divisões e adições/subtrações).

Sistema triangular superior

Um sistema linear triangular superior tem a forma algébrica

\begin{matrix}a_{1,1}x_{1}&+a_{1,2}x_{2}&+\cdots&+a_{1,n}x_{n}&=b_{1}\\ &a_{2,2}x_{2}&+\cdots&+a_{2,n}x_{n}&=b_{2}\\ &&&\vdots&\vdots\\ &&&a_{n,n}x_{n}&=b_{n}\end{matrix}

(3.31)

Pode ser diretamente resolvido de baixo para cima, i.e.

	$\displaystyle x_{n}=\frac{b_{n}}{a_{n,n}}$		(3.32)
	$\displaystyle x_{n-1}=\frac{b_{n-1}-a_{n-1,n}x_{n}}{a_{n-1,n-1}}$		(3.33)
	$\displaystyle\text{}\quad\vdots$		(3.34)
	$\displaystyle x_{n}=\frac{b_{1}-a_{1,1}x_{1}-a_{1,2}x_{2}-\ldots-a_{1,n-1}x_{n% -1}}{a_{1,1}}$		(3.35)

Exemplo 3.1.2.

Vamos resolver o sistema triangular superior

\begin{matrix}2x_{1}&-x_{2}&+2x_{3}&=7\\ &2x_{2}&-x_{3}&=-3\\ &&3x_{3}&=3\end{matrix}

(3.36)

Da última equação, temos

	$\displaystyle x_{3}=\frac{3}{3}$		(3.37)
	$\displaystyle x_{3}=1$		(3.38)

Então, da segunda equação do sistema, obtemos

	$\displaystyle 2x_{2}-x_{3}=-3$		(3.39)
	$\displaystyle x_{2}=\frac{x_{3}-3}{2}$		(3.40)
	$\displaystyle x_{2}=\frac{1-3}{2}$		(3.41)
	$\displaystyle x_{2}=-1$		(3.42)

Por fim, da primeira equação

	$\displaystyle 2x_{1}-x_{2}+2x_{3}=7$		(3.43)
	$\displaystyle x_{1}=\frac{7+x_{2}-2x_{3}}{2}$		(3.44)
	$\displaystyle x_{1}=\frac{7-1-2}{2}$		(3.45)
	$\displaystyle x_{1}=2$		(3.46)

Concluímos que a solução do sistema é $\bm{x}=(2,-1,1)$ .

Código 5: solSisTriaSup.py

⬇

1import numpy as np

3def solSisTriaSup(A, b):

4 n = b.size

5 x = np.zeros_like(b)

6 for i in range(n-1,-1,-1):

7 x[i] = b[i]

8 for j in range(n-1,i,-1):

9 x[i] -= A[i,j]*x[j]

10 x[i] /= A[i,i]

11 return x

14# mat coefficientes

15A = np.array([[2., -1., 2.],

16 [0., 2., -1.],

17 [0., 0., 3.]])

19# vet termos constantes

20b = np.array([7., -3., 3.])

22# sol sis lin

23x = solSisTriaSup(A, b)

24print(x)

Observação 3.1.2.(Número de operações em ponto flutuante)

A computação da solução um sistema $n\times n$ triangular superior requer $O(n^{2})$ operações em ponto flutuante (multiplicações/divisões e adições/subtrações).

3.1.2 Decomposição LU

O procedimento da decomposição LU é equivalente ao método de eliminação gaussiana. Consideramos uma matriz $A=[a_{ij}]_{i,j=1}^{n,n}$ , com $a_{1,1}\neq 0$ , e começamos denotando esta matriz por $U^{(0)}=[u_{i,j}^{(0)}]_{i,j=1}^{n,n}=A$ e tomando $L^{(0)}=I_{n\times n}$ . A eliminação abaixo do pivô $u_{1,1}^{(0)}$ , pode ser computada com as seguintes operações equivalentes por linha

U^{(0)}=\begin{bmatrix}u_{1,1}^{(0)}&u_{1,2}^{(0)}&u_{1,3}^{(0)}&\ldots&u_{1,n% }^{(0)}\\ u_{2,1}^{(0)}&u_{2,2}^{(0)}&u_{2,3}^{(0)}&\ldots&u_{2,n}^{(0)}\\ u_{3,1}^{(0)}&u_{3,2}^{(0)}&u_{3,3}^{(0)}&\ldots&u_{3,n}^{(0)}\\ \vdots&\vdots&\vdots&\ldots&\vdots\\ u_{n,1}^{(0)}&u_{n,2}^{(0)}&a_{n,3}^{(0)}&\ldots&u_{n,n}^{(0)}\end{bmatrix}% \begin{array}[]{l}U_{1}^{(1)}\leftarrow U_{1}^{(0)}\\ U_{2}^{(1)}\leftarrow U_{2}^{(0)}-l_{2,1}^{(1)}U_{1}^{(0)}\\ U_{3}^{(1)}\leftarrow U_{3}^{(0)}-l_{3,1}^{(1)}U_{1}^{(0)}\\ \vdots\\ U_{n}^{(1)}\leftarrow U_{n}^{(0)}-l_{n,1}^{(1)}U_{1}^{(0)}\end{array}

(3.47)

onde, $l_{i,1}^{(1)}=u_{i,1}^{(0)}/u_{1,1}^{(0)}$ , $i=2,3,\dotsc,n$ .

Destas computações, obtemos uma nova matriz da forma

U^{(1)}=\begin{bmatrix}u_{1,1}^{(1)}&u_{1,2}^{(1)}&u_{1,3}^{(1)}&\ldots&u_{1,n% }^{(1)}\\ 0&u_{2,2}^{(1)}&u_{2,3}^{(1)}&\ldots&u_{2,n}^{(1)}\\ 0&u_{3,2}^{(1)}&u_{3,3}^{(1)}&\ldots&u_{3,n}^{(1)}\\ \vdots&\vdots&\vdots&\ldots&\vdots\\ 0&u_{n,2}^{(1)}&u_{n,3}^{(1)}&\ldots&u_{n,n}^{(1)}\end{bmatrix}

(3.48)

E, denotando

L^{(1)}=\begin{bmatrix}1&0&0&\ldots&0\\ l_{2,1}^{(1)}&1&0&\ldots&0\\ l_{3,1}^{(1)}&0&1&\ldots&0\\ \vdots&\vdots&\vdots&\ldots&\vdots\\ l_{n,1}^{(1)}&0&0&\ldots&1\end{bmatrix}

(3.49)

temos

A=L^{(1)}U^{(1)}.

(3.50)

No caso de $u^{(1)}_{2,2}\neq 0$ , podemos continuar com o procedimento de eliminação gaussiana com as seguintes operações equivalentes por linha

U^{(1)}=\begin{bmatrix}u_{1,1}^{(1)}&u_{1,2}^{(1)}&u_{1,3}^{(1)}&\ldots&u_{1,n% }^{(1)}\\ 0&u_{2,2}^{(1)}&u_{2,3}^{(1)}&\ldots&u_{2,n}^{(1)}\\ 0&u_{3,2}^{(1)}&u_{3,3}^{(1)}&\ldots&u_{3,n}^{(1)}\\ \vdots&\vdots&\vdots&\ldots&\vdots\\ 0&u_{n,2}^{(1)}&u_{n,3}^{(1)}&\ldots&u_{n,n}^{(1)}\end{bmatrix}\begin{array}[]% {l}U_{1}^{(2)}\leftarrow U_{1}^{(1)}\\ U_{2}^{(2)}\leftarrow U_{2}^{(1)}\\ U_{3}^{(2)}\leftarrow U_{3}^{(1)}-l_{3,2}^{(2)}U_{2}^{(1)}\\ \vdots\\ U_{n}^{(2)}\leftarrow U_{n}^{(1)}-l_{n,2}^{(2)}U_{2}^{(1)}\\ \end{array}

(3.51)

onde, $l_{i,2}^{(2)}=u_{i,2}^{(1)}/u_{2,2}^{(1)}$ , $i=3,4,\ldots,n$ . Isto nos fornece o que nos fornece

U^{(2)}=\begin{bmatrix}u_{1,1}^{(2)}&u_{1,2}^{(2)}&u_{1,3}^{(2)}&\ldots&u_{1,n% }^{(2)}\\ 0&u_{2,2}^{(2)}&u_{2,3}^{(2)}&\ldots&u_{2,n}^{(2)}\\ 0&0&u_{3,3}^{(2)}&\ldots&u_{3,n}^{(2)}\\ \vdots&\vdots&\vdots&\ldots&\vdots\\ 0&0&u_{n,3}^{(2)}&\ldots&u_{n,n}^{(2)}\end{bmatrix}.

(3.52)

Além disso, denotando

L^{(2)}=\begin{bmatrix}1&0&0&\ldots&0\\ l_{2,1}^{(1)}&1&0&\ldots&0\\ l_{3,1}^{(1)}&l_{3,2}^{(2)}&1&\ldots&0\\ \vdots&\vdots&\vdots&\ldots&\vdots\\ l_{n,1}^{(1)}&l_{n,2}^{(2)}&0&\ldots&1\end{bmatrix}

(3.53)

temos

A=L^{(2)}U^{(2)}.

(3.54)

Continuando com este procedimento, ao final de $n-1$ passos teremos obtido a decomposição

A=LU,

(3.55)

onde $L$ é a matriz triangular inferior

L=L^{(n-1)}=\begin{bmatrix}1&0&0&\ldots&0\\ l_{2,1}^{(1)}&1&0&\ldots&0\\ l_{3,1}^{(1)}&l_{3,2}^{(2)}&1&\ldots&0\\ \vdots&\vdots&\vdots&\ldots&\vdots\\ l_{n,1}^{(1)}&l_{n,2}^{(2)}&l_{n,3}^{(3)}&\ldots&1\end{bmatrix}

(3.56)

e $U$ é a matriz triangular superior

U=U^{(n-1)}=\begin{bmatrix}u_{1,1}^{(0)}&u_{1,2}^{(0)}&u_{1,3}^{(0)}&\ldots&u_% {1,n}^{(0)}\\ 0&u_{2,2}^{(1)}&u_{2,3}^{(1)}&\ldots&u_{2,n}^{(1)}\\ 0&0&u_{3,3}^{(2)}&\ldots&u_{3,n}^{(2)}\\ \vdots&\vdots&\vdots&\ldots&\vdots\\ 0&0&0&\ldots&u_{n,n}^{(n-1)}\end{bmatrix}.

(3.57)

Exemplo 3.1.3.

Consideramos a seguinte matriz

A=\begin{bmatrix}-1&2&-2\\ 3&-4&1\\ 1&-5&3\end{bmatrix}.

(3.58)

Para obtermos sua decomposição $L U$ começamos com

	$\displaystyle L^{(0)}$	$\displaystyle:=\begin{bmatrix}1&0&0\\ 0&1&0\\ 0&0&1\end{bmatrix},$		(3.59)
	$\displaystyle U^{(0)}$	$\displaystyle:=\begin{bmatrix}-1&2&-2\\ 3&-4&1\\ 1&-5&3\end{bmatrix}.$		(3.60)

Então, observando que a eliminação abaixo do pivô $u_{1,1}=-1$ pode ser feita com as seguintes operações equivalentes por linha

U^{(0)}=\begin{bmatrix}-1&2&-2\\ 3&-4&1\\ 1&-5&3\end{bmatrix}\begin{array}[]{l}U_{1}^{(1)}\leftarrow U_{1}^{(0)}\\ U_{2}^{1}\leftarrow U_{2}^{(0)}-\frac{3}{-1}U_{1}^{(0)}\\ U_{3}^{1}\leftarrow U_{3}^{(0)}-\frac{1}{-1}U_{1}^{(0)}\end{array},

(3.61)

temos

	$\displaystyle L^{(1)}$	$\displaystyle:=\begin{bmatrix}1&0&0\\ -3&1&0\\ -1&0&1\end{bmatrix},$		(3.62)
	$\displaystyle U^{(1)}$	$\displaystyle:=\begin{bmatrix}-1&2&-2\\ 0&2&-5\\ 0&-3&1\end{bmatrix}.$		(3.63)

Agora, para eliminarmos abaixo do pivô $u_{2,2}=2$ , usamos as operações

U^{(1)}:=\begin{bmatrix}-1&2&-2\\ 0&2&-5\\ 0&-3&1\end{bmatrix}\begin{array}[]{l}U_{1}^{(2)}\leftarrow U_{1}^{(1)}\\ U_{2}^{(2)}\leftarrow U_{2}^{(1)}\\ U_{3}^{(2)}\leftarrow U_{3}^{(1)}-\frac{-3}{2}U_{2}^{(1)}\end{array},

(3.64)

donde

	$\displaystyle L^{(2)}=\begin{bmatrix}1&0&0\\ -3&1&0\\ -1&-1.5&1\end{bmatrix},$		(3.65)
	$\displaystyle U^{(1)}=\begin{bmatrix}-1&2&-2\\ 0&2&-5\\ 0&0&-6.5\end{bmatrix}.$		(3.66)

Isso completa a decomposição, sendo $L:=L^{(3)}$ e $U:=U^{(3)}$ .

Código 6: lu.py

⬇

1import numpy as np

3def lu(A):

4 # num. de linhas

5 n = A.shape[0]

7 # inicialização

8 U = A.copy()

9 L = np.eye(n)

11 # decomposição

12 for i in range(n-1):

13 for j in range(i+1,n):

14 L[j,i] = U[j,i]/U[i,i]

15 U[j,i:] -= L[j,i]*U[i,i:]

17 return L, U

19# matriz

20A = np.array([[-1., 2., -2.],

21 [3., -4., 1.],

22 [1., -5., 3.]])

24L, U = lu(A)

Observação 3.1.3.(Número de operações em ponto flutuante)

A decomposição LU de um sistema $n\times n$ requer $O(n^{3})$ operações em ponto flutuante (multiplicações/divisões e adições/subtrações).

3.1.3 Resolução do sistema com decomposição LU

Consideramos o sistema linear

A\bm{x}=\bm{b}.

(3.67)

Para resolvê-lo com o Método LU, fazemos

1.

Computamos a decomposição LU

$\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}A=LU.$ (3.68)
2.

Resolvemos o sistema triangular inferior

$\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}L\bm{y}=\bm{b}.$ (3.69)
3.

Resolvemos o sistema triangular superior

$\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}U\bm{x}=\bm{y}.$ (3.70)

Exemplo 3.1.4.

Vamos resolver o seguinte sistema linear

	$\displaystyle-x_{1}+2x_{2}-2x_{3}=6$		(3.71)
	$\displaystyle 3x_{1}-4x_{2}+x_{3}=-11$		(3.72)
	$\displaystyle x_{1}-5x_{2}+3x_{3}=-10.$		(3.73)

No exemplo anterior (Exemplo 3.1.4), vimos que a matriz de coeficientes $A$ deste sistema admite a seguinte decomposição LU

\underbrace{\begin{bmatrix}-1&2&-2\\ 3&-4&1\\ 1&-5&3\end{bmatrix}}_{A}=\underbrace{\begin{bmatrix}1&0&0\\ -3&1&0\\ -1&-1.5&1\end{bmatrix}}_{L}\underbrace{\begin{bmatrix}-1&2&-2\\ 0&2&-5\\ 0&0&-6.5\end{bmatrix}}_{U}

(3.74)

Daí, iniciamos resolvendo o seguinte sistema triangular inferior $Ly=b$ , i.e.

	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% y_{1}=6$		(3.75)
	$\displaystyle-3y_{1}+y_{2}=-11$		(3.76)
	$\displaystyle\text{}\qquad\Rightarrow\color[rgb]{0,0,1}\definecolor[named]{% pgfstrokecolor}{rgb}{0,0,1}y_{2}=7,$		(3.77)
	$\displaystyle-y_{1}-1.5y_{2}+y_{3}=-10$		(3.78)
	$\displaystyle\text{}\qquad\Rightarrow\color[rgb]{0,0,1}\definecolor[named]{% pgfstrokecolor}{rgb}{0,0,1}y_{3}=6.5.$		(3.79)

Por fim, computamos a solução $x$ resolvendo o sistema triangular superior $Ux=y$ , i.e.

	$\displaystyle-6.5x_{3}=6.5$		(3.80)
	$\displaystyle\text{}\qquad\Rightarrow\color[rgb]{0,0,1}\definecolor[named]{% pgfstrokecolor}{rgb}{0,0,1}x_{3}=-1,$		(3.81)
	$\displaystyle 2x_{2}-5x_{3}=7$		(3.82)
	$\displaystyle\text{}\qquad\Rightarrow\color[rgb]{0,0,1}\definecolor[named]{% pgfstrokecolor}{rgb}{0,0,1}x_{2}=1$		(3.83)
	$\displaystyle-x_{1}+2x_{2}-2x_{3}=6$		(3.84)
	$\displaystyle\text{}\qquad\Rightarrow\color[rgb]{0,0,1}\definecolor[named]{% pgfstrokecolor}{rgb}{0,0,1}x_{1}=-2.$		(3.85)

⬇

1import numpy as np

3# mat coefs

4A = np.array([[-1., 2., -2.],

5 [3., -4., 1.],

6 [1., -5., 3.]])

7# vet termos const

8b = np.array([6., -11., -10.])

10# 1. LU

11L, U = lu(A)

13# 2. Ly = b

14y = solSisTriaInf(L, b)

16# 3. Ux = y

17x = solSisTriaSup(U, y)

3.1.4 Fatoração LU com pivotamento parcial

O algoritmo estudando acima não é aplicável no caso de o pivô ser nulo, o que pode ser corrigido através de permutações de linhas. Na fatoração LU com pivotamento parcial fazemos permutações de linha na matriz de forma que o pivô seja sempre aquele de maior valor em módulo. Por exemplo, suponha que o elemento $a_{3,1}$ seja o maior valor em módulo na primeira coluna da matriz $A=U^{(0)}$ com

U^{(0)}=\begin{bmatrix}u_{1,1}^{(0)}&u_{1,2}^{(0)}&u_{1,3}^{(0)}&\ldots&u_{1,n% }^{(0)}\\ u_{2,1}^{(0)}&u_{2,2}^{(0)}&u_{2,3}^{(0)}&\ldots&u_{2,n}^{(0)}\\ \bm{u_{3,1}^{(0)}}&u_{3,2}^{(0)}&u_{3,3}^{(0)}&\ldots&u_{3,n}^{(0)}\\ \vdots&\vdots&\vdots&\ldots&\vdots\\ u_{n,1}^{(0)}&u_{n,2}^{(0)}&a_{n,3}^{(0)}&\ldots&u_{n,n}^{(0)}\end{bmatrix}.

(3.86)

Neste caso, o procedimento de eliminação na primeira coluna deve usar $u_{3,1}^{(0)}$ como pivô, o que requer a permutação entre as linhas $1$ e $3$ ( $U_{1}^{(0)}\leftrightarrow U_{3}^{(0)}$ ). Isto pode ser feito utilizando-se da seguinte matriz de permutação

P=\begin{bmatrix}0&0&1&\ldots&0\\ 0&1&0&\ldots&0\\ 1&0&0&\ldots&0\\ \vdots&\vdots&\vdots&\ldots&\vdots\\ 0&0&0&\ldots&1\end{bmatrix}.

(3.87)

Com essa, iniciamos o procedimento de decomposição LU com $PA=L^{(0)}U^{(0)}$ , onde $L^{(0)}=I_{n\times n}$ e $U^{(0)}=PA$ . Caso sejam necessárias outras mudanças de linhas no decorrer do procedimento de decomposição, a matriz de permutação $P$ deve ser atualizada apropriadamente.

Exemplo 3.1.5.

Vamos fazer a decomposição LU com pivotamento parcial da seguinte matriz

A=\begin{bmatrix}-1&2&-2\\ 3&-4&1\\ 1&-5&3\end{bmatrix}

(3.88)

Começamos, tomando

	$\displaystyle P^{(0)}=\begin{bmatrix}1&0&0\\ 0&1&0\\ 0&0&1\end{bmatrix},$		(3.89)
	$\displaystyle L^{(0)}=\begin{bmatrix}1&0&0\\ 0&1&0\\ 0&0&1\end{bmatrix},$		(3.90)
	$\displaystyle U^{(0)}=\begin{bmatrix}-1&2&-2\\ 3&-4&1\\ 1&-5&3\end{bmatrix}$		(3.91)

O candidato a pivô é o elemento $u^{(0)}_{2,1}$ . Então, fazemos as permutações de linhas

	$\displaystyle P_{1}^{((0)}$	$\displaystyle\leftrightarrow P_{2}^{(0)},$		(3.92)
	$\displaystyle U_{1}^{(0)}$	$\displaystyle\leftrightarrow U_{2}^{(0)}$		(3.93)

e, na sequência, as operações elementares por linhas

U_{2:3}^{(1)}\leftarrow U_{2:3}^{(0)}-m^{(0)}_{2:3,1}U_{1}^{(0)},

(3.94)

donde obtemos

	$\displaystyle P^{(1)}=\begin{bmatrix}0&1&0\\ 1&0&0\\ 0&0&1\end{bmatrix},$		(3.95)
	$\displaystyle L^{(1)}=\begin{bmatrix}1&0&0\\ -0,\overline{3}&1&0\\ 0,\overline{3}&0&1\end{bmatrix},$		(3.96)
	$\displaystyle U^{(1)}=\begin{bmatrix}3&-4&1\\ 0&0,\overline{6}&-1,\overline{6}\\ 0&-3,\overline{6}&2,\overline{6}\end{bmatrix}$		(3.97)

Agora, o candidato a pivô é o elemento $u^{(1)}_{3,2}$ . Assim, fazemos as permutações de linhas

	$\displaystyle P^{(1)}_{2}\leftrightarrow P^{(1)}_{3},$		(3.98)
	$\displaystyle U^{(1)}_{2}\leftrightarrow U^{(1)}_{3}$		(3.99)

e análogo para os elementos da coluna 1 de $L$ . Então, fazemos a operação elementar por linha

U^{(2)}_{3}\leftarrow U^{(1)}_{3}-m^{(1)}_{3,2}U^{(1)}_{2}

(3.100)

. Com isso, obtemos

	$\displaystyle P^{(2)}=\begin{bmatrix}0&1&0\\ 0&0&1\\ 1&0&0\end{bmatrix},$		(3.101)
	$\displaystyle L^{(2)}=\begin{bmatrix}1&0&0\\ 0,\overline{3}&1&0\\ -0,\overline{3}&-0,\overline{18}&1\end{bmatrix},$		(3.102)
	$\displaystyle U^{(2)}=\begin{bmatrix}3&-4&1\\ 0&-3,\overline{6}&2,\overline{6}\\ 0&0&-1,\overline{18}\end{bmatrix}$		(3.103)

Por fim, temos obtido a decomposição LU de $A$ na forma

PA=LU,

(3.104)

com $P=P^{(2)}$ , $L=L^{(2)}$ e $U=U^{(2)}$ .

Código 7: lup.py

⬇

1import numpy as np

3def lup(A):

4 # num. de linhas

5 n = A.shape[0]

7 # inicialização

8 U = A.copy()

9 L = np.eye(n)

10 P = np.eye(n)

12 # decomposição

13 for i in range(n-1):

14 # permutação de linhas

15 p = i + np.argmax(np.fabs(U[i:,i]))

16 P[[i,p]] = P[[p,i]]

17 U[[i,p]] = U[[p,i]]

18 L[[i,p],:i] = L[[p,i],:i]

19 # eliminação gaussiana

20 for j in range(i+1,n):

21 L[j,i] = U[j,i]/U[i,i]

22 U[j,i:] -= L[j,i]*U[i,i:]

24 return P, L, U

26# matriz

27A = np.array([[-1., 2, -2],

28 [3, -4, 1],

29 [1, -5, 3]])

31P, L, U = lup(A)

Exemplo 3.1.6.

Vamos computar a solução do seguinte sistema linear com o Método da Decomposição LU com Pivotamento Parcial.

	$\displaystyle-x_{1}+2x_{2}-2x_{3}=6$		(3.105)
	$\displaystyle 3x_{1}-4x_{2}+x_{3}=-11$		(3.106)
	$\displaystyle x_{1}-5x_{2}+3x_{3}=-10.$		(3.107)

No exemplo anterior (Exemplo 3.1.5), vimos que a matriz de coeficientes $A$ deste sistema admite a seguinte decomposição LU

PA=LU

(3.108)

com

	$\displaystyle P=\begin{bmatrix}0&1&0\\ 0&0&1\\ 1&0&0\end{bmatrix},$		(3.109)
	$\displaystyle L=\begin{bmatrix}1&0&0\\ 0,\overline{3}&1&0\\ -0,\overline{3}&-0,\overline{18}&1\end{bmatrix},$		(3.110)
	$\displaystyle U=\begin{bmatrix}3&-4&1\\ 0&-3,\overline{6}&2,\overline{6}\\ 0&0&-1,\overline{18}\end{bmatrix}$		(3.111)

Multiplicando o sistema a esquerda pela matriz $P$ , obtemos

	$\displaystyle PA\bm{x}=P\bm{b}$		(3.112)
	$\displaystyle LU\bm{x}=P\bm{b}$		(3.113)
	$\displaystyle L(U\bm{x})=P\bm{b}$		(3.114)

Com isso, resolvemos

L\bm{y}=P\bm{b}

(3.115)

donde obtemos

\bm{y}=(-11,-6,\overline{3},1,\overline{18})

(3.116)

Então, resolvemos

U\bm{x}=y

(3.117)

donde obtemos a solução

\bm{x}=(-2,1,-1).

(3.118)

⬇

1# matriz

2A = np.array([[-1., 2., -2.],

3 [3., -4., 1.],

4 [1., -5., 3.]])

5b = np.array([6,-11,-10])

7P, L, U = lup(A)

9y = solSisTriaInf(L, P@b)

10x = solSisTriaSup(U, y)

3.1.5 Exercícios

E. 3.1.1.

Seja a matriz

A=\begin{bmatrix}-1&2&-2\\ 3&4&1\\ -4&-5&3\end{bmatrix}

(3.119)

a)

Compute sua decomposição LU sem pivotamento parcial.
b)

Compute sua decomposição LU com pivotamento parcial.

	$\displaystyle L=\begin{bmatrix}1&0&0\\ -3&1&0\\ 4&0&1\end{bmatrix}$		(3.120)
	$\displaystyle U=\begin{bmatrix}-1&2&2\\ 0&10&-5\\ 0&0&4.5\end{bmatrix}$		(3.121)

	$\displaystyle P=\begin{bmatrix}0&0&1\\ 1&0&0\\ 0&1&0\end{bmatrix}$		(3.122)
	$\displaystyle L=\begin{bmatrix}0&0&0\\ 0.25&1&0\\ -0.75&7.692\mathrm{e}\!-\!2&1\end{bmatrix}$		(3.123)
	$\displaystyle U=\begin{bmatrix}-4&-5&3\\ 0&3.25&-2.75\\ 0&0&3.462\end{bmatrix}$		(3.124)

E. 3.1.2.

Use o Método da Decomposição LU para resolver o sistema linear

	$\displaystyle-x_{1}+2x_{2}-2x_{3}=-1$		(3.125)
	$\displaystyle 3x_{1}-4x_{2}+x_{3}=-4$		(3.126)
	$\displaystyle-4x_{1}-5x_{2}+3x_{3}=20$		(3.127)

usando LU.

$x_{1}=-3$ , $x_{2}=-1$ , $x_{3}=1$

E. 3.1.3.

Compute a decomposição LU da matriz

A=\begin{bmatrix}-1&0&-2&1\\ 3&0&0&-2\\ 1&-1&0&-1\\ 0&2&-3&0\end{bmatrix}.

(3.128)

	$\displaystyle P=\begin{bmatrix}0&1&0&0\\ 0&0&0&1\\ 1&0&0&0\\ 0&0&1&0\end{bmatrix}$		(3.129)
	$\displaystyle L=\begin{bmatrix}1&0&0&0\\ 0&1&0&0\\ -0.\overline{3}&0&1&0\\ 0.\overline{3}&-0.5&0.75&1\end{bmatrix}$		(3.130)
	$\displaystyle U=\begin{bmatrix}3&0&0&-2\\ 0&2&-3&0\\ 0&0&-2&0.\overline{3}\\ 0&0&0&-0.58\overline{3}\end{bmatrix}$		(3.131)

E. 3.1.4.

Use o Método da Decomposição LU para resolver o sistema linear

	$\displaystyle-x_{1}-2x_{3}+x_{4}=1$		(3.132)
	$\displaystyle 3x_{1}-2x_{4}=-7$		(3.133)
	$\displaystyle x_{1}-x_{2}-x_{4}=-3$		(3.134)
	$\displaystyle 2x_{2}-3x_{3}=-3$		(3.135)

$x=(-1,0,1,2)$

E. 3.1.5.

A matriz de Vandermonde¹¹1Alexandre-Théophile Vandermonde, 1735 - 1796, matemático francês. Fonte: Wikipédia: Alexandre-Theóphile Vandermonde. é definida por

V(x_{1},\ldots,x_{n})=\begin{bmatrix}x_{1}^{n-1}&\cdots&x_{1}^{2}&x_{1}&1\\ x_{2}^{n-1}&\cdots&x_{2}^{2}&x_{2}&1\\ \vdots&&\vdots&\vdots&\vdots\\ x_{n}^{n-1}&\cdots&x_{n}^{2}&x_{n}&1\\ \end{bmatrix}

(3.136)

Compute a decomposição LU com pivotamento parcial das matrizes de Vandermonde²²2Alexandre-Théophile Vandermonde, 1735 - 1796, matemático francês. Fonte: Wikipédia: Alexandre-Theóphile Vandermonde..

a)

$V(1,2,3)$
b)

$V(-2,-1,0,1)$
c)

$V(0.1,0.25,0.5,0.75)$
d)

$V(-0.1,0.5,1,2,10)$

Dica: A matriz de Vandermonde pode ser alocada com o seguinte código:

⬇

1import numpy as np

2x = np.array([1.,2,3])

3n = x.size

4V = np.ones((n,n))

5for i in range(n-1):

6 V[:,i] = x**(n-1-i)

7print(V)

E. 3.1.6.

Use o Método da Decomposição LU para computar a matriz inversa de cada uma das seguintes matrizes:

a)

$A=\begin{bmatrix}-1&2&-2\\ 3&4&1\\ -4&-5&3\end{bmatrix}$ (3.137)
b)

$B=\begin{bmatrix}-1&0&-2&1\\ 3&0&0&-2\\ 1&-1&0&-1\\ 0&2&-3&0\end{bmatrix}.$ (3.138)

Dica: $AA^{-1}=I$ .

a)

$A^{-1}=\begin{bmatrix}-0.3778&-0.0889,&-0.2222\\ 0.2889&0.2444&0.1111\\ -0.0222&0.2889&0.2222\end{bmatrix}$ (3.139)
b)

$B^{-1}=\begin{bmatrix}0.8571&1&-1.1429&-0.5714\\ -0.4286&0&-0.4286&0.2857\\ -0.2857&0&-0.2857&-0.1429\\ 1.2857&1.&-1.7143&-0.8571\end{bmatrix}$ (3.140)

Envie seu comentário

Aproveito para agradecer a todas/os que de forma assídua ou esporádica contribuem enviando correções, sugestões e críticas!

Este texto é disponibilizado nos termos da Licença Creative Commons Atribuição-CompartilhaIgual 4.0 Internacional. Ícones e elementos gráficos podem estar sujeitos a condições adicionais.

Política de uso de dados

Política de uso de dados

Matemática numérica I

3.1 Método da decomposição LU

3.1.1 Sistemas triangulares

Sistema triangular inferior

Exemplo 3.1.1.

Observação 3.1.1.(Número de operações em ponto flutuante)

Sistema triangular superior

Exemplo 3.1.2.

Observação 3.1.2.(Número de operações em ponto flutuante)

3.1.2 Decomposição LU

Exemplo 3.1.3.

Observação 3.1.3.(Número de operações em ponto flutuante)

3.1.3 Resolução do sistema com decomposição LU

Exemplo 3.1.4.

3.1.4 Fatoração LU com pivotamento parcial

Exemplo 3.1.5.

Exemplo 3.1.6.

3.1.5 Exercícios

E. 3.1.1.

E. 3.1.2.

E. 3.1.3.

E. 3.1.4.

E. 3.1.5.

E. 3.1.6.

Envie seu comentário