Física Computacional - FSC-5705

só um divisor

Prova de Física Computacional,

I semestre de 2013

Antes de iniciar a prova é necessário que o estudante crie uma pasta denominada prova_I_matricula_do_aluno dentro desta pasta serão feitos cada um dos programas correspondente aos problemas da prova. O nome dos programas serão problema1.py, problema2.py, etc. Ao concluir a prova o estudante tem de criar um arquivo comprimido com tudo, isto é criar um arquivo do tipo .tar.gz. Este arquivo devera ser enviado por email para o endereço evy.fisica@gmail.com e o assunto do email deve ser prova 1 - matricula do aluno. Depois disso deverá chamar ao professor a fim de copiar para um pendrive o mesmo arquivo.

No mínimo escolha tantos problemas como necessários para completar 10 pts. A nota pode ser maior do que 10.

Exercício 01 (3 pts)

Um diodo é um dispositivo não linear que permite o passo de uma corrente elétrica quando uma diferencia de potencial positiva é colocada entre seus terminais caso contrario a corrente elétrica é quase nula. Por esta razão se utilizam diodos como retificadores. A equação da corrente que passa por um diodo ideal está dada por \[ I = I_s e^{(qV)/(nk_BT)} \nonumber \] onde

escreva um programa para calcular a corrente num diodo de silício onde se varia a diferencia de potencial $V$ desde $-0.20$ até $0.55$ volts, com intervalos de $0.05$ volts. Faça o gráfico de $I$ em função de $V$.

Solução

Este problema é direto, só deve prestar atenção à condição $I=0$ se $V<=0$

 
      #/usr/bin/env python
      #-*- coding: utf-8 -*-

      import numpy as np
      import matplotlib.pyplot as plt

      Is = 0.01e-6
      q  = 1.6e-19
      n  = 2.0
      K  = 1.38e-23
      T  = 300.0

      V  = np.arange(-0.20, 0.55, 0.05)
      I  = np.zeros_like(V)

      for i in range(len(V)):
        if (V[i] > 0.0):
          I[i] = Is * np.exp(q*V[i] / (n*K*T))

      plt.title(u"Exercício 01", fontsize=30)
      plt.xlabel("V", fontsize=25)
      plt.ylabel("I", fontsize=25)
      plt.tight_layout()
      plt.plot(V, I, ':', marker="o")
      plt.savefig("exe010.png")
      plt.show()
      

Exercício 02 (3 pts)

Um corpo com velocidade $\vec{v}_1$ e massa $m_1$ colide com um corpo com velocidade $\vec{v}_2$, massa $m_2$ no ponto $p$. Crie um programa que pergunte esses dados iniciais e calcule as velocidades finais dos corpos. O programa deve perguntar ao usuário se a colisão é perfeitamente elástica o totalmente inelástica e realizar os cálculos em base a essas opções. Utilize visual Python para representar os vetores antes e depois da colisão.

Solução

Neste problema era necessário relembrar de conceitos básicos de Física I. O exercício apresentava duas situação:

A solução matemática é fácil de obter, mas digitando no seu buscador preferido você poderia utilizar a solução dada pelo professor Canzian a qual vou reproduzir a continuação

Colisão elástica

Seja $m_1$, $\vec{v}_{1i}$ e $\vec{v}_{1f}$ a massa, velocidade antes da colisão e velocidade depois da colisão do corpo $1$, respectivamente; e $m_2$, $\vec{v}_{2i}$ e $\vec{v}_{2f}$ a massa, velocidade antes da colisão e velocidade depois da colisão do corpo $2$, respectivamente, as lei de conservação nos dizem que: \[ \begin{align} m_1\vec{v}_{1i} + m_2\vec{v}_{2i} &= m_1\vec{v}_{1f} + m_2\vec{v}_{2f} \label{p1}\\ \frac{1}{2}m_1\vec{v}_{1i}^2 + \frac{1}{2}m_2\vec{v}_{2i}^2 &= \frac{1}{2}m_1\vec{v}_{1f}^2 + \frac{1}{2}m_2\vec{v}_{2f}^2\label{e1} \end{align} \] cancelando os $1/2$ da expressão $\ref{e1}$ e agrupando em ambas equações os termos com massa iguais: \[ \begin{align} m_1\left( \vec{v}_{1i} - \vec{v}_{1f} \right) &= m_2\left( \vec{v}_{2f} - \vec{v}_{2i} \right) \label{p2}\\ m_1\left( \vec{v}^2_{1i} - \vec{v}^2_{1f} \right) &= m_2\left( \vec{v}^2_{2f} - \vec{v}^2_{2i} \right) \label{e2} \end{align} \] observe que podemos reescrever a equação $\ref{e2}$ da seguinte forma \[ m_1\left( \vec{v}_{1i} - \vec{v}_{1f} \right)\left( \vec{v}_{1i} + \vec{v}_{1f} \right) = m_2\left( \vec{v}_{2f} - \vec{v}_{2i} \right)\left( \vec{v}_{2f} + \vec{v}_{2i} \right) \nonumber \] dividindo esta última equação pela equação $\ref{p2}$ \[ \frac{m_1\bcancel{\left( \vec{v}_{1i} - \vec{v}_{1f} \right)}\left( \vec{v}_{1i} + \vec{v}_{1f} \right)} {\bcancel{m_1\left( \vec{v}_{1i} - \vec{v}_{1f} \right)}}= \frac{\bcancel{m_2\left( \vec{v}_{2f} - \vec{v}_{2i} \right)}\left( \vec{v}_{2f} + \vec{v}_{2i} \right)} {\bcancel{m_2\left( \vec{v}_{2f} - \vec{v}_{2i} \right)}} \nonumber \] de onde \[ \begin{align} \vec{v}_{1i} + \vec{v}_{1f} =& \vec{v}_{2f} + \vec{v}_{2i} \nonumber\\ \vec{v}_{1f} =& \vec{v}_{2f} + \vec{v}_{2i} - \vec{v}_{1i} \label{pe1} \end{align} \] colocando em $\ref{p2}$ \[ \begin{align} m_1\left( \vec{v}_{1i} - \vec{v}_{2f} - \vec{v}_{2i} + \vec{v}_{1i} \right) &= m_2\left( \vec{v}_{2f} - \vec{v}_{2i} \right) \nonumber\\ m_1\left( 2\vec{v}_{1i} - \vec{v}_{2f} - \vec{v}_{2i} \right) &= m_2\left( \vec{v}_{2f} - \vec{v}_{2i} \right) \nonumber\\ \vec{v}_{2f} \left( m_2 + m_1 \right) &= \vec{v}_{2i} \left( m_2 - m_1 \right) + 2m_1\vec{v}_{1i}\nonumber\\ \vec{v}_{2f} &= \frac{m_2 - m_1}{m_2 + m_1} \vec{v}_{2i} + \frac{2m_1}{m_2 + m_1} \vec{v}_{1i} \end{align} \] substituindo em $\ref{pe1}$: \[ \begin{align} \vec{v}_{1f} =& \frac{m_2 - m_1}{m_2 + m_1} \vec{v}_{2i} + \frac{2m_1}{m_2 + m_1} \vec{v}_{1i} + \vec{v}_{2i} - \vec{v}_{1i}\nonumber\\ \vec{v}_{1f} =& \frac{m_1 - m_2}{m_2 + m_1} \vec{v}_{1i} + \frac{2m_2}{m_2 + m_1} \vec{v}_{2i} \end{align} \]

Com base no resultado anterior, o script que resolve este problema é

 
      #!/usr/bin/env python
      #!/usr/bin/env python
      #-*- coding: utf-8 -*-

      import numpy as np
      import visual as vis


      m1  = input('digite a massa da partícula 1: ')
      m2  = input('digite a massa da partícula 2: ')
      while(True):
        v1i = input('digite a velocidade da partícula 1: ')
        v2i = input('digite a velocidade da partícula 2: ')
        if (np.sign(v1i) == np.sign(v2i)):
          print "Para ter colisão as velocidades devem ter sentidos opostos"
        else:
          break

      while (True):
        tipo = input('digite o tipo de colisão: 1)Elástica, 2) inelástica: ')
        if (tipo == 1 or tipo ==2):
          break

      mtot = m1 + m2
      if (tipo == 1):
        v1f = 2.0*m1*v1i / mtot + v2i * (m2 - m1) / mtot
        v2f = 2.0*m2*v2i / mtot + v1i * (m1 - m2) / mtot
      else:
        v1f = (m1*v1i + m2*v2i) / mtot
        v2f = v1f

      print 'velocidade da partícula 1 apos a colisão = ', v1f
      print 'velocidade da partícula 2 apos a colisão = ', v2f

      #--------------------- Desenho dos vetores --------------------------

      centro      = vis.vector(0.0,0.0,0.0)
      centroCena  = centro
      larguraCena = vis.vector(2.0, 2.0, 2.0)

      cena1 = vis.display(title="Antes da colisão", x=0, y=0,\
        width=400, height=400)

      vecV1i = vis.arrow(pos=centro, axis=(v1i,0,0), color=(1,0,0), \
                        shaftwidth=0.025, display=cena1)
      vecV2i = vis.arrow(pos=centro, axis=(v2i,0.0,0), color=(0,1,0), \
                        shaftwidth=0.035, display=cena1, opacity=0.1)

      cena2 = vis.display(title="Depois da colisão", x=410, y=0,\
        width=400, height=400, background=(1,1,1))
        
      vecV1f = vis.arrow(pos=centro, axis=(v1f,0,0), color=(1,0,0), \
                        shaftwidth=0.025, display=cena2)
      vecV2f = vis.arrow(pos=centro, axis=(v2f,0.0,0), color=(0,1,0), \
                        shaftwidth=0.035, display=cena2, opacity=0.1)
      

Exercício 03 (4 pts)

A função \[ f(x) = \frac{1.0}{1.0 - x} \label{eq1} \] pode ser aproximada pela serie \[ f(x) \approx x^0 + x^1 + x^2 + x^3 + \ldots + x^n = \sum_{i=0}^n x^i \label{eq2} \] desde que $|x| < 1.0$. Gostaríamos de avaliar quantos termos da serie são necessários para tornar a diferencia entre a equação $\ref{eq1}$ e a equação $\ref{eq2}$ menor que um dada tolerância, $\epsilon$.

Por exemplo, para um $\epsilon < 1.0\times 10^{-1}$ e $x = 0$ é necessário apenas um termo da série, enquanto que para $x=1/2$: \[ f(1.0/2.0) = \frac{1.0}{1.0 - 1.0/2.0} = 2.0 \nonumber \] e \[ f(1.0/2.0) \approx 1.0 + \left(\frac{1.0}{2.0}\right) + \left(\frac{1.0}{2.0}\right)^2 + \left(\frac{1.0}{2.0}\right)^3 + \left(\frac{1.0}{2.0}\right)^4 = 1.9375 \nonumber \] de onde vemos que a diferencia entre esses valores é \[ \left|2.0 - 1.9375\right| = 0.0625 < \epsilon \nonumber \]

Obviamente, o programa deve perguntar o valor de $x$ a ser avaliado e responder o numero mínimo de interações necessárias que verifique a condição previamente explicada.s

Solução

definimos \[ df = \left| f - f_e \right| \nonumber \] como sendo a diferencia entre o valor exato, $f_e$, da função no ponto $x$ e o valor obtido pela serie, $f$. A ideia é calcular num "loop" infinito cada termo da serie é ir adicionando ao valor de $f$ que inicialmente deve ser colocado como sendo igual a zero, a cada calculo de um novo elemento a equação anterior deve ser avaliada e quando a condição de saída do "loop" for confirmada $df \leq \epsilon$ então o "loop" é abandonado.

 
      #!/usr/bin/env python
      #-*- coding: utf-8 -*-

      import numpy as np
      import matplotlib.pyplot as plt


      x  = input('digite número menor do que 1: ')
      e  = input('digite a tolerância: ')
      fe = 1.0 / (1.0 - x)

      df = 10000.0
      i  = 1
      f  = 0.0
      while(df > e):
        f = f + x**i
        i = i + 1
        df = np.abs(f - fe)
      

Exercício 04 (6 pts)

Considere o seguinte sistema ideal, um camnhião de massa $M_c$ que carrega um quantidade de areia $m_a$. Devido a que a caçamba não foi fechada o camnhião perde a areia que transporta numa taxa linear, assim depois de passados $T$ minutos não há mais areia na caçamba. Admita que o motor do camnhião é um motor ideal que equivale a supor que o camnhião está sendo arrastrado com um força $F$ constante durante todo o percurso.
Faça um programa que pergunte $M_c$, $m_a$, $T$, $F$ e com esses dados construa um gráfico da velocidade em função do tempo com 20 pontos, desde $t=0$ até $t=1.5T$. O programa também deve guardar num arquivo chamado camnhiao_pos.dat os 20 dados de $t$ vs $x_c$ e camnhiao_ace.dat os 20 dados de $t$ vs $a_c$, onde $x_c$ e $a_c$ são a posição e aceleração do camnhião.

Solução

A ideia por trás deste problema é reduzir a massa total $M_t = M_c+m_a$ de forma linear segundo a equação \[ m = M_c+m_a - \frac{m_a}{T}t \nonumber \] onde $T$ é o tempo necessário para esvaziar a caçamba do caminhão e t é o tempo transcorrido desde o inicio do percurso do caminhão.

      #!/usr/bin/env python
      #-*- coding: utf-8 -*-

      import numpy as np
      import matplotlib.pyplot as plt

      mc = input('massa do caminhão: ')
      ma = input('massa de areia: ')
      T  = input('em quanto tempo o caminhão se esvazia (em minutos):  ')
      F  = input('qual é a força que atua sobre o caminhão: ')

      T  = T * 60.0
      tf = 1.5 * T
      vo = 0.0
      xo = 0.0

      t  = np.linspace(0, tf, 20)
      print t

      a = np.zeros_like(t)
      v = np.zeros_like(t)
      x = np.zeros_like(t)
      m = np.zeros_like(t)

      pos   = open('camnhiao_pos.dat', 'w')
      ace   = open('camnhiao_ace.dat', 'w')
      massa = open('camnhiao_mas.dat', 'w')

      for i in range(len(t)):
        if (t[i] > T):
          m[i] = mc
        else:
          m[i] = mc + ma - ma/T * t[i]

        a[i] = F / m[i]
        v[i] = vo + a[i] * t[i]
        x[i] = xo + vo * t[i] + 0.5 * a[i] * t[i]**2

        print >>pos, t[i], x[i]
        print >>ace, t[i], a[i]
        print >>massa, t[i], m[i]

      plt.plot(t, v)
      plt.show()
      

Exercício 05 (10 pts)

Atualmente sabemos que com toda certeza que o sistema solar é um sistema heliocêntrico similar ao proposto por Nicolás Copérnico graças às observação posteriores de Tycho Brahe, Johannes Kepler e por fim a teoria da gravitação universal de Newton. Mas essa concepção nem sempre foi universalmente aceito. Quase 2000 anos antes, Tolomeo de Alexandria propôs um modelo geocêntrico o qual explicava o movimento aparente dos planetas no céu chamados de epiciclos. No modelo tolemaico supõe-se que além da Terra estar no centro do sistema solar (e do universo) os planetas giravam entorno de um ponto central e era este ponto central que girava em volta da Terra. Mesmo está teoria sendo totalmente ineficaz em predizer qualquer coisa relativa ao movimento dos planetas se manteve válida por 2000 anos dada a sua enorme cargamento filosófico que dava sustentação (por isso que a física é uma ciência experimental). O nossos modelos atuais permitem entender sem muito esforço o porque desse movimentos aparente dos corpos celestes, a razão é que a Terra não é um sistema inercial e qualquer observação que realizemos a partir deste planeta deve levar em consideração o movimento deste em relação ao Sol.

O objetivo deste problema é reproduzir o movimento aparente de Marte visto desde a Terra. Para isso você vai precisar saber
a massa do Sol $M_S=1,99\times 10^{30}Kg$,
massa da Terra $M_T = 5,98\times 10^{24}kg$,
massa de Marte $M_M = 0,64\times 10^{24}kg$ e
a constante de gravitação universal $G=6,673\times 10^{-11}\frac{Nm^2}{kg^2}$

Como vocês sabem o movimento dos planetas é elíptico e a equação que descreve a posição do planeta em relação ao centro é

\[ r = \frac{a\left(1-e^2\right)}{1 + e\,\cos\theta} \nonumber \] onde $e$ é a excentricidade da orbita do planeta que para a Terra a excentricidade é \[ e_T = 0.0167, \nonumber \] e par Marte é \[ e_M = 0.0935 \nonumber \] O semi eixo maior da órbita da Terra \[ a_T = 149.60\times 10^9\,m, \nonumber \] e a órbita de Marte é \[ a_M = 227.92\times 10^9\,m, \nonumber \]

Para o calculo considere que a frequência angular, $\omega$, é constante (coisa que não é verdade) para ambos dos planetas, utilizando para isso que o período de rotação da Terra \[ T_T = 365.256\, dias, \nonumber \] e de Marte é \[ T_M = 686.980\; dias \nonumber \]

Bonus (5 pts)

Crie uma animação de como se veria o movimento de Marte desde a Terra.

Solução

Devido à falta de dados o problema foi eliminado, mas de qualquer forma coloco a disposição a ideia original. O objetivo é perceber que o movimento dos epiciclos observados de Marte tem origem no movimento relativo deste planeta em relação à Terra, por tanto podemos calcular ele via a seguinte equação \[ \vec{R}_{M/T} = \vec{R}_{M} - \vec{R}_T \nonumber \]

 
      #!/usr/bin/env python
      #-*- coding: utf-8 -*-

      import numpy as np
      import matplotlib.pyplot as plt
      import visual as vis

      G  = 6.673e-11
      Ms = 1.99e30
      Mt = 5.98e24
      Mm = 6.4e23

      et = 0.0167
      em = 0.0935

      Tm = 686.98 * (24.0*3600.0)
      Tt = 365.256* (24.0*3600.0)

      wm = 2.0 * np.pi / Tm
      wt = 2.0 * np.pi / Tt

      at = 149.60e6
      am = 227.92e6

      t     = np.linspace(0, 50.0*Tm, 5000)
      theta_t = wt * t
      theta_m = wm * t
      rm = am * (1.0 - em) / (1.0 + em * np.cos(theta_m))
      rt = at * (1.0 - et) / (1.0 + et * np.cos(theta_t))

      xm = rm * np.cos(theta_m)
      xt = rt * np.cos(theta_t)

      ym = rm * np.sin(theta_m)
      yt = rt * np.sin(theta_t)

      xmt = xm - xt
      ymt = ym - yt

      xmax = np.max(xmt)
      xmin = np.max(xmt)
      ymax = np.max(ymt)
      ymin = np.max(ymt)

      #plt.plot(xm, ym)
      #plt.plot(xt, yt)
      #plt.show()

      #--------------  Animação ----------------

      larguraCena = vis.vector(1.5*xmax, 1.5*ymax, 1.0)
      centro      = vis.vector(0,0,0)

      cena1 = vis.display(title="movimento dos planetas", x=0, y=0,\
        width=600, height=600, range=larguraCena)

      cena2 = vis.display(title="movimento de Marte desde Terra", x=610, y=0,\
        width=600, height=600, range=larguraCena)

      radio = xmax / 20.0
      Marte = vis.sphere(pos=(0, 0, 0), radius=radio, display=cena1, color=(1,0,0))
      Terra = vis.sphere(pos=(0, 0, 0), radius=1.25*radio, display=cena1, color=(0,0,1))
      MarteDesdeTerra = vis.sphere(pos=(0, 0, 0), radius=radio, display=cena2, color=(1,0,0))

      orbita = vis.curve( color = vis.color.cyan, radius=radio/4.0, display=cena2)


      for i in range(len(xmt)):
        vis.rate(30)
        novoCentro = vis.vector(xmt[i], ymt[i], 0)
        orbita.append( pos=novoCentro,  retain=150 )
        MarteDesdeTerra.pos = novoCentro
        
        Marte.pos = vis.vector(xm[i], ym[i], 0)
        Terra.pos = vis.vector(xt[i], yt[i], 0)Centro
      

O resultado esperado para o movimento relativo de Marte visto desde a Terra é