Rotina python p/ Inversão de Matrizes pelo método de fatoração QR

Aqui está postado o código em python 2.6.5 que faz a inversão de uma matriz inversível.

# -*- coding: cp1252 -*-

# Inverção da matriz [K] pelo método de fatoração QR

from numpy import *

# Matriz [K]

K = array([[3, 1, 0],[1, 3, 1],[0, 1, 3]],float)

print 'Matriz [K] = '

print K

# Ordem da matriz

n=len(K) # conta o número de colunas da matriz [K]

# Inicia Pg

Pg=zeros([n,n],float)

for i in range(n):

____Pg[i][i]=1

# Fatoração QR (Givens)

for i in range(n-1):

____for h in range(i+1,n):

________if (K[i][i])**2+(K[h][i])**2==0:

____________cs=1

____________sn=0

________else:

____________# Determina cossenos e senos

____________cs=(K[i][i])/((K[i][i])**2+(K[h][i])**2)**0.5

____________sn=(K[h][i])/((K[i][i])**2+(K[h][i])**2)**0.5

________for j in range(n):

____________k1,k2 = 0, 0

____________k1,k2 = K[i][j], K[h][j]

____________# Cálculo da matriz ~K ( vai transformando K em Sg)

____________K[i][j] = k1*cs + k2*sn

____________K[h][j] = -k1*sn + k2*cs

____________# Cálcula Pg transposto

____________p1,p2 = 0, 0

____________p1,p2 = Pg[i][j], Pg[h][j]

____________Pg[i][j] = p1*cs + p2*sn

____________Pg[h][j] = -p1*sn + p2*cs

# Obtém a matriz inversa de [K] armazenando cada vetor coluna ci resultante da solução

#do sistema Sg.ci=li onde li é o vetor linha de Pg.

# No código abaixo K foi tranformado em Sg.

R=zeros([n,1],float)

for k in range(n):

____for i in range(n):

________R[i][0]=Pg[i][k]

____U=zeros([n,1],float)

____for i in range(n):

________SS=0

________for j in range(n):

____________SS=SS+(K[n-1-i][n-1-j]*U[n-1-j])

________U[n-1-i]=(R[n-1-i]-SS)/K[n-1-i][n-1-i]

____for i in range(n):

________Pg[i][k]=U[i][0]# Matriz inversa [K]-1

print Pg

# Obs: Ao copiar o arquivo, antes de rodar a rotina é necessário editar o código substituindo cada underline por um espaço.

Para que esta rotina funcione é necessário que vc tenha instalado em seu computador o compilador python e o modulo de ferramentas numpy. Ambos podem ser baixados do site

http://www.python.org.

Posteriormente, vou editar esta publicação adicionando a teoria que fundamenta o código desta rotina.

Rotina python para solução de equações de 1° e 2° grau

Bom!

Está postagem é a primeira que faço no meu blog.

Digamos que ela é o marco inicial de uma série de postagens que compõem parte do conhecimento técnico que existe dentro da minha cabeça.

Tá certo! vamos direto ao assunto.

Logo abaixo vou postar uma rotina em Python que resolve equações de 1° e 2° grau.

A versão do compilador Python que vou utilizar é a 2.6.5.

1° Passo: Abra o compilador Python 2.6.5

2° Passo: Abra uma nova janela clicando em File > New Window

3° Passo: Escreva o seguinte código

# -*- coding: cp1252 -*-

##################### Solução de equações de 1° e 2° grau ######################

print 'Digite 1 se a equação é de 1°grau e 2 se é de 2° grau'

eq = raw_input('Entre com o número que define o tipo de equação: ')

if eq=='1' or eq=='2':

____eq=int(eq) #Converte o texto em valor numérico inteiro

____if eq==1:

________print 'Formato geral da equação do 1° grau ax+b=0'

________a=raw_input('Digite o valor de a: ')

________a=float(a) #Converte o texto num valor numérico real

________b=raw_input('Digite o valor de b: ')

________b=float(b) #Converte o texto num valor numérico real

________x = -b/a

________print 'x = ', x

____elif eq==2:

________print 'Formato geral da equação do 2° grau ax²+bx+c=0'

________a=raw_input('Digite o valor de a: ')

________a=float(a) #Converte o texto num valor numérico real

________b=raw_input('Digite o valor de b: ')

________b=float(b) #Converte o texto num valor numérico real

________c=raw_input('Digite o valor de c: ')

________c=float(c) #Converte o texto num valor numérico real

________# Cálculo do valor delta

________delta = b**2 - 4*a*c

________if delta>=0:

____________x1=(-b-delta**0.5)/(2*a)

____________x2=(-b+delta**0.5)/(2*a)

____________print 'x1 = ', x1

____________print 'x2 = ', x2

________elif delta<0:

____________raizdelta=complex(0,abs(delta))

____________x1=(-b-raizdelta)/(2*a)

____________x2=(-b+raizdelta)/(2*a)

____________print 'Raízes complexas'

____________print 'x1 = ', x1

____________print 'x2 = ', x2

else:

____print 'Entrada não válida'

Obs.: Ao copiar a rotina, antes de usá-la, substitua cada underline por um espaço.

A rotina e seus resultados tem o seguinte formato: