Python – Aula 3

Publicado em por

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))
In [43]:

plot_taylor_compara(tan(x),x,(-1.5,1.5),(-3,3),[1,3,5,7,9]) # especificando tudo
(tan(x),0,x,x33+x,2×515+x33+x,17×7315+2×515+x33+x,(x,−1.5,1.5))

Fica uma confusão para ler essses gráficos;

Vamos incluir cores:

In [44]:
def plot_taylor_compara_color(func, var=x, rang=(-3.14,3.14),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,)
    gr=sym.plot(*arg_plot,ylim=yrang,title= "Grafico "+string0(func) + "e suas expansoes nas ordens "+string0(ordem),show=False)
    for i in range(1,len(ordem)+1):
        rgb=((1./random_integers(1,15)),(1./random_integers(1,15)),(1./random_integers(5,15)))
        gr[i].line_color = rgb  
    gr.show()
    print 'A função ' + string0(func) + ' sempre é azul!'

Exercício 7: Qual a serventia do random_integers?

Vamos ver o resultado do nosso trabalho…

In [45]:
plot_taylor_compara_color(exp(x)) # isso não ajudou muito!
A função exp(x) sempre é azul!
In [46]:
plot_taylor_compara_color(exp(x),x,(-8,20),(-10,200),range(1,8,1)) # especificando tudo
A função exp(x) sempre é azul!
In [47]:
plot_taylor_compara_color(exp(x),x,(-8,20),(155,165),range(1,8,1)) # Vamos mexer no yrange para ver melhor o que ocorre x=5;
A função exp(x) sempre é azul!
In [48]:
plot_taylor_compara_color(log(x+1),x,(0,2),(-0.75,1.5),range(2,15,1)) # Só está aqui pq é legal!
A função log(x + 1) sempre é azul!

Um pouco mais sobre gráficos do sympy

Como dissemos antes, para gráficos realmente profissionais devemos usar o MATPLOTLIB.

Mas mesmo com suas limitações, os gráficos do sympy são bastante versátis e com eles é possível esboçar muita coisa;

Gráficos implícitos;

In [49]:
sym.plot_implicit(Eq(x**2 + y**2, 4))
Out[49]:
<sympy.plotting.plot.Plot at 0x69da0d0>
In [50]:
plot_implicit(y > x**2 ,(x,-1,1)) # podemos pintar áreas
Out[50]:
<sympy.plotting.plot.Plot at 0x79d1390>
In [51]:
plot_implicit(And(y > x**2 , y < sqrt(x)),(x,-1,1.5),(y,-1,2),title='Intersecao de duas curvas') # podemos fazer algumas operações lógicas
Out[51]:
<sympy.plotting.plot.Plot at 0x6729750>
In [52]:
plot_implicit(Or(And(y > x**2 , y < sqrt(x)),And(y < -x**2 , y > -sqrt(-x))),(x,-1,1.5),(y,-1,2)) 
# Um pouco mais de operações lógicas
Out[52]:
<sympy.plotting.plot.Plot at 0x7d442d0>

Gráficos paramétricos

In [53]:
from sympy.plotting import plot_parametric
In [54]:
plot_parametric(cos(x), sin(x), (x, 0, 2*pi))
Out[54]:
<sympy.plotting.plot.Plot at 0x79d1110>
In [55]:
plot_parametric(x*cos(x), x*sin(x), (x, 0, 10*pi))
Out[55]:
<sympy.plotting.plot.Plot at 0x6217390>
In [56]:
plot_parametric(cos(x)*cos(x), cos(x)*sin(x), (x, 0, 2*pi))
Out[56]:
<sympy.plotting.plot.Plot at 0x7a6bf90>
In [57]:
plot_parametric(cos(20*x)*cos(x), cos(10*x)*sin(x), (x, 0, 2*pi),axis=False) # Até um pouco de arte da para fazer!
Out[57]:
<sympy.plotting.plot.Plot at 0xb155e50>
In [58]:
pl=plot_parametric(cos(20*x)*cos(x), cos(10*x)*sin(x), (x, 0, 2*pi),axis=False,show=False)
pl[0].line_color= lambda a: log(0.001*a)
pl.show()

Exercício 8 : Descubra se sympy faz gráficos 3D, em caso afirmativo mostre exemplos.