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.
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.