Apresentação

Destacado

 Notícias para os alunos:

Notícia do Colegiado do Curso de Física :

Sem notícias essa semana.

Horário de Atendimento aos alunos (inclusive Colegiado do Curso de Física) :

Terça – Feiras: 15h – 18h;

Quarta – Feira: 13h-15h (Apenas assuntos de aula ou entrega de documentos ou agendamentos feitos pelo coordenador);

Quintas – Feiras: 15h – 18h;

Acompanhe o blog semanalmente, pois em virtude de reuniões os dias e horários de atendimento pode sofrer alterações em algumas semanas.

Contato por email:

ramon.campos(arroba)ufes.br

troque o (arroba) pelo @, essa é uma medida para evitar os robozinhos que enviam mensagens aleatórias para meu email de trabalho.

Disciplinas do semestre 2015/2:

Veja abaixo a data da última atualização.

Mecânica Clássica – 20 de Outubro

Seminários de Física Contemporânea  – Em Construção

E + abaixo temos notícias…

 

Motor Boxer e a Física

Olá a todos,

hoje vamos mostrar o Motor Boxer e como podemos entender vários conceitos físicos estudando esse motor. Centro de Gravidade (igual ao Centro de Massa nesse caso),  Vibrações entre outros conceitos;

O vídeo apesar de ser em inglês, é muito ilustrativo e auto explicativo;

Disciplinas do semestre 2015/1

Notícias para os alunos:

Em breve notícia sobre o motor boxer e seu desempenho campeão.

Horário de Atendimento aos alunos (inclusive Colegiado do Curso de Física) :

Semana do 6 a 10 de Julho:

Quarta – Feiras: 9:30h – 11:00h;

Quarta – Feira: 14:30h-16:30h;

Quintas – Feiras: 12:30h – 15:30h;

Sexta – Feira: 9h – 11h;

Acompanhe o blog semanalmente, pois em virtude de reuniões os dias e horários de atendimento pode sofrer alterações em algumas semanas.

Nada te impede de me procurar em outros dias e horários, mas nestes eu estou a disposição para tende-los.

Contato por email:

ramon.campos(arroba)ufes.br

troque o (arroba) pelo @, essa é uma medida para evitar os robozinhos que enviam mensagens aleatórias para meu email de trabalho.

Disciplinas do semestre 2015/1:

Veja abaixo a data da última atualização.

Eletromagnetismo 1 – 03 de Julho

Física Básica –  03 de Julho

E + abaixo temos notícias…

 

Transmissão de motocicletas

Esse aqui é para entender as diferenças de funcionamento, vantagens e desvantagem dos diferentes tipos de transmissão das motocicletas.

Vídeo do Youtube sobre o assunto.

Algumas motos para pensar no assunto dos custos de aquisição, e ver que não é só a transmissão que manda aqui, motor, acessórios e sistemas como ABS fazem muita diferença.

Yamaha – XT 660 – Motor 48 cv por R$ 28.870,00, ou um BMW G 650 motor de 50 cv por R$29.000,00. Mas poderia-se pagar R$33.680 em uma Ténéré 660 (Com ABS) ou R$ 30.990,00  na XL 700V Transalp (sem ABS) mas com motor 60 cv – Todas elas com corrente.

Por usa vez Hoda Shadow 750 – Motor 45,5 cv – sai a R$29.900,00 (Sem abs), mas com eixo cardan;

Finalmente se pode comprar uma Harley Davison Iron 883 – motor de ~51cv por exatos R$34.900,00 e ela vem com correia dentada, mas é uma Harley e no geral sempre vem com esse sistema .

Python – Aula 4

Computação Simbólica no IPython Notebook usando sympy

Disciplina “Computação no Ensino de Física” (CEF) da UFES/Alegre

Curso de Licenciatura em Física – Turma 2014/2

Autoria : Prof. Ramón Giostri Campos – 2014/2.

Veremos abaixo como usar o IPython notebook e o pacote matplotlib para fazer gráficos profissionais

Carregar o IPyhon notebook e os pacotes iniciais:

Digite $ipython notebook

no terminal do seu sistema operacional (shell no linux e prompt no windowns);

Abaixo vamos chamar os principais pacotes, separando eles por contexto

In [1]: %pylab inline

Populating the interactive namespace from numpy and matplotlib

%pylab inline = chama o pylab e coloca os graficos embutidos;

In [2]:  import numpy as np
In [3]:  import scipy as sp
In [4]:  import matplotlib as mpl
In [5]:  import matplotlib.pyplot as plt
In [6]:  #ion() #isso liga o modo interativo do matplotlib, usamos só se o inline nao for usado no primeiro comando!

Exercício 1

Escreva um comentário ilustrativo sobre cada uma das linhas de input, assim como foi feito abaixo da primeira (e na ultima) linha da sequência de comandos;

# %pylab inline = chama o pylab e coloca os graficos embutidos;
(...)
#ion() #isso liga o modo interativo do matplotlib, usamos so se o inline nao for usado no primeiro comando!

Pintando curvas, pontos, e figuras geométricas;

Um gráfico bruto de uma reta, dado as coordenadas;

In [7]:   plt.plot([1.0,3.4],[4.7,5.9])
Out[7]:  [<matplotlib.lines.Line2D at 0x3dd5ad0>]

A ideia do matplotlib é usar uma interpolação linear dos pontos para gerar os gráficos, isso para retas não tem problema, mas para curvas…

In [8]:
x1=np.arange(-2,3) #domínio
y1=(x1)**2 #Imagem
plt.plot(x1,y1)
Out[8]:  [<matplotlib.lines.Line2D at 0x41241d0>]

Claramente temos de almentar o número de pontos no intervalo para ter uma melhor interpolação…+ é bom não exagerar! Além disso podemos fazer uma coisa sistematizada mais elegante!

In [9]:
xi=-2. # limite inferior do domínio
xf=2. #liminte superior do domínio
numd=100. #número de divisões
dx=(xf-xi)/numd #intervalo das divisões
x1=np.arange(xi,xf+dx,dx) # Gera o vetor domínio
y1=(x1)**2 # Gera o vetor imagem
plt.plot(x1,y1)
Out[9]:  [<matplotlib.lines.Line2D at 0x433d550>]

Poderíamos encurtar a coisa usando o comando…

linspace

 

Vamos aprender mais um pouco sobre esse comando…

In [13]:  help(linspace)
numpy.core.function_base:

linspace(start, stop, num=50, endpoint=True, retstep=False)
    Return evenly spaced numbers over a specified interval.

    Returns `num` evenly spaced samples, calculated over the
    interval [`start`, `stop` ].

    The endpoint of the interval can optionally be excluded.

    Parameters
    ----------
    start : scalar
        The starting value of the sequence.
    stop : scalar
        The end value of the sequence, unless `endpoint` is set to False.
        In that case, the sequence consists of all but the last of ``num + 1``
        evenly spaced samples, so that `stop` is excluded.  Note that the step
        size changes when `endpoint` is False.
    num : int, optional
        Number of samples to generate. Default is 50.
    endpoint : bool, optional
        If True, `stop` is the last sample. Otherwise, it is not included.
        Default is True.
    retstep : bool, optional
        If True, return (`samples`, `step`), where `step` is the spacing
        between samples.

    Returns
    -------
    samples : ndarray
        There are `num` equally spaced samples in the closed interval
        ``[start, stop]`` or the half-open interval ``[start, stop)``
        (depending on whether `endpoint` is True or False).
    step : float (only if `retstep` is True)
        Size of spacing between samples.

    See Also
    --------
    arange : Similar to `linspace`, but uses a step size (instead of the
             number of samples).
    logspace : Samples uniformly distributed in log space.

    Examples
    --------
    >>> np.linspace(2.0, 3.0, num=5)
        array([ 2.  ,  2.25,  2.5 ,  2.75,  3.  ])
    >>> np.linspace(2.0, 3.0, num=5, endpoint=False)
        array([ 2. ,  2.2,  2.4,  2.6,  2.8])
    >>> np.linspace(2.0, 3.0, num=5, retstep=True)
        (array([ 2.  ,  2.25,  2.5 ,  2.75,  3.  ]), 0.25)

    Graphical illustration:

    >>> import matplotlib.pyplot as plt
    >>> N = 8
    >>> y = np.zeros(N)
    >>> x1 = np.linspace(0, 10, N, endpoint=True)
    >>> x2 = np.linspace(0, 10, N, endpoint=False)
    >>> plt.plot(x1, y, 'o')
    [<matplotlib.lines.Line2D object at 0x...>]
    >>> plt.plot(x2, y + 0.5, 'o')
    [<matplotlib.lines.Line2D object at 0x...>]
    >>> plt.ylim([-0.5, 1])
    (-0.5, 1)
    >>> plt.show()

De qualquer forma isso é uma questão de gosto!

Vamos voltar aos nosso gráficos e tentar ver os pontos de um e do outro para enteder melhor o que se passa!

In [12]:
#Código para gerar os gráficos sobrepostos
#
xi=-2. # limite inferior do domínio
xf=2. #liminte superior do domínio
#
# 5 Pontos, a subdivisão e os vetores de pontos correspondente
#
numd=5. #número de divisões
dx=(xf-xi)/numd #intervalo das divisões
x1=np.arange(xi,xf+dx,dx)
y1=(x1)**2
#
# 20 Pontos a subdivisão e vetor de pontos correspondente
#
numd=10. #número de divisões
dx=(xf-xi)/numd #intervalo das divisões
x2=np.arange(xi,xf+dx,dx)
y2=(x2)**2
plt.plot(x1,y1,'r-',x2,y2,'go',x2,y2,'b-')
#Isso sim é código mastigado!!!
Out[12]:
[<matplotlib.lines.Line2D at 0x4a05ed0>,
 <matplotlib.lines.Line2D at 0x4a093d0>,
 <matplotlib.lines.Line2D at 0x4a09890>]

Exercício 2:

Tente descobri o que significam as opções ‘r-‘, ‘go’ e ‘b-‘ no gráfico acima.

Vamos as figuras geométricas agora:

Triângulos:

In [13]:  plt.plot([-1,0.5,1,-1],[0,2,1,0])
Out[13]: <matplotlib.lines.Line2D at 0x4c33390>]

Você provavelmente notou que a sintax do python é diferente do que estamos acostumados lá da Geometria Analítica:

Em GA temos pares (x,y) e com 3 desses pares (não colineares) fazemos o triângulo.

No Matplotlib é fornecido duas listas, uma com todos os x’s e outra com todos os y’s;

Mas com um pouco de programação podemos tornar o matplotlib mais intuitivo.

In [21]:
#Script para fazer triângulos usandos os pares ordenados.
#
#Escreva os pares ordenados abaixo:
pares=([-1,0],[0.5,2],[1,1])
print pares
pares_f=pares+(pares[0],) # essa operação inclui um termo adicional igual ao primeiro termo na ultima posição;
print pares_f
pares_t=np.transpose(pares_f) # essa operação é a transposição da matriz (uma tupla do python que se comporta como matriz;
print pares_t
plt.plot(*pares_t)
([-1, 0], [0.5, 2], [1, 1])
([-1, 0], [0.5, 2], [1, 1], [-1, 0])
[[-1.   0.5  1.  -1. ]
 [ 0.   2.   1.   0. ]]
Out[21]: [<matplotlib.lines.Line2D at 0x599a390>]
In [22]:
#Script para fazer triângulos usandos os pares ordenados.
#
#Escreva os pares ordenados abaixo:
pares=([-1,0],[0.5,2],[1,1],[2,1])
print pares
pares_f=pares+(pares[0],) # essa operação inclui um termo adicional igual ao primeiro termo na ultima posição;
print pares_f
pares_t=np.transpose(pares_f) # essa operação é a transposição da matriz (uma tupla do python que se comporta como matriz;
print pares_t
plt.plot(*pares_t)
([-1, 0], [0.5, 2], [1, 1], [2, 1])
([-1, 0], [0.5, 2], [1, 1], [2, 1], [-1, 0])
[[-1.   0.5  1.   2.  -1. ]
 [ 0.   2.   1.   1.   0. ]]
Out[22]:  [<matplotlib.lines.Line2D at 0x5b76a50>]

Aqui também podemos tentar algo mais elaborado, com cores e legenda.

In [23]:
from __future__ import unicode_literals # isso permite incluir acentos nos gráficos. E era um exercício em outra aula.
In [24]:
#Script para fazer triângulos usandos os pares ordenados e com legendas.
#
#Escreva os pares ordenados abaixo:
pares=([-1,0],[0.5,2],[1,1])
#
#
pares_f=pares+(pares[0],) # essa operacao inclui um termo adicional igual ao primeiro termo na ultima posicao;
pares_t=np.transpose(pares_f) # essa operacao e a transposicao da matriz (uma tupla do python que se comporta como matriz;
plt.plot(pares_t[0],pares_t[1],'ro', label='Vértices')
plt.plot(pares_t[0],pares_t[1],label='Lados')
plt.title('Triângulo')
plt.ylim((-0.1,2.1)) #expande a area visivel do eixo y
plt.xlim((-1.1,1.1)) #expande a area visivel do eixo y
plt.legend(loc='upper left') # legenda inclusive especificando a posição na figura
plt.show()

Pintando pontos, e outros símbolos;

Se você passou por isso tudo, acredito que já tenha entendido como plotar pontos.

Vamos resumir para facilitar:

1 – Tome as duas listas de mesmo tamanho, uma representando os valores x’dos pontos e outra para os y’dos pontos (evidente que só pode ter número na lista).

2 – Insira as duas listas em sequência no comando PLOT (no nosso CONTEXTO plt.plot).

3 – Acrescente ‘o’ para o estilo de marcação de pontos, “-” para linha e “–” para linha tracejada.

Vamos mais  exemplos;

Pontos aleatórios:

In [29]:
num_p=100 #total de pontos
x=np.random.random_sample(num_p) #faz uma amostragem aleatório dos pontos para o x
y= np.random.random_sample(num_p) #idem para o y
plt.plot(x,y,'o')
Out[29]: [<matplotlib.lines.Line2D at 0x6583e50>]

Exercício 3:

Se em lugar de bolinhas azuis, tivéssemos sinais ‘+’ vermelhos, como faríamos isso?

Pontos de figuras regulares:

In [35]:
num=15 #numero de lados da figura
dt= 2*np.pi/num
#usamos a matriz de rotação sobre o vetor (1,1) para gerar os pares x e y
lista_x=map(lambda t:(np.cos(t)+np.sin(t)),arange(0,2*pi+dt,dt))
lista_xf=append(lista_x,lista_x[0])
lista_y=map(lambda t:(-np.sin(t)+np.cos(t)),arange(0,2*pi+dt,dt))
lista_yf=append(lista_y,lista_x[0])
#expande a area visivel do eixo y
plt.ylim((-1.5,1.5)) #usamos o valor de 1.5, pois o maior valor que pode ser alcançado aqui é raiz de 2;
plt.xlim((-1.5,1.5))
#
plt.plot(lista_x,lista_y) #linhas
plt.plot(lista_x,lista_y,'o',label='Vértices') # vertices
plt.legend(loc='upper left')
# acerta a proporção dos eixos, você lembra que isso foi um problema na aula passada usando SymPy!
plt.axes().set_aspect('equal')
In [21]:
num=5#numero de lados da figura
dt= 2*np.pi/num
#usamos a matriz de rotação sobre o vetor (1,1) para gerar os pares x e y
lista_x=map(lambda t:(np.cos(t)+np.sin(t)),arange(0,2*pi+dt,dt))
lista_xf=append(lista_x,lista_x[0])
lista_y=map(lambda t:(-np.sin(t)+np.cos(t)),arange(0,2*pi+dt,dt))
lista_yf=append(lista_y,lista_x[0])
#expande a area visivel do eixo y
plt.ylim((-1.5,1.5)) #usamos o valor de 1.5, pois o maior valor que pode ser alcançado aqui é raiz de 2;
plt.xlim((-1.5,1.5))
#
plt.plot(lista_x,lista_y)
# acerta a proporção dos eixos, você lembra que isso foi um problema na aula passada usando SymPy!
plt.axes().set_aspect('equal')

Movimento Browniano – Caminho aletório – Andar do Bêbado. Usando Scatter Plot;

Comecemos com um caso unidimensional que não precisa de gráfico.

In [41]:
num_p= 1000 #número de passos
pa=0
lista_p=[0]
for i in range(num_p-1):
    pn = 1.0*random_integers(0,1) # isso indica que ele pode dar um passo a frente(1) ou ficar parado (0)
    #alternativamente pode-se colocar -1 no lugar do ZERO, assim dizemos que o bêbado pode voltar um passo;
    pa=pn+pa
    lista_p=append(lista_p,pa)
print 'A fração do caminho percorrido ao final de', +len(lista_p), 'passos foi',+ lista_p[-1]/len(lista_p)
A fração do caminho percorrido ao final de 1000 passos foi 0.499
Agora vamos construir duas listas usando esse procedimento e fazer o problema bidimensional;
In [47]:
num_p= 200 #número de passos
pax=0
pay=0
lista_px=[0]
lista_py=[0]
for i in range(num_p-1):
    pnx = 1.0*random_integers(0,1) # isso indica que ele pode dar um passo a frente(1) ou ficar parado (0)
    #alternativamente pode-se colocar -1 no lugar do ZERO, assim dizemos que o bêbado pode voltar um passo;
    pax=pnx+pax
    lista_px=append(lista_px,pax)
    pny = 1.0*random_integers(-1,1)
    pay=pny+pay
    lista_py=append(lista_py,pay)
plt.scatter(lista_px,lista_py)
Out[47]:
<matplotlib.collections.PathCollection at 0x7201350>

Note que o código é muito parecido com o anterior, apenas duplicamos as funções do loop e no lugar de mostrar o print, fazemos o gráfico, no caso usando um scatter plot (genuino gráfico de pontinhos);

In [52]:
num_p= 10000 #número de passos
pax=0
pay=0
lista_px=[0]
lista_py=[0]
lista=[0]
for i in range(num_p-1):
    lista=append(lista,i)
    pnx = 1.0*random_integers(-1,1) # isso indica que ele pode dar um passo a frente(1) ou ficar parado (0)
    #alternativamente pode-se colocar -1 no lugar do ZERO, assim dizemos que o bêbado pode voltar um passo;
    pax=pnx+pax
    lista_px=append(lista_px,pax)
    pny = 1.0*random_integers(-1,1)
    pay=pny+pay
    lista_py=append(lista_py,pay)
plt.scatter(lista_px,lista_py,s=50,c=lista) #pinta todos os pontos seguindo a distribuição do "i" para as cores;
plt.scatter(lista_px[0],lista_py[0],s=120, c= 'g') #pinta o ultimo ponto maior e em vermelho
plt.scatter(lista_px[-1],lista_py[-1],s=120, c= 'r') #pinta o ultimo ponto maior e em vermelho
plt.show()

Gráficos de funções, de uma ou mais variáveis

Como vimos acima, tudo gira em torno de criar as listas corretas, uma de domínio (x) e uma para a imagem (f(x));

Usar em sequência para plotar os gráficos, vamos ver alguns exemplos:

Funções trigonométricas

In [53]:
xi=-2*pi # limite inferior do domínio
xf=2*pi #liminte superior do domínio
#
numd=500. #número de divisões
dx=(xf-xi)/numd
x_lista=np.arange(xi,xf+dx,dx)

Seno e cosseno

In [54]:
y_sen=np.sin(x_lista)
y_cos=np.cos(x_lista)
plt.plot(x_lista,y_cos,'r--', label='Cos(x)')
plt.plot(x_lista,y_sen,'b-',label='Sen(x)')
plt.title('Trigonometria')
plt.xlabel('x')
plt.ylabel('f (x)')
plt.legend(loc='upper right') # legenda inclusive especificando a posição na figura
plt.show()
In [55]:
y_sen=np.sin(10.0*x_lista)*np.cos(x_lista)
y_cos=np.cos(x_lista)
plt.plot(x_lista,y_cos,'r--', label='+ - Sen(x)')
plt.plot(x_lista,y_sen,'b-',label='Sen(10x)*Cos(x)')
plt.plot(x_lista,-y_sen,'b-') # sem o label = sem aparecer na legenda!
plt.title('Trigonometria')
plt.xlabel('x')
plt.ylabel('f (x)')
plt.legend(loc='upper right') # legenda inclusive especificando a posição na figura
plt.show()

Gráficos paramétricos

Curiosamente, dada a forma como os gráficos são feitos em Python isso se torna trivial…

In [56]:
ti=0 # limite inferior do parâmetro
tf=10*pi #liminte superior do parâmetro
#
numd=500. #número de divisões
dt=(tf-ti)/numd
t_lista=np.arange(ti,tf+dt,dt) #Lista dos pontos do parâmetro
x_cost=t_lista*np.cos(t_lista) #como x se relaciona com o parâmetro
y_sent=t_lista*np.sin(t_lista) #como y se relaciona com o parâmetro
plt.plot(x_cost,y_sent)
plt.axes().set_aspect('equal')
Out[56]:
[<matplotlib.lines.Line2D at 0x7228f90>]
In [64]:
# Usando teoria de coordenadas polares;
ti=0 # limite inferior do parâmetro
tf=10*pi #liminte superior do parâmetro
#
numd=1000. #número de divisões
dt=(tf-ti)/numd
t_lista=np.arange(ti,tf+dt,dt) #Lista dos pontos do parâmetro
r_polar = cos(3*t_lista) #Raio das coordenadas pelares
x_cost=r_polar*np.cos(t_lista) # x em termos de r e t
y_sent=r_polar*np.sin(t_lista) # y em termosde r e t
plt.plot(x_cost,y_sent)
plt.axes().set_aspect('equal')
Out[64]:
[<matplotlib.lines.Line2D at 0x88e4b50>]
In [68]:
# Usando teoria de coordenadas polares;
ti=0 # limite inferior do parâmetro
tf=10*pi #liminte superior do parâmetro
#
numd=1000. #número de divisões
dt=(tf-ti)/numd
t_lista=np.arange(ti,tf+dt,dt) #Lista dos pontos do parâmetro
r_polar = cos(t_lista)*sin(t_lista) #Raio das coordenadas pelares
x_cost=r_polar*np.cos(t_lista) # x em termos de r e t
y_sent=r_polar*np.sin(t_lista) # y em termosde r e t
plt.plot(x_cost,y_sent)
plt.axes().set_aspect('equal')
Out[68]:
[<matplotlib.lines.Line2D at 0x923e110>]

Outros gráficos 2D interessantes:

Regiões:

In [69]:
x = np.linspace(0, 4,100)
y = np.sin(4*np.pi*x) * np.exp(-0.5 * x)
plt.fill(x,y, 'b')
plt.show()

Podemos inclusive misturar as coisas…

In [70]:
x = np.linspace(0, 6,100)
y = np.sin(4*np.pi*x) * np.exp(-0.5 * x)
plt.fill(x,y, 'r', label='Oscilação')
y_e= np.exp(-0.5 * x)
plt.plot(x,y_e, 'b', label='+ - Envoltória')
plt.plot(x,-y_e,'b')
plt.legend(loc='upper right')
plt.show()

Histograma:

In [71]:
mu = 100 # media da disbribuicao
sigma = 15 # desvio padrao
x = mu + sigma * np.random.randn(10000) # amostra de 10 mil numeros
num_bins = 50
# esse comando gera 3 resutados;
n, bins, patches = plt.hist(x, num_bins, normed=1, facecolor='green', alpha=0.4)
# fazendo o ajuste da curva normal ao histograma
y = mlab.normpdf(bins, mu, sigma)
plt.plot(bins, y, 'r--')
#Titulos
plt.xlabel('Valores Prováveis')
plt.ylabel('Probabilidade')
plt.title(r'Histograma de IQ: $\mu=100$, $\sigma=15$')
# chega os graficos para o lado
plt.subplots_adjust(left=0.15)
plt.show()

 

Visualização de figuras de mais variáveis;

 

In [72]:
from mpl_toolkits.mplot3d import Axes3D
#para fazer gráficos 3D
In [74]:
mpl.rcParams['legend.fontsize'] = 10
fig = plt.figure()
ax = fig.gca(projection='3d')
theta = np.linspace(-4 * np.pi, 4 * np.pi, 100)
z = np.linspace(-2, 2, 100)
x =  np.sin(theta)
y =  np.cos(theta)
ax.plot(x, y, z, label='Curva Espiral')
ax.legend()
plt.show()
In [77]:
mpl.rcParams['legend.fontsize'] = 10
fig = plt.figure()
ax = fig.gca(projection='3d')
theta = np.linspace(0, 4 * np.pi, 200)
r = theta
x =  r*np.sin(4*theta)
y =  r*np.cos(4*theta)
z = theta**0.5
ax.plot(x, y, z, label='Tornado')
ax.legend()
plt.show()
In [79]:
mpl.rcParams['legend.fontsize'] = 10
fig = plt.figure()
ax = fig.gca(projection='3d')
theta = np.linspace(0, 4 * np.pi, 200)
r = theta
x =  r*np.sin(4*theta)
y =  r*np.cos(4*theta)
z = theta
ax.plot(x, y, z, label='Tornado')
ax.plot(0.5*x, 0.5*y, 3*z, label='Outro Tornado')
ax.legend()
plt.show()
In [91]:
mpl.rcParams['legend.fontsize'] = 10
fig = plt.figure()
ax = fig.gca(projection='3d')
theta = np.linspace(0, 4 * np.pi, 200)
r = theta
x =  r*np.sin(4*theta)
y =  r*np.cos(4*theta)
z = theta
ax.plot(x, y, z, label='Tornado')
plt.plot(x, y, label='Sombra do Tornado')
ax.plot(x,0.*y+15,z, label='Outra Sombra do Tornado')
#ax.legend()
plt.show()

Campo Vetorial

In [99]:
Y, X = np.mgrid[-3:3:100j, -3:3:100j]
U = -1 - X**2 + Y
V = 1 + X - Y**2
speed = np.sqrt(U*U + V*V)
plt.streamplot(X, Y, U, V, color=U, linewidth=2, cmap=plt.cm.autumn)
plt.colorbar()
f, (ax1, ax2) = plt.subplots(ncols=2)
ax1.streamplot(X, Y, U, V, density=[0.5, 1]) 
lw = 5*speed/speed.max()
ax2.streamplot(X, Y, U, V, density=0.6, color='k', linewidth=lw)
plt.show()

 

In [39]:
#Campo elétrico da carga elétrica pontual!
Y, X = np.mgrid[-3:3:100j, -3:3:100j] #Pontos do vetor r
U = X/(X**2 +Y**2)**1.5 # Coordenada X do Campo
V =  Y/(X**2 +Y**2)**1.5 # Coordenada Y do Campo
plt.streamplot(X, Y,U,V, color=U, linewidth=2, cmap=plt.cm.autumn)
plt.show()
In [94]:
#Campo magnético do fio retilíneo
Y, X = np.mgrid[-3:3:100j, -3:3:100j] #Pontos do vetor r
U = -Y/(X**2 +Y**2)**1.5 # Coordenada X do Campo
V =  X/(X +Y**2)**1.5 # Coordenada Y do Campo
plt.streamplot(X, Y,U,V, color=U, linewidth=2, cmap=plt.cm.autumn)
plt.show()

Animação

As vezes é importante rodar as coisas em um script separado, abaixo temos um exemplo!

In [95]:
#Retire os comentarios das linhas 2 abaixo apenas nos scripts. Nao coloque acentos no comentario do script.
#import numpy as np
#import matplotlib.pyplot as plt

import matplotlib.animation as animation
fig = plt.figure()
def f(x, y):
    return np.sin(x) + np.cos(y)
x = np.linspace(0, 2 * np.pi, 120)
y = np.linspace(0, 2 * np.pi, 100).reshape(-1, 1)
im = plt.imshow(f(x, y), cmap=plt.get_cmap('jet'))
def updatefig(*args):
    global x,y
    x += np.pi / 15.
    y += np.pi / 20.
    im.set_array(f(x,y))
    return im,
ani = animation.FuncAnimation(fig, updatefig, interval=50, blit=True)
plt.show()
In [ ]:

Outra coisa que pode ficar melhor de ver com script separado…

In [111]:
# Retire os comentarios das linhas 3 abaixo apenas nos scripts.
#import numpy as np
#from mpl_toolkits.mplot3d import Axes3D
#import matplotlib.pyplot as plt
# Nao coloque acentos no comentario do script.
fig = plt.figure()
ax = fig.gca(projection='3d')
num_p= 1000 #numero de passos
pax=0
pay=0
paz=0
lista_px=[0]
lista_py=[0]
lista_pz=[0]
lista=[0]
for i in range(num_p-1):
    lista=np.append(lista,i)
    #Numeros do x
    pnx = 1.0*np.random.random_integers(-1,1) # isso indica que ele pode dar um passo a frente(1) ou ficar parado (0)
    #alternativamente pode-se colocar -1 no lugar do ZERO, assim dizemos que o bebado pode voltar um passo;
    pax=pnx+pax
    lista_px=np.append(lista_px,pax)
    #Numeros do y
    pny = 1.0*np.random.random_integers(-1,1)
    pay=pny+pay
    lista_py=np.append(lista_py,pay)
    #Numeros do z
    pnz = 1.0*np.random.random_integers(-1,1)
    paz=pnz+paz
    lista_pz=np.append(lista_pz,paz)
ax.scatter(lista_px,lista_py,lista_pz,s=15,c=lista,marker='+') #pinta todos os pontos seguindo a distribuicao do "i" para as cores;
ax.scatter(lista_px[0],lista_py[0],lista_pz[0],s=400, c= 'g') #pinta o ultimo ponto maior e em vermelho
ax.scatter(lista_px[-1],lista_py[-1],lista_pz[-1],s=400, c= 'r') #pinta o ultimo ponto maior e em vermelho
plt.show()

:

Python – Aula 3

Computação Simbólica no IPython Notebook usando sympy

Disciplina “Computação no Ensino de Física” (CEF) da UFES/Alegre

Curso de Licenciatura em Física – Turma 2013/2

Autoria : Prof. Ramón Giostri Campos – 2014/1.

Veremos abaixo como usar o IPython notebook para realizar cálculos simbólicos, usando o pacote sympy, no processo vamos exercitar a definição de funções (def) em python e gráficos usando o próprio sympy.

Carregar o IPyhon notebook e os pacotes iniciais:

Digite $ipython notebook

no terminal do seu sistema operacional (shell no linux e prompt no windowns);

Abaixo vamos chamar os principais pacotes

In [1]:
%pylab inline
Populating the interactive namespace from numpy and matplotlib

É importante usar os pacotes compatíveis, como estamos usando o ipython novo, devemos usar o sympy novo.

No ubuntu uma solução é fazer download do sympy mais novo no sítio:

http://sympy.org/pt/download.html

Dai descompactamos e colocamos o arquivo que queremos executar (file.ipynb) dentro do diretório do recém descompactado.

In [2]:
import sympy as sym # importa as funções do sympy em um contexto diferenciado (no caso sym)
from sympy import *

A mudança de contexto presente na primeira linha é importante para não sobrepor as funções vindas do sympy com as funções nativas do python e do pylab.

Executamos agora um comando para ver qual a versão do sympy, se tudo foi feito corretamente, o número que aparece aqui é o mesmo que está no pacote baixado do site e descompactado posteriormente;

In [3]:
from sympy import __version__
__version__
Out[3]:
'0.7.4.1'

Primeiro estágio: Operações básicas e considerações iniciais

In [4]:
3*pi/2 + exp(I*x) / (x**2 + y)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-4-057106a93d7b> in <module>()
----> 1 3*pi/2 + exp(I*x) / (x**2 + y)

NameError: name 'x' is not defined

É importante dizer para o sympy quais as objetos são SÍMBOLOS;

In [5]:
x, y, z, t = symbols("x y z t") # Define símbolos para variáveis genericas
k, m, n = symbols("k m n", integer=True)  # Define símbolos para variáveis inteiras
f, g, h = map(Function, 'fgh') # define símbolos para funções
In [6]:
3*pi/2 + sym.exp(I*x) / (x**2 + y)
Out[6]:
3*pi/2 + exp(I*x)/(x**2 + y)

Agora o Python entendeu o que fazer, pois x e y são símbolos.

Mas essa forma 1D é difícil de ler e muito feia…vamos melhorar isso.

Para isso devemos carregar o comando abaixo:

In [7]:
init_printing()

Vamos repetir o nosso teste;

In [8]:
3*pi/2 + sym.exp(I*x) / (x**2 + y)
#Muito mais bonito!
Out[8]:
3π2+eixx2+y

Algumas ações que podem ser feitas com expressões simbólicas: expandir (expand), simplificar(simplify) e resolver(solve);

 

In [9]:
eq = ((x+y)**2 * (x+1))
eq
Out[9]:
(x+1)(x+y)2
In [10]:
expand(eq)
Out[10]:
x3+2x2y+x2+xy2+2xy+y2
In [11]:
a = 1/x + (x*sin(x) - 1)/x
print(a) #assim fica feito
display(a) #assim fica bonito
In [12]:
simplify(a)
Out[12]:
sin(x)
In [13]:
eq2 = Eq(x**3 - 3*x**2 + 4*x + 8., 0)
display(eq2)
x3−3×2+4x+8.0=0

Notem que a equação gerada usa o comando “Eq”, para definir a igualdade;

In [14]:
solve(eq2,x)
Out[14]:
[−1.0,2.0−2.0i,2.0+2.0i]

Exercícios 1: Tente resolver equações de ordem maior que 3 e equações transcendentais com o solve para ver o que ocorre;

Terceiro esstágio: Física B;

Vamos começar com o sistema massa mola: F=-kx; A segunda lei de Newtom para ele nos gera ma+k*x=0, que em termo de derivadas e já dividindo pela massa fica:

In [15]:
w0=symbols("w_0",real=True)

Atentem para a Flag “real”, isso significa que o símbolo w0 é REAL;

In [16]:
eqn = Eq(Derivative(x(t),t,t) + w0**2*x(t), 0) # já trocamos k/m por w0**2
display(eqn)
dsolve(eqn, x(t))

E se houvesse atrito, como ficaria.

In [17]:

b=symbols("b",real=True)
In [18]:
eqn2 = Eq(Derivative(x(t),t,t) +b*Derivative(x(t),t)+ w0**2*x(t), 0) # já trocamos k/m por w0**2
display(eqn2)
dsolve(eqn2, x(t))
Out[18]:

Não se assunste, isso é feio assim pois o sympy não sabe quem é maior, b ou w0

Se b e w0 forem conhecidos teremoss algo relativamente simples;

In [19]:
eqn2_crit = Eq(Derivative(x(t),t,t) +4*Derivative(x(t),t)+ (2**2)*x(t), 0) # já trocamos w0=2 e b=4
display(eqn2_crit)
dsolve(eqn2_crit, x(t))
Out[19]:
In [20]:
eqn2_super = Eq(Derivative(x(t),t,t) +5*Derivative(x(t),t)+ (2**2)*x(t), 0) # já trocamos w0=2 e b=5
display(eqn2_super)
dsolve(eqn2_super, x(t))
Out[20]:
In [21]:
eqn2_sub = Eq(Derivative(x(t),t,t) +3*Derivative(x(t),t)+ (2**2)*x(t), 0) # já trocamos w0=2 e b=3
display(eqn2_sub)
dsolve(eqn2_sub, x(t))
Out[21]:

E se houvesse uma força externa, por exemplo senoidal!

In [22]:
wp, A=symbols("w_p A",real=True)
In [23]:
eqn3 = Eq(Derivative(x(t),t,t) + w0**2*x(t)-A*sin(wp*t), 0) # já trocamos k/m por w0**2
display(eqn3)
sol3=dsolve(eqn3, x(t))
display(sol3)

Exercício 2: Descubra se é possível reaproveitar o que é gerado do solve e do dsolve sem necessariamente digitar as equações manualmente;

 

Terceiro estágio: Cálculo B;

Primeiro vejamos uma derivada, já usamos ela acima nas equações diferencias.

Usamos os comandos “Derivative” para escrever a derivada e o comando “diff” para realizar a derivação;

In [24]:
Derivative(cos(m*x**k)**10,x)
Out[24]: Esse mostra a equação
In [25]:
diff(cos(m*x**k)**10,x)
Out[25]: Esse resolve a equação

Note que o comando “Derivative” é esperto, e se tivermos uma expressão com duas VARIÁVEIS, ele imprime na tela uma derivada parcial em lugar da derivada ordinária;

In [26]:
Derivative(cos(y*x**k),x)
Out[26]:
∂cos(xky)/∂x

Agora vamos a integral, não necessariamente uma simples.

Usamos o comando “Integral” para escrever a INTEGRAL sem resolve-la e o comando “integrate” para INTEGRAR e ver o resultado;

In [27]:
Integral(sin(k*x)*sin(m*x),(x,-pi/2,pi/2))
Out[27]: Esse mostra a equação
In [28]:
integrate(sin(k*x)*sin(m*x),(x,-pi/2,pi/2))
Out[28]: Esse resolve a equação

Essa é uma das integrais da base da Série de Fourier, usada nas séries de Fourier, muito importante para resolver problemas de MECÂNICA CLÁSSICA e ELETROMAGNETISMO.

Exercício 3: Pesquise quais as outras integrais dessa base e use o sympy para resolve-las;

Vamos a um exemplo mais trabalhoso, vamos investigar as Séries de Taylor:

Séries de Taylor (e outros tipos) são uma formas de aproximar funções ao redor de um ponto de avaliação.

In [29]:
series(f(x),x,0,6)
Out[29]:

Ainda é possível fazer de outra forma…

In [30]:
f(x).series(x,0,6)
Out[30]:

Adiantes vamos estudar graficamente a questão da aproximação das seríes de Taylor.

Tendo o SYMPY em funcionamento, é possivel fazer gráficos nele usando um comando tipo PLOT, porém o sym.plot é menos rico que o plt.plot (do MATPLOTLIB) que usamos na aula passada. Mas ele quebra bem o galho e tem uma sintax super simples;

Plot do sympy básico

In [31]:
sym.plot(cos(x))
Out[31]:
<sympy.plotting.plot.Plot at 0x5d261d0>

Note que não houve necessidade de gerar o dois vetores de números (usando NUMPY ARANGE) para gerar o gráfico (na verdade essa sintax é absurdamente simples).

In [32]:
sym.plot(cos(x),sin(x),log(x)) # vários gráficos
Out[32]:
<sympy.plotting.plot.Plot at 0x61eb810>
In [33]:
sym.plot(cos(x),log(x),(x,0.1,2*pi)) # gráficos com intervalos semelhantes definidos
Out[33]:
<sympy.plotting.plot.Plot at 0x6203490>
In [34]:

sym.plot((cos(x),(x,-pi,2*pi)),(log(x),(x,0.1,2*pi))) # gráficos com intervalos diferentes definidos
Out[34]:
<sympy.plotting.plot.Plot at 0x6478950>

Para inspecionar coisas simples, esse comando funciona muito bem, para gráficos realmente profissionais devemos usar o MATPLOTLIB. A a sintax simples é muito limitada e mesmo ações simples como mudar a cor dos gráficos já deixa deixa a simplicidade de lado;

In [35]:
p=sym.plot((cos(x),(x,-pi,2*pi)),(log(x),(x,0.1,2*pi)), show=False)
p[0].line_color = 'red'
p[1].line_color = 'green'
p.title = 'Figuras Coloridass'
p.show()
In [36]:
p=sym.plot((cos(x),(x,-pi,2*pi)),(log(x),(x,0.1,2*pi)), show=False)
p[0].line_color = (0.5,1,1) #aqui usando RGB;
p[1].line_color = lambda a : log(a) # Aqui usando uma função para mapear a cor;
p.title = 'Figuras Coloridass'
p.show()

Exercício 4: Pesquise se é possível colocar o nome dos eixos nos gráficos do sympy; Se sim, mostre um exemplo;

Quarto estágio: O que você não faria facilmente;

Vamos agora fazer uma função que faça a comparação da função e da séries;

In [37]:
def f_taylor_ordem(func, var=x, ordem=5):
    ftaylor = func.series(var,0,ordem).removeO()
    return ftaylor

Definimos uma função que faz a série de Taylor de uma função genérica (variável func), expandida na variável específica (var), que por omissão é escolhida variável x (var=x) e na ordem que quisermos (ordem), que por omissão é a ordem 5 (ordem=5);

Atente para o remove0() do final, ele serve para transformar a sérei em polinômio.

In [38]:
f_taylor_ordem(sinh(x*cos(t))) #tudo funcionando por omissão
Out[38]:
x36cos3(t)+xcos(t)
In [39]:
f_taylor_ordem(sinh(x*cos(t)),t,5)#aqui espeficicando a variável e a ordem
Out[39]:
t4(x28sinh(x)+x24cosh(x))−t2x2cosh(x)+sinh(x)

Vamos agora unir essa função a estrutura que gera gráficos do sympy;

In [40]:
def plot_taylor_compara(func, var=x, rang=(-1,1),yrang=(-1.2,1.2),ordem=[1,5]):
    xrang = (var,)+rang
    arg_plot=(func,)+tuple(map(lambda num:f_taylor_ordem(func,var,num),ordem))+(xrang,)
    display(arg_plot)
    sym.plot(*arg_plot,ylim=yrang)

Existem 4 soluções importantes usadas aqui que merecem destaque:

1 – Operações com TUPLAS: Note que a forma de definir o intervalo do plot é uma tupla, fazermos uma operação de soma com duas tuplas para gerar o formato esperado pelo plot;

2 – Uso do MAP, map mapeia em uma função os elementos listados;

map(f(x),[1,2,3]) = [f(1),f(2),f(3)]

Mas o map mapeia em uma função e devemos avisar para ele qual é a variável que ele deve mapear, para isso usamos…

3 – Função LAMBDA, essa função basicamente deixa aberto o espaço da variável declarada na função em que está aplicada. O map mapeará nesse espaço aberto;

O MAP + LAMBDA poderiam ser substituidos por um loop, mas em geral o loop é mais lento que o MAP + LAMBDA;

4 – Usar o “” para aplicar uma tupla como argumento de uma função do python; Apenas substituir a tupla não funciona, umando nome_da_tupla, os n elementos são aplicados um após o outro como argumentos da função;

Exercício 5: Como verificamos que a forma de definir intervalo no sym.plot é uma tupla?

Exercício 6: Faça o loop correspondente ao MAP abaixo;

In [41]:
map(lambda x:sin(pi*x/8),range(10))

Vamos testar nossa criação: plot_taylor_compara

In [42]:
plot_taylor_compara(tan(x)) # com omissão dos parâmetros
(tan(x),0,x33+x,(x,−1,1))