====== Física Computacional - 2017.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 =====

  * [[http://www.ebah.com.br/content/ABAAAgUS4AF/gnuplot-manual-simplificado-iniciantes|Guia para iniciantes]]

  * [[http://www.notasemcfd.com/2010/03/gnuplot-tipos-de-pontos-e-linhas.html|Estilo de pontos e linhas]]

  * {{:nuno:gnuplot-ajuste.pdf|Ajustes numéricos}}

  * [[http://complex.if.uff.br/marcio/semana2009|Ferramentas de Computação em gnuplot - página do prof. Marcio Argollo de Menezes]]

  * {{:nuno:gnuplot_aula.pdf|Aula de gnuplot}}

===== Referências básicas de Programação =====

  * Manuais de Fortan 77:
  - [[http://www.idris.fr/data/cours/lang/fortran/f90/F77.html|Fortan 77 for beginners]]
  - [[http://web.stanford.edu/class/me200c/tutorial_77/|Fortran 77 tutorial]]

  * Manuais de Fortran 90:
  - [[http://www.cenapad.unicamp.br/servicos/treinamentos/apostilas/apostila_fortran90.pdf|PDF - Curso da Unicamp]]

  * Manuais de C/C++:
  - [[http://www.cenapad.unicamp.br/servicos/treinamentos/apostilas/apostila_C.pdf|Curso de C da Unicamp]]
  - [[http://www.mtm.ufsc.br/~azeredo/cursoC/c.html|Curso de C da UFMG]]
  - [[http://imada.sdu.dk/~svalle/courses/dm14-2005/mirror/c/|C Programming]]



===== Listas de Exercícios =====

<note important>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</note>

**Construção de Histogramas**


  * Assim como discutimos em sala, baixe {{:nuno:tempos_de_relaxacao.txt|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? 

  * Construa o mapa de bifurcações para 0.0 < λ < 4.0





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







===== Notas =====

  * {{:nuno:notas_fiscomp_2017_1.xls|Notas}}



















===== 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)