Sistemas Lineares (projeto 1)

Neste projeto serão mostrados dois métodos de resolução de sistemas lineares Ax = b (método da eliminação de Gauss e usando a função solve do programa), e comparar os resultados obtidos para diferentes valores de N, onde N: ordem da matriz.

Os métodos serão aplicados em quatro itens com diferentes matrizes para serem resolvidas.

Para cada imagem exibida, é possível ampliá-la clicando na imagem.

i)

Usando o método solve: 

Código usado para gerar as matrizes A e b, e a solução:

Usando o método de eliminação de Gauss:

Código usado para gerar as matrizes A e b, e a solução:

O código usado para calcular o determinante de A e a razão entre o menor e o maior autovalores de A em módulo serão os mesmos para os dois métodos. Assim, o código é dado por:

Para N= 10 e 20 os códigos serão os mesmos, mudando apenas o valor da variável N, e a posição do maior autovalor da matriz A. Assim, será exibido apenas a solução do sistema.

Nas imagens seguintes serão exibidas os gráficos, e estarão identificados a ordem da matriz A (N), o valor do determinante da matriz A (detA).

Usando o método solve:

Usando o método gauss:

Sendo a variável "a" a razão entre o maior e o menor autovalores da matriz A em módulo, temos que:

N = 5: a = 476607.25024156051

N = 10: a = 16024663066292.234

N = 20: a = 1.103140868798481e+18

Conforme a ordem da matriz A aumenta, a razão entre o maior e o menor autovalores da matriz aumenta de forma muito acentuada, enquanto que o determinante da matriz diminui. A partir dos gráficos do vetor solução observa-se que o método solve apresenta uma oscilação muito grande, demonstrando uma maior instabilidade quando aumenta-se a ordem da matriz, enquanto que usando o método de eliminação de Gauss apresenta gráficos mais comportados.

Para os itens seguintes,o que será modificado serão apenas as matrizes A. O código para gerar o determinante de A e a razão entre o menor e o maior autovalores em módulo serão os mesmos. Assim, para os itens seguintes serão exibidos apenas o novo código de geração da matriz A, e a solução do sistema e interpretação.

ii)

Usando o método solve: 

Código usado para gerar as matrizes A e b, e a solução:

Usando o método de eliminação de Gauss:

Código usado para gerar as matrizes A e b, e a solução:

Nas imagens seguintes serão exibidas os gráficos, e estarão identificados a ordem da matriz A (N), o valor do determinante da matriz A (detA).

Usando o método solve:

Usando o método gauss:

Sendo a variável "a" a razão entre o maior e o menor autovalores da matriz A em módulo, temos que:

N = 5: a = 476607.25024156051

N = 10: a = 16024633473180.574

N = 20: a = 2.2595038042659722e+18

Mais uma vez, o aumento na ordem da matriz A influencia diretamente na razão entre o maior e o menor autovalores da matriz, e influencia inversamente no determinante da sua matriz. O método solve apresenta um aumento na oscilação do gráfico conforme a ordem da matriz A aumenta, apresentando instabilidade para a resolução do sistema. O método da eliminação de Gauss também apresenta oscilação, porém a instabilidade é menor, com um pico apenas no final do gráfico.

iii)

Usando o método solve: 

Código usado para gerar as matrizes A e b, e a solução:

Usando o método de eliminação de Gauss:

Código usado para gerar as matrizes A e b, e a solução:

Nas imagens seguintes serão exibidas os gráficos, e estarão identificados a ordem da matriz A (N), o valor do determinante da matriz A (detA).

Usando o método solve:

Usando o método gauss:

Sendo a variável "a" a razão entre o maior e o menor autovalores da matriz A em módulo, temos que:

N = 5: a = 9.969544242846343

N = 10: a = 7.0795177357028924

N = 20: a = 20.251387959149991

Neste item não é possível fazer uma conclusão em relação à mudança da ordem da matriz A, já que no enunciado foi exigido que a geração da matriz A fosse feita a partir de um comando que gerasse valores aleatórios entre 0 e 1. Tal problema está evidenciado nos valores obtidos no determinante e na razão entre o maior e menor autovalores da matriz A, onde os valores não demonstram nenhum padrão, aumentando e diminuindo. O comportamento dos gráficos também são inconclusivos, também por causa da geração aleatória dos valores da matriz A.

vi)

Usando o método solve: 

Código usado para gerar as matrizes A e b, e a solução:

Apesar do código para resolver o sistema estar presente, a implementação do código não é possível, assim o método solve não permite a obtenção da resolução de um sistema linear.

Usando o método gauss: 

Código usado para gerar as matrizes A e b, e a solução:

A solução do sistema linear usando o método de eliminação de Gauss é possível, e está representada nos seguintes gráficos:

Os gráficos de resolução do sistema mostram um certo padrão no seu comportamento, com pico apenas no fim do gráfico, e antes disso uma grande estabilidade. A razão entre o maior e menor autovalores não foi calculada pois existe o valor zero, o que faz com que o cálculo seja impossível (zero no denominador).

Autômatos Celulares Não Determinísticos (projeto 9)

Utilizando o Modelo de Ising, simular a evolução de sistemas com dimensões mínimas de 300 x 300 ao longo de 400 gerações considerando configurações iniciais aleatórias e beta inicial = 0.

Em cada um dos três casos a seguir, a variável beta deve ser incrementada em cada geração de:

• a) beta = 0.0001

• b) beta = 0.0008

• c) beta = 0.002

O que acontece em cada caso? Gere animações do processo de modo a ilustrar a dinâmica da evolução através de imagens do processo. Plote a configuração final obtida em cada caso.

Deve ser computada a entropia do sistema H(t) em cada geração e plotado um gráfico que expressa a variação de H(t) no tempo (gerações). 

H(t) = - [p0*log(p0) + p1*log(p1)] 

p0 = h0/N e p1 = h1/N, sendo h0 e h1 as componentes do histograma (números de 0's e 1's na configuração atual), e N a área do sistema (300*300)

Programa e sua execução para os 3 casos:

(clique na imagem para ampliá-la)

 Observação: infelizmente não foi possível transferir a animação para o blog, então estou colocando a imagem final da animação para os 3 casos.
 

Figura 1: primeira parte dos códigos

Figura 2: segunda e última parte dos códigos

Figura 3: Executando o programa para acréscimo de Beta = 0.0001

Figura 4: Gráfico de entropia das gerações para o acréscimo 0.0001

Figura 5: Último frame da animação das gerações com acréscimo de 0.0001 em Beta

Figura 6: Entropia das gerações, sob acréscimo de 0.0008

Figura 7: Último quadro da animação da simulação com acréscimo de 0.0008

Figura 8: Entropia sob acréscimo de 0.002

Figura 9: Final da animação da simulação sob acréscimo em Beta de 0.002 a cada geração

Podemos observar que conforme o valor de Beta aumenta, certa "organização" vai aparecendo e o "caos" torna-se pequeno, as células agrupam-se com maior facilidade, porque beta representa a "importância" que se dá à vizinhança de uma célula para calcular seu futuro estado. Por se tratar de um modelo não determinístico, toda vez que o programa for executado novamente (inclusive com mesmo valor de acréscimo) uma animação e gráfico diferentes serão gerados.

 Código para ser copiado:

#importando bibliotecas necessárias ao programa
import numpy as np
import math as math
import matplotlib.pyplot as plt
import matplotlib.animation as animation

MAX = 400 #número de gerações
SIZE = 300 #dimensões do sistema
beta = 0 #beta inicial
#pedir ao usuário o valor de quanto beta cresce em cada geração
crescimento_beta = float(input('Quanto Beta deve crescer a cada geração?( 0.0001; 0.0008; 0.002 ) '))

geracao = np.random.randint(0, 2, (SIZE,SIZE)) #gerando matriz inicial com 0's ou 1's aleatórios
nova_geracao = np.zeros((SIZE, SIZE)) #matriz onde será armazenada a proxima geração, preenchida com zeros
matriz_evolucao = np.zeros((SIZE, SIZE, MAX)) #matriz 3D onde será gravado todo o processo
entropia = [] #lista vazia onde será armazenado o valor da entropia em cada geração

# laço principal, criará todas as gerações
for k in range(MAX):

      print("Processando geração %d..." %k)#mostra a geração em processo
      matriz_evolucao[:,:,k] = geracao #k-ésima geração será feita
     
      for i in range(1, SIZE-1): #desconsidera as bordas da matriz, p/ que exista vizinhança das células nas extremidades
            for j in range(1, SIZE-1): #confere as 8 células vizinhas a uma célula central

                  # seleciona janela 3X3 (vizinhança: 1 central com 8 ao seu redor)
                  bloco = geracao[i-1:i+2, j-1:j+2]

                  # soma elementos da janela para verificar o nº de 0s e 1s ao redor da célula
                  numero_ums = sum(sum(bloco)) - geracao[i,j]#subtrai o elemento central(queremos apenas a vizinhança dele)
                  numero_zeros = 8 - numero_ums #obtendo o nº de zeros nessa vizinhança

                  #probabilidade que usa a vizinhança (ou não se beta=0) para definir se a célula em questão se tornará 0 ou 1
                  p = np.exp(beta*numero_zeros)/(np.exp(beta*numero_zeros)+np.exp(beta*numero_ums))

                  # X é nº aleatório entre 0 e 1
                  x = np.random.random()

                  # aplicar decisão, com auxilio da probabilidade e do X
                  if x<=p:
                        nova_geracao[i,j] = 0
                  else:
                        nova_geracao[i,j] = 1
      geracao = nova_geracao.copy() #após calculada a nova geração, ela passa a ser a atual
      beta = beta + crescimento_beta #atualização do beta

      #cálculo da entropia do sistema
      p1 = (np.sum(nova_geracao[:,:]))/(SIZE*SIZE) #proporção de células em estado 1, na geração atual, ou seja,
                                                    #soma todas os valores da matriz e divide pelo seu tamanho
      p0 = 1 - p1 #proporção de células em estado 0, na atual geração
      h = -( p0*(np.log10(p0)) + p1*(np.log10(p1)) ) #valor da entropia nessa geração
      entropia.append(h) #acrescentando a entropia atual na lista das entropias
     
     

# gerando animação
fig = plt.figure(1)
lista = [] #lista vazia que conterá os frames da animação
for i in range(MAX):
      im = plt.imshow(matriz_evolucao[:,:,i], cmap='gray')#cada frame da animação representa uma geração
      lista.append([im]) #acrescenta o novo frame a lista
anime = animation.ArtistAnimation(fig, lista, interval=50, blit=True, repeat_delay=1000)
plt.show()

    
#gerando o gráfico dos valores da entropia em cada geração
plt.plot(entropia)
plt.show()
  
#fim

Simulação de trending topics (projeto 10)

Neste projeto o objetivo consiste em, dado um arquivo com inúmeras citações famosas em inglês, encontrar a distribuição das palavras do texto, ou seja, a frequência com que cada uma ocorre no texto. Não devem ser considerados artigos, verbo to be, preposições, números, símbolos, etc. (a, an, in, on, if, it, is, I, you, we, the, but), (0, 1, 2,..., 9), (!, *, &, ", '). De forma geral, palavras com menos de 3 letras devem ser descartadas.

Deve ser montado um ranking, com as top 10 palavras mais usadas nas citações (de modo a simular os trending topics do Twitter).

Arquivos podem ser encontrados em:

http://www.textfiles.com/humor/TAGLINES/quotes-1.txt

http://www.textfiles.com/humor/TAGLINES/quotes-2.txt

http://www.textfiles.com/humor/TAGLINES/quotes2.txt

http://www.textfiles.com/humor/TAGLINES/quotes.frt

http://www.textfiles.com/humor/TAGLINES/bonehead.txt

http://www.textfiles.com/humor/TAGLINES/quotes.jok

http://www.textfiles.com/humor/TAGLINES/ansafone.txt

Programa e sua execução:

(clique na imagem para ampliá-la)

 Código para ser copiado:

##### CRIANDO UM RANKING DE PALAVRAS QUE MAIS APARECEM EM UM TEXTO #####


#OBS: o arquivo .txt/.jok/.frt/"etc" que deseja utilizar, deve estar no mesmo diretorio que este .py

import string #importando biblioteca que contém caracteres especiais não desejados

#armazenando o aquivo de texto
arquivo = open(input('Informe o nome do arquivo (com sua extenção): ')) #usuário digita o nome do arquivo
texto = arquivo.read() #deixando o texto 'pronto para leitura' numa string
texto = texto.lower() #deixando todas as letras minúsculas

for c in string.punctuation: #eliminando caracteres indesejáveis (!?:;., etc)
    texto = texto.replace(c,'')

texto = texto.split() #quebra a string numa lista de palavras

#criando um dicionário a partir da lista
d={}
for p in texto:
    if p not in d:
        d[p] = 1
    else:
        d[p] += 1

arquivo.close() #fechando o arquivo para não ocupar memória do computador

#criando ranking das palavras que mais aparecem
lista_palavras = sorted(d, key=d.get, reverse = True) #ordenar a lista de forma decrescente
i = 0 #contador para mostrar apenas 10 palavras
for palavra in lista_palavras:
    if len(palavra)>3 and i<10: #selecionando palavras com mais de 3 letras, mostrar até 10 palavras
        print('%s : %d' %(palavra, d[palavra])) #mostrando a palavra e seu respectivo nº de aparições
        i += 1 #atualização do contador de palavras mostradas
       
#fim

Aproximação númerica por séries (projeto 3)

Sabe-se que diversas funções matemáticas podem ser aproximadas através de séries infinitas. Utilizando as definições a seguir,desenvolva um (ou vários) programa(s) em Python para calcular o valor das seguintes funções num ponto x. Tanto x quanto n (número de termos) devem ser conhecidos.

• a) exp(x) = 1 + x + x^2/2! + x^3/3! + x^4/4! + ...

• b) cos(x) = 1 - x^2/2! + x^4/4! - x^6/6! + x^8/8! - ...

• c) sen(x) = x - x^3/3! + x^5/5! - x^7/7! + x^9/9! - ...

• d) integral de 0 a x exp{-t^2} dt = x - x^3/3*1! + x^5/5*2! - x^7/7*3! + x^9/9*4! - ...

Sabendo disso compute a integral de

  0 a 1 de 1/sqrt(pi) exp{-t^2} dt

o que significa isso em termos de probabilidade (se t é uma VA normal)? 

RESOLUÇÃO

CÓDIGOS INSERIDOS NO PYTHON + COMENTÁRIO E EXECUÇÃO

(os códigos podem ser copiados no final dessa postagem)

 

Comandos ao lado esquerda das imagens e execução ao direito. Clique na imagem para ampliá-la.

 

Figura 1: Fatorial e exponencial

 

Figura 2: Cosseno

 

Figura 3: Seno

Figura 4: Integral da função

Figura 5: Cálculo da integral definida

 

COMANDOS PARA SEREM COPIADOS:

 

import math as math #importando a biblioteca necessária para o programa

#criando a função para calcular fatorial
def fat(n):
    if n <= 1:
        return 1
    return fat(n-1) * n

#################################################################   
#criando a função chamada "expoente" para calcular e^x

def expoente():

    #pedindo ao usuário o valor de x e o nº de termos para a série do cálculo
    x = float(input('informe o valor do numero real X: '))
    n = int(input('informe o número inteiro de termos usados no cálculo: '))

    #definindo como 0 o valor inicial da série, para que a soma da série
    #nao entre num loop infinito durante seu cálculo
    s = 0
   
    if (x==0): #caso seja pedido e^0
        print('O resultado é 1')
    else:

        #calculando e^x
        for i in range(n):
            s = x**i/ fat(i) + s
            #aqui é calculado cada termo da série e somado ao anterior,
            #sucessivamente até atingir os N termos
           
    #expressando o resultado final, com algarismos depois da vírgula
    print('O resultado é %5f' %s)
    return

#################################################################
# criando uma função que calcula cos(x)

def cosseno():

    x = float(input('informe o valor do numero real X: '))
    n = int(input('informe o número inteiro de termos usados no cálculo: '))

    s = 0
   
    #calculando o valor de cos(x)
    for i in range(n):
        s = ( (x**(2*i))*((-1)**(i)) )/(fat(i*2)) + s
               
    print('O resultado é %5f' %s)
    return

##################################################################
# criando uma função para calcular sen(x)

def seno():

    x = float(input('informe o valor do numero real X: '))
    n = int(input('informe o número inteiro de termos usados no cálculo: '))

    s = 0
   
    #calculando o valor de cos(x)
    for i in range(n):
        s = ( (x**(2*i+1))*((-1)**(i)) )/(fat(2*i+1)) + s
               
    print('O resultado é %5f' %s)
    return
   
###################################################################
# criando uma função que retorna a integral de 0 a X de exp{-t^2}dt

def integral(x,n):#aqui o usuário deverá informar o X e N logo que
                  #chamar a função
       
    s = 0
   
    #calculando o valor da integral
    for i in range(n):
        s = ( (x**(2*i+1))*((-1)**(i)) )/((2*i+1)*fat(i)) + s

    #aqui a função retorna diretamente o valor do resultado final           
    return(s)

###################################################################
#computando a integral de 0 a 1 de (1/sqrt(pi))*exp{-t^2} dt
                                               y = (1/math.sqrt(math.pi))*integral(1,500)



#em termos de probabilidade, se t é um variável aleatória da distribuição Normal
#esta possui média=0 e variância=1/2, e então
#o valor y significa a probabilidade de um valor estar entre 0 e 1

#fim

"Quem ensina aprende duas vezes." (Joseph Joubert)