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