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.
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
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
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@mail.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é 10/04/17)
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. Se você desejar, 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é 28/04/17)
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é 22/05/17)
Resolva numericamente a equação do pêndulo simples através do mapa iterativo gerado pelo uso do método de Euler, com os dados 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 1o 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 1, 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é 23/06/17)
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?
Difusão de Calor e Caminhadas Aleatórias (entregar até 18/07/17)
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 que essa constante vale √(2D). Verifique também essa relação dessa constante de proporcionalidade para D=0,1 e 0,3.
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é 18/07/17)
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; (iv) a=16807, m=2147483647. Calcule o período destes geradores, considerando como semente r=27. Observe que você deve gerar uma quantidade de números aleatórios suficiente para cada caso, dado que o período destes geradores vai variar bastante.
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 (iv) 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 (iv), 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.10 < y < 3.18.
Método de Monte Carlo e o modelo de Ising (entregar até 18/07/17)
Simule o modelo de Ising na rede quadrada L x L, com condições de contorno periódicas, considerando L=128. Inicie com uma configuração aleatória de spins, com 50% de spins +1 e 50% de spins -1. Plote curvas da magnetização por spin em função do tempo. Você deve ver a magnetização flutuar e variar no tempo, até atingir o equilíbrio. Considere diferentes temperaturas, algumas abaixo de Tc=2.27 (por exemplo T=1.0, 1.5, 1.8, 2.0 e 2.2) e outras acima de Tc (por exemplo 2.5, 2.7, 3.0 e 3.5). Não é necessário fazer médias sobre simulações diferentes.
Simule o modelo de Ising na rede quadrada L x L, com condições de contorno periódicas, considerando L=16, 24, 32, 64 e 128. Plote as curvas das quantidades de interesse no equilíbrio (magnetização por spin média, energia média, susceptibilidade e calor específico) em função da temperatura para os 5 tamanhos de rede. Lembre de descartar alguns passos de Monte Carlo para garantir que o sistema esteja no equilíbrio. Faça médias temporais sobre uma quantidade razoável de passos de Monte Carlo. Não é necessário fazer médias sobre simulações diferentes. Verifique e discuta os efeitos de tamanho finito que falamos em sala, assim como o comportamento esperado para as quantidades de interesse.
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)
