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
%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.
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;
from sympy import __version__ __version__
'0.7.4.1'
Primeiro estágio: Operações básicas e considerações iniciais
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;
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
3*pi/2 + sym.exp(I*x) / (x**2 + y)
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:
init_printing()
Vamos repetir o nosso teste;
3*pi/2 + sym.exp(I*x) / (x**2 + y) #Muito mais bonito!
Algumas ações que podem ser feitas com expressões simbólicas: expandir (expand), simplificar(simplify) e resolver(solve);
eq = ((x+y)**2 * (x+1)) eq
expand(eq)
a = 1/x + (x*sin(x) - 1)/x print(a) #assim fica feito display(a) #assim fica bonito
simplify(a)
eq2 = Eq(x**3 - 3*x**2 + 4*x + 8., 0) display(eq2)
Notem que a equação gerada usa o comando “Eq”, para definir a igualdade;
solve(eq2,x)
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:
w0=symbols("w_0",real=True)
Atentem para a Flag “real”, isso significa que o símbolo w0 é REAL;
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.
b=symbols("b",real=True)
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))
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;
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))
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))
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))
E se houvesse uma força externa, por exemplo senoidal!
wp, A=symbols("w_p A",real=True)
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;
Derivative(cos(m*x**k)**10,x)
diff(cos(m*x**k)**10,x)
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;
Derivative(cos(y*x**k),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;
Integral(sin(k*x)*sin(m*x),(x,-pi/2,pi/2))
integrate(sin(k*x)*sin(m*x),(x,-pi/2,pi/2))
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.
series(f(x),x,0,6)
Ainda é possível fazer de outra forma…
f(x).series(x,0,6)
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
sym.plot(cos(x))
<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).
sym.plot(cos(x),sin(x),log(x)) # vários gráficos
<sympy.plotting.plot.Plot at 0x61eb810>
sym.plot(cos(x),log(x),(x,0.1,2*pi)) # gráficos com intervalos semelhantes definidos
<sympy.plotting.plot.Plot at 0x6203490>
sym.plot((cos(x),(x,-pi,2*pi)),(log(x),(x,0.1,2*pi))) # gráficos com intervalos diferentes definidos
<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;
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()
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;
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.
f_taylor_ordem(sinh(x*cos(t))) #tudo funcionando por omissão
f_taylor_ordem(sinh(x*cos(t)),t,5)#aqui espeficicando a variável e a ordem
Vamos agora unir essa função a estrutura que gera gráficos do sympy;
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;
map(lambda x:sin(pi*x/8),range(10))
Vamos testar nossa criação: plot_taylor_compara
plot_taylor_compara(tan(x)) # com omissão dos parâmetros
plot_taylor_compara(tan(x),x,(-1.5,1.5),(-3,3),[1,3,5,7,9]) # especificando tudo
Fica uma confusão para ler essses gráficos;
Vamos incluir cores:
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…
plot_taylor_compara_color(exp(x)) # isso não ajudou muito!
A função exp(x) sempre é azul!
plot_taylor_compara_color(exp(x),x,(-8,20),(-10,200),range(1,8,1)) # especificando tudo
A função exp(x) sempre é azul!
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!
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;
sym.plot_implicit(Eq(x**2 + y**2, 4))
<sympy.plotting.plot.Plot at 0x69da0d0>
plot_implicit(y > x**2 ,(x,-1,1)) # podemos pintar áreas
<sympy.plotting.plot.Plot at 0x79d1390>
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
<sympy.plotting.plot.Plot at 0x6729750>
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
<sympy.plotting.plot.Plot at 0x7d442d0>
Gráficos paramétricos
from sympy.plotting import plot_parametric
plot_parametric(cos(x), sin(x), (x, 0, 2*pi))
<sympy.plotting.plot.Plot at 0x79d1110>
plot_parametric(x*cos(x), x*sin(x), (x, 0, 10*pi))
<sympy.plotting.plot.Plot at 0x6217390>
plot_parametric(cos(x)*cos(x), cos(x)*sin(x), (x, 0, 2*pi))
<sympy.plotting.plot.Plot at 0x7a6bf90>
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!
<sympy.plotting.plot.Plot at 0xb155e50>
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()