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()
: