Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

Reproduzir o algoritmo da seção 2.5.3 do capítulo 2 do livro Scientific Computing: An Introductory Survey - Michael T Heath

0 votos
10 visitas
perguntada Nov 19 em Ciência da Computação por JOAO PAULO (1 ponto)  
editado Nov 19 por JOAO PAULO

O exercício 2.12 pede para reproduzirmos o código da seção 2.5.3 e testar esse código em alguns sistemas. Também pergunta a alteração no algoritmo caso seja necessário aplicar "partial pivoting". O último ponto questiona-se o que altera no algoritmo caso a matriz seja positiva definida e seja realizada a decomposição de Cholesky.

Compartilhe

1 Resposta

0 votos
respondida Nov 19 por JOAO PAULO (1 ponto)  
editado Nov 20 por JOAO PAULO

Uma matriz banda é definida como uma matriz uma matriz esparsa cujas entradas diferentes de zero estão confinadas a uma banda diagonal, compreendendo a diagonal principal e zero ou mais diagonais em cada lado. Formalmente temos que:

A imagem será apresentada aqui.

A matriz tridiagonal é um caso particular da matriz banda, onde β=1. Ou seja, uma matriz tridiagonal é uma matriz que possui elementos diferentes de zero apenas na sua diagonal principal e nas diagonais acima e abaixo da diagonal principal. Todos os elementos restantes da matriz resultante são 0.
Uma característica interessante dessa matriz é que com esse tipo de matriz temos a possibilidade de fazermos a decomposição LU de maneira muito mais eficiente do que faríamos com outros tipos de matrizes não singulares. Essa eficiência advêm do fato de não precisarmos “levar em conta” os inúmeros zeros ao fazermos a decomposição LU. Esse é exatamente o sentido do algoritmo da seção 2.5.3. Que é definido da seguinte maneira:

A imagem será apresentada aqui.

Abaixo eu mostro a implementação desse código no Python.

import numpy as np

def lu_trid(A):

    n = A.shape[0]

    U = np.zeros((n, n), dtype=np.double)
    L = np.eye(n, dtype=np.double)

    U[0][0]=A[0][0]
    for k in range(n-1):

        L[k+1,:(k+1)] = A[k+1,:(k+1)]/ U[k, k]

        U[k+1,k+1] = A[(k+1),(k+1)] - L[k+1:,(k)] @ A[k,(k+1):]
        U[k,k+1:] = A[k,k+1:]



    return L, U

if __name__ == '__main__':
#    A = np.array([[1,4,0,0], [3,4,1,0], [0,2,3,4], [0,0,1,3]])
    A = np.array([[2,3,0], [1,2,3], [0,1,2]])
    L, U = lu_trid(A)
    print("matriz A", A)
    print("matriz L",L)
    print("matriz U",U)
    print("total check",L@U)

No próximo passo, o exercício pergunta como o código será alterado caso consideremos “partial pivoting”. Geralmente as matrizes tridiagonais são positivas definidas ou diagonalmente definidas, logo partial pivoting não é necessário aplicar nesses casos, a não ser que o sistema seja instável. No entanto, se o sistema não positivo definido ou ter alguma instabilidade, será necessário considerarmos partial pivoting. Nesse caso, o algoritmo fica mais complexo pois a matriz A inicial pode deixar de ser banda. Por exemplo: considere a seguinte matriz tridiagonal:

\[A= \begin{bmatrix}2&3&0\\ 3&2&3\\ 0&1&2\end{bmatrix}\]

Ao aplicarmos a matriz P de ”partial pivoting” à matriz A, a matriz resultante PA será:

\[PA=\begin{bmatrix}3&2&3\\ 2&3&0\\ 0&1&2\end{bmatrix}\]

O que claramente deixa de ser tridiagonal. Nesse caso podemos fazer uma decomposição LU tradicional com partial pivoting. É importante ressaltar que esse algoritmo PLU será menos eficiente do que o algoritmo da seção anterior. Esse código pode ser visto abaixo:

import numpy as np

def PLU(A):

    n = A.shape[0]

    P = np.eye(n, dtype=np.double)
    U = np.zeros((n, n), dtype=np.double)
    L = np.eye(n, dtype=np.double)

    for k in range(n-1):
        max_row_index = np.argmax(abs(A[k:n,k]))+k
        P[[k,max_row_index]]=P[[max_row_index,k]]
        A[[k,max_row_index]]=A[[max_row_index,k]]


    for k in range(n):

        U[k, k:] = A[k, k:] - L[k,:k] @ U[:k,k:]
        L[(k+1):,k] = (A[(k+1):,k] - L[(k+1):,:] @ U[:,k]) / U[k, k]

    return P, L, U

if __name__ == '__main__':    
    A = np.array([[2,3,0], [3,2,3], [0,1,2]])
    P, L, U= PLU(A)
    print("matriz P", P)
    print("matriz U", U)
    print("matriz L", L)
    print("matriz PLU", P@L@U)

Em relação ao último ponto, a decomposição de Cholesky pode ser aplicada se a matriz em questão for simétrica e positiva definida. Para uma matriz simétrica qualquer, a decomposição de Cholesky é mais eficiente do que a decomposição LU padrão (pois na decomposição de Cholesky somente é necessário encontrar uma das matrizes L ou U). No entanto, se estivermos tratando de uma matriz tridiagonal, como é o caso, podemos ter um algoritmo mais eficiente ainda do que o algoritmo considerado na primeira parte do exercício. Abaixo, eu mostro um código que faz essa decomposição (na verdade, é uma versão ligeiramente diferente do código da parte 1 do exercício)

import numpy as np

def chol_trid(A):

    n = A.shape[0]
    L = np.zeros((n, n), dtype=np.double)


    L[0,0]=np.sqrt(A[0,0])
    for k in range(n-1):

        L[k+1,k] = A[k,(k+1)]/ L[k, k]
        L[k+1,k+1] = np.sqrt(A[k+1,k+1]- L[k+1,k]**2)

    return L

if __name__ == '__main__':
    A=  np.array([[4,2,0,0,0], [2,5,2,0,0], [0,2,5,2,0],[0,0,2,5,2],[0,0,0,2,5]])
    L = chol_trid(A)
    print(A)
    print(L)
    print("total",L@L.T)
...