Essa é uma revisão anterior do documento!


Física Computacional - 2015.1

Neste curso estudaremos alguns métodos computacionais voltados para a solução de problemas de Física. Sempre que possível iremos tratar sistemas que possuem solução analítica, de forma a testarmos a eficiência dos métodos numéricos. Os tópicos vão desde a análise de sistemas físicos básicos descritos por equações diferenciais ordinárias ou parciais, até problemas de difusão, caos e modelos estatísticos.

Estrutura do Curso

  • Introdução à Modelagem Computacional em Física, Ajustes de Resultados Numéricos, Manipulação de Dados, Construção de Histogramas
  • Integração de Equações Diferenciais e Aplicações: Leis de Newton, Resfriamento de Fluidos, Decaimento Nuclear, Queda Livre, Pêndulo Simples e Amortecido, Difusão de Calor, Dinâmicas Populacionais, etc
  • Introdução ao Caos, Mapa Logístico
  • Geração de Números Pseudo-aleatórios
  • Estatística e Simulação de Caminhadas Aleatórias, Caminhadas Aleatórias Modificadas (persistentes, auto-excludentes, restritas, etc), Tempo de Primeira Passagem, Aplicações
  • Método de Monte Carlo, Modelo de Ising e o Algoritmo de Metropolis, Análise de Tamanho Finito, Expoentes Críticos
  • Tópico Avançado Opcional: Autômatos Celulares e Fractais (introdução)

Tutoriais de Gnuplot

Referências básicas de Programação

Listas de Exercícios

Como falei na 1a aula, faça um pdf com todos os seus resultados (gráficos, valores dos parâmetros de ajuste, discussão, etc) e me envie por email: nuno@if.uff.br

Construção de Histogramas

  • Assim como discutimos em sala, baixe este arquivo, que contém tempos de relaxação de um modelo estatístico. Utilize inicialmente um intervalo entre pontos delta=100, e calcule quantos valores dos mencionados tempos ocorrem e plote o histograma no gnuplot ou no seu programa preferido de gráficos. Se for no gnuplot, vc pode gerar um arquivo do tipo PS ou EPS e me enviar por email. Depois, experimente outros valores de delta como delta=10, 50 e 200, e faça comparações entre os resultados. ESTE EXERCÍCIO NÃO FAZ PARTE DA AVALIAÇÃO!

Decaimento Nuclear (entregar até 27/03/15)

  • Utilize o mapa N(t+1) = N(t)*[1-α] e construa um gráfico do número de núcleos ainda radioativos em função do tempo. Considere como condição inicial N(0)=1000 e utilize o valor de α do rubídio 82, que é 0,00924 1/s. Para simplificar, imprima no arquivo de dados apenas os valores obtidos a cada 10s, e comece com um tempo máximo de 100 s.
  • Construa o mesmo gráfico anterior, com os mesmos dados, adotando desta vez uma escala logarítmica (com log na base neperiana) no eixo y. Isto pode ser feito de 2 formas: (i) calculado no seu programa o ln dos valores obtidos para N(t) e plotando numa 3a coluna do arquivo de dados, ou (ii) utilizando os comandos do gnuplot para calcular o ln a partir dos dados originais de N(t). Seu gráfico tem o formato linear desta vez?
  • Meça o coeficiente angular da reta obtida no problema anterior e compare com a constante α do rubídio 82 (0,00924 1/s). Neste caso, mostre que a forma analítica que descreve o problema é N(t)=N(0)*exp(-αt).
  • Considere diferentes valores do tempo máximo de aplicação do mapa, como 100, 500 e 1000, e compare as estimativas obtidas para o coeficiente angular da reta.

Integração Numérica de Equações Diferenciais (entregar até 17/04/15)

  • Calcule a derivada numérica da função f(x)=[sen(x^2) exp(x/3)]/[√(x^2+4)] no ponto x=3 usando o método de Euler simples com passo de discretização entre Δt=10^{-8} e Δt=1 (ou seja, 10^{-8}, 10^{-7}, 10^{-6},…, 1) e faça um gráfico mostrando a diferença relativa entre a derivada numérica e o valor exato da derivada. Entre quais valores de Δt a aproximação é aceitável?
  • Integre a equação diferencial dy/dt - y = -(1/2)exp(t/2)sen(5t) + 5exp(t/2)cos(5t) com a condição inicial y(0)=1 com o método de Euler simples para diferentes valores de passo Δt=0.1, Δt=0.05, Δt=0.01, Δt=0.005 e Δt=0.001 e compare cada solução em t=1,2,3,4 e 5 com a solução exata y(t)=exp(t) + exp(t/2)sen(5t).
  • Faça um gráfico com a solução para Δt=0.05 entre t=0 e t=5 e a solução exata.
  • Repita o item anterior utilizando o método de Runge-Kutta de 2a ordem e de 4a ordem, ou seja, faça um gráfico com a solução para Δt=0.05 entre t=0 e t=5 (para os 2 métodos de Runge-Kutta) e a solução exata.

Pêndulo (entregar até 18/05/15)

  • Considere um pêndulo simples de comprimento L=2,00 m cuja velocidade inicial é v_o=1,00 m/s na posição vertical. O valor da aceleração local da gravidade é g=9,8 m/s^2. Dados estes valores, justifica-se ou não a aproximação de pequenas oscilações neste caso?
  • Resolva numericamente a equação do pêndulo simples através do mapa iterativo gerado pelo uso do método de Euler, com os mesmos dados anteriores (L=2,00 m e v_o=1,00 m/s) e plote θ em função do tempo t. Escolha um intervalo Δt adequado à precisão compatível com a dos dados apresentados. Além de θ_o=0, você necessitará conhecer a priori o valor de θ_1, a partir da condição inicial dθ/dt=v_o/L. Basta notar que a derivada de θ(t) no instante t=0 pode ser escrita como dθ(t=0)/dt = θ_1/Δt.
  • Mostre que o período do pêndulo vale T = 2π√(L/g) dentro da aproximação de pequenas oscilações, e verifique este resultado no caso do gráfico do problema anterior. Você pode fazer isso plotando os dados do item anterior e uma reta vertical com o resultado analítico T = 2π√(L/g).
  • Compare os resultados obtidos no 2o item com a aproximação de pequenas oscilações, ou seja, plote os resultados analítico e numérico juntos no mesmo gráfico.
  • Repita o problema do item 2, desta vez usando uma velocidade inicial 5 vezes menor, v_o=0,200 m/s. Compare o período T e a amplitude θ_A obtidos com os valores previstos pela aproximação de pequenas oscilações. O acordo é bom? por quê?
  • Repita novamente, desta vez usando uma velocidade inicial 5 vezes maior, v_o=5,00 m/s. Compare o período T e a amplitude θ_A obtidos com os valores previstos pela aproximação de pequenas oscilações. O acordo é bom? por quê?
  • Repita o problema com diferentes valores de v_o, todos menores que o valor crítico v_c=2√(Lg), e meça o período T e a amplitude θ_A em cada caso. Construa um gráfico de T em função de θ_A.
  • Se construíssemos um pêndulo real para experimentos equivalentes aos que tratamos numericamente, não poderíamos admitir que o atrito da resistência do ar é desprezível. Modelando este atrito como sendo proporcional à velocidade, a EDO a ser resolvida é d^2 θ/dt^2 + bdθ/dt + (g/L)senθ=0, onde b é uma constante que representa a intensidade do atrito. Este é o chamado Pêndulo Amortecido. Resolva esta equação numericamente usando o método de Euler, novamente com L=2,00 m e v_o=1,00 m/s, considerando um atrito relativamente pequeno, b=1,00 1/s. Você poderá verificar que o movimento ainda é oscilatório, mas não é mais periódico. Portanto, não há mais período e sim o tempo decorrido durante 1 oscilação completa, que poder ser medido. Meça-o e compare-o com o período obtido sem atrito.
  • Repita o problema anterior para diferentes valores de b, e analise os resultados. O movimento é sempre oscilatório? Adote diferentes valores de velocidade inicial, inclusive maiores que v_c.

Mapa Logístico e Caos (entregar até o fim do período letivo, 12/07/15)

  • Considere o mapa logístico x(t+1)=λx(t)[1-x(t)]. Começando com uma condição inicial x(0)=0,0001, mostre que para λ=1,01 o sistema se estabiliza no ponto fixo x*=0,00990 depois de cerca de 20min. Neste escala de tempo, lembre-se que cada iteração do mapa corresponde a 1s.
  • Construa o gráfico de x' contra x para λ=2, considerando o mapa logístico na forma x'=λx(1-x). Comece a escada a partir de x=0,01. Plote a escada, a parábola λx(1-x) e a reta x'=x no mesmo gráfico, e estime o valor do atrator x*. Este valor está de acordo com o resultado analítico? Repita toda a análise para λ=0,9 e 1,5, escolhendo valores adequados do valor inicial de x para cada caso.
  • Construa o gráfico de x' contra x para λ=3,1. Comece a escada a partir de x=0,2. Plote a escada, a parábola λx(1-x) e a reta x'=x no mesmo gráfico. Você consegue identificar o valor do atrator x*? Você deverá verificar que há mais de 1 atrator neste caso. Determine os valores destes atratores usando o método discutido em sala. Repita para λ=3,47 e λ=3,55.
  • Itere umas dez mil vezes o mapa logístico para λ=0,9. Comece em t=0 com um valor inicial arbitrário, x=0,1, e armazene num arquivo os dados x e t. Plote este gráfico e observe que você obterá o decaimento temporal de x(t). Que função é esta que descreve o decaimento no tempo?
  • Itere umas dez mil vezes o mapa logístico na situação crítica, λ=1. Comece em t=0 com um valor inicial arbitrário, x=0,1, e armazene num arquivo os dados x e t. Plote este gráfico e observe que você obterá o decaimento temporal de x(t). Que função é esta que descreve o decaimento no tempo?
  • Construa o mapa de bifurcações para 0.0 < λ < 4.0

Difusão de Calor e Caminhadas Aleatórias (entregar até o fim do período letivo, 12/07/15)

  • Uma barra metálica com temperatura inicialmente uniforme é submetida a um choque térmico localizado no seu centro. Escolha a barra entre x=0 e x=+L, ou entre x=-L/2 e x=+L/2, a seu critério. Considere a condição inicial u(x,0) nula em todos os pontos da barra, exceto no centro, onde temos u=1. Considere Δx=Δt=1, em unidades arbitrárias nas quais o coeficiente de difusão vale D=0,2. Determine numericamente através do método de Euler a distribuição de temperaturas ao longo da barra, ou seja, a função u_x,t para t=1,2,4,8,16,32,64,128,256,512 e 1024.
  • Com os dados do problema anterior, calcule a dispersão Δ(t) nos mesmos instantes anteriores e mostre que Δ(t) ~ t^{1/2}. Quanto vale a constante de proporcionalidade da relação anterior? Dê o resultado numérico e verifique também o valor dessa constante para D=0,1 e 0,3. Você deverá estimar a relação entre essa constante de proporcionalidade e o coeficiente de difusão D baseado nestes 3 valores.
  • Faça a simulação do problema do caminhante aleatório em 1 dimensão utilizando p=q=1/2. Cada caminhante deve dar 100.000 passos, e faça médias sobre 10.000 caminhantes independentes. Plote num gráfico os os valores médios <x> e <x^2> em função do tempo, e em outro gráfico a dispersão Δ(t) = <x^2> - <x>^2. Como estes valores crescem com o tempo?
  • Estime também a distribuição final de posições dos 10.000 caminhantes, após os 100.000 passos. Coloque na escala log-linear, e você deverá observar uma parábola. Que distribuição é essa? ela mudará seu formato se você utilizar valores diferentes de p?

Geradores de Números Aleatórios (entregar até o fim do período letivo, 12/07/15)

  • Vamos considerar o gerador linear congruencial r = (a*r) mod m. Considere as escolhas: (i) a=85, m=256; (ii) a=899, m=32768; (iii) a=16807, m=4294967295. Calcule o período destes geradores, considerando como semente r=27.
  • Construa histogramas dos números aleatórios gerados pelos geradores (ii) e (iii) do item anterior, para diferentes valores do total de números aleatórios gerados (10^{2}, 10^{3}, 10^{4}, 10^{5}, 10^{6} e 10^{7}).
  • Dada um probabilidade p (faça separadamente para p=0.3, 0.5 e 0.9), gere N números aleatórios usando o gerador (iii) acima, normalizando estes números para ficarem no intervalo [0,1]. Conte o número n deles que se situam entre 0 e p. A fração n/N deve reproduzir o valor de p fornecido. Considerar N=10^{3}, 10^{4}, 10^{5}, 10^{6}, 10^{7} e 10^{8}. Para verificar o resultado, plote a fração n/N contra N, e juntamente plote uma reta horizontal com o valor de p. Utilize a escala log para o eixo x, ficará mais fácil de visualizar.
  • Gerar N pares (x,y) de números aleatórios com o gerador (iii), todos os números no intervalo [0,1]. Todos estes pontos estarão localizados no interior de um quadrado de lado unitário. Os pontos que satisfazem à condição x^{2} + y^{2} < 1 também estarão localizados no interior de 1/4 de uma circunferência de raio unitário. Conte este número n de pontos dentro do círculo e obtenha a estimativa para π através de π=4n/N, conforme discutido em sala. Faça para diferentes valores de N, tais como N=10^{3}, 10^{4}, 10^{5}, 10^{6}, 10^{7}, 10^{8} e 10^{9}, e plote o valor obtido π=4n/N contra N juntamente com uma reta horizontal com o valor π=3,141592. Utilize a escala log para o eixo x, ficará mais fácil de visualizar, e plote o eixo y no intervalo 3.12 < y < 3.18.

Método de Monte Carlo (entregar até o fim do período letivo, 12/07/15)

Monitoria

  • O aluno de doutorado Allan Rodrigues Vieira é o monitor do curso. Se precisar tirar dúvidas, entre em contato com ele por email: allanrv@if.uff.br ou procure por ele no gabinete A3-25
  • O horário de atendimento será: Terças feiras de 14h às 16h

Bibliografia

  • Paulo Murilo Castro de Oliveira, Suzana M. Moss de Oliveira, Física em Computadores, Coleção Tópicos de Física do CBPF, Editora Livraria da Física
  • Harvey Gould, Jan Tobochnik, An Introduction do Computer Simulation Methods: Applications to Physical Systems, 2nd Edition, Addison-Wesley Publishing Company
  • Tânia Tomé, Mário José de Oliveira, Dinâmica Estocástica e Irreversibilidade, Editora da Universidade de São Paulo (edUSP)
nuno/fiscomp_2015_1.1434379073.txt.gz · Última modificação: 2015/06/15 11:37 por nuno
CC Attribution-Share Alike 3.0 Unported
www.chimeric.de Valid CSS Driven by DokuWiki do yourself a favour and use a real browser - get firefox!! Recent changes RSS feed Valid XHTML 1.0