CoCalc Shared Filessympy2.ipynb
Authors: Ney Lemke, Taina Reis
Views : 6
Description: Jupyter notebook sympy2.ipynb

# Sympy & Scipy

In [2]:
%matplotlib inline

In [3]:
from sympy import *
init_printing()


## Plotting

In [4]:
from sympy import *
x, t, z, nu = symbols('x t z nu')
diff(sin(x)*exp(x), x)


In [5]:
from sympy import symbols
from sympy.plotting import plot,plot3d
from sympy.abc import *


In [6]:
x = symbols('x')
plot(x,(x,-5,5))

<sympy.plotting.plot.Plot at 0x7fa33b0c6cf8>
In [7]:
from sympy import symbols
from sympy.plotting import plot
p1 = plot(x*x,xlabel='teste',line_color=0.0)
p2 = plot(x)
p1.extend(p2)
p1

<sympy.plotting.plot.Plot at 0x7fa369101860>
In [8]:
x+1


In [9]:
 plot3d((x**2 + y**2, (x, -5, 5), (y, -5, 5)),(x*y, (x, -3, 3), (y, -3, 3)))

<sympy.plotting.plot.Plot at 0x7fa368795d68>
In [10]:
plot3d(sin(x*10)*cos(y*5) - x*y,(x,-10,10),(y,-10,10))

<sympy.plotting.plot.Plot at 0x7fa33b16db70>

## ODE

In [12]:
import numpy as np
import scipy.integrate as spi
import matplotlib.pyplot as plt


In [13]:
m = 1.  # particle's mass
k = 1.  # drag coefficient
g = 9.81  # gravity acceleration

In [14]:
# We want to evaluate the system on 30 linearly
# spaced times between t=0 and t=3.
t = np.linspace(0., 3., 30)
# We simulate the system for different values of k.
for k in np.linspace(0., 1., 5):
# We simulate the system and evaluate $v$ on
# the given times.
v = spi.odeint(f, v0, t, args=(k,))
# We plot the particle's trajectory.
plt.plot(v[:,0], v[:,1], 'o-',
mew=1, ms=8, mec='w',
label='k={0:.1f}'.format(k))
plt.legend()
plt.xlim(0, 12)
plt.ylim(0, 6)

--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-14-d5648b085ab5> in <module>() 6 # We simulate the system and evaluate $v$ on 7 # the given times. ----> 8 v = spi.odeint(f, v0, t, args=(k,)) 9 # We plot the particle's trajectory. 10 plt.plot(v[:,0], v[:,1], 'o-', NameError: name 'v0' is not defined 
In [16]:
# The initial position is (0, 0).
v0 = np.zeros(4)
# The initial speed vector is oriented
# to the top right.
v0[2] = 4.
v0[3] = 10.

In [ ]:
def f(v, t0, k):
# v has four components: v=[u, u'].
u, udot = v[:2], v[2:]
# We compute the second derivative u'' of u.
udotdot = -k/m * udot
udotdot[1] -= g
# We return v'=[u', u''].
return np.r_[udot, udotdot]


In [17]:
m = 1.  # particle's mass
b = 1.  # drag coefficient
k= 1. #spring constant

In [19]:
# The initial position is (0, 0).
v0 = np.zeros(2)
# The initial speed vector is oriented
# to the top right.
v0[0] = 1.
v0[1] = 0.

In [20]:
def f(v, t0, *args):
# v has four components: v=[u, u'].
b,k=args
u, udot = v[0], v[1]
# We compute the second derivative u'' of u.
udotdot = -b/m*udot-k/m*u
# We return v'=[u', u''].
return np.array([udot,udotdot])

In [21]:
t = np.linspace(0., 10., 300)
v = spi.odeint(f, v0, t, args=(4,2))

In [22]:
plt.plot(t,v[:,0], '-',
mew=1, ms=8, mec='w')
plt.ylabel('$x$')
plt.xlabel('$t$')

<matplotlib.text.Text at 0x7fa330160a20>
In [23]:
plt.xlabel('$x$')
plt.ylabel('$\dot{x}$')
plt.plot(v[:,0],v[:,1], '-',
mew=1, ms=8, mec='w')


[<matplotlib.lines.Line2D at 0x7fa3300d1320>]

## VÃ¡rios casos para a mesma condiÃ§Ã£o inicial

In [24]:
from random import random
from math import cos,sin

In [14]:
t = np.linspace(0., 30., 300)
for theta in np.linspace(0, 6.28, 100):
v0=np.array([float(3*cos(theta)),float(3*sin(theta))])
# We simulate the system and evaluate $v$ on
# the given times.
v = spi.odeint(f, v0, t, args=(2,1))
# We plot the particle's trajectory.
plt.plot(v[:,0], v[:,1], '-',mew=1, ms=8, mec='w')
plt.legend()
plt.xlabel('$x$')
plt.ylabel('$\dot{x}$')
plt.title('$b=2$')
plt.xlim(-2, 2)
plt.ylim(-2, 2)

/Users/neylemke/.pyenv/versions/anaconda3-2.4.0/lib/python3.5/site-packages/matplotlib/axes/_axes.py:545: UserWarning: No labelled objects found. Use label='...' kwarg on individual plots. warnings.warn("No labelled objects found. "
(-2, 2)

## Diagrama de Fase para vÃ¡rios parÃ¢metros

In [25]:
t = np.linspace(0., 30., 300)
nmax=6
fig, axarr = plt.subplots(nmax, nmax,figsize=(16,12))
for i in range(nmax):
for j in range(nmax):
for theta in np.linspace(0, 6.28, 100):
v0=np.array([float(3*cos(theta)),float(3*sin(theta))])
# We simulate the system and evaluate $v$ on
# the given times.
v = spi.odeint(f, v0, t, args=(i+1,nmax-j))
# We plot the particle's trajectory.
axarr[i,j].plot(v[:,0], v[:,1], '-',mew=1, ms=8, mec='w')
if (4*(nmax-j)<=(i+1)**2):
color='yellow'
else:
color='orange'
axarr[nmax-i-1,nmax-j-1].set_facecolor(color)
axarr[nmax-i-1,nmax-j-1].set_xlim(-2,2)
axarr[nmax-i-1,nmax-j-1].set_ylim(-2,2)
axarr[nmax-i-1,nmax-j-1].set_title('$b$={:2.1f} $k$={:2.1f}'
.format((i+1)*1.0,(nmax-j)*1.0))
plt.setp([axarr[i, j].get_xticklabels() for i in range(nmax) for j in range(nmax)], visible=False)
plt.setp([axarr[i, j].get_yticklabels() for i in range(nmax) for j in range(nmax)], visible=False)
axarr[nmax-1,2].set_xlabel('                   Atrito',size=24)
axarr[3,0].set_ylabel('           Constante da Mola',size=24)