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)
Out[4]:
$$e^{x} \sin{\left (x \right )} + e^{x} \cos{\left (x \right )}$$
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))
Out[6]:
<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
Out[7]:
<sympy.plotting.plot.Plot at 0x7fa369101860>
In [8]:
x+1
Out[8]:
$$x + 1$$
In [9]:
 plot3d((x**2 + y**2, (x, -5, 5), (y, -5, 5)),(x*y, (x, -3, 3), (y, -3, 3)))
Out[9]:
<sympy.plotting.plot.Plot at 0x7fa368795d68>
In [10]:
plot3d(sin(x*10)*cos(y*5) - x*y,(x,-10,10),(y,-10,10))
Out[10]:
<sympy.plotting.plot.Plot at 0x7fa33b16db70>

ODE

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

Lançamento de Projétil

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]

Oscilador Amortecido

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$')
Out[22]:
<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')
Out[23]:
[<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. "
Out[14]:
(-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)
plt.suptitle('Osciladores Amortecidos',size=24)
plt.show()

Oscilações Forçadas

In [26]:
def ff(v, t0, *args):
    # v has four components: v=[u, u'].
    b,k,A,omega=args
    u, udot = v[0], v[1]
    # We compute the second derivative u'' of u.
    udotdot = -b/m*udot-k/m*u+A*cos(omega*t0)/m
    # We return v'=[u', u''].
    return np.array([udot,udotdot])
In [27]:
t = np.linspace(0., 30., 300) 
v = spi.odeint(ff, v0, t, args=(4,2,1,1))
In [94]:
plt.plot(t,v[:,0], '-',
             mew=1, ms=8, mec='w')
plt.ylabel('$x$')
plt.xlabel('$t$')
Out[94]:
<matplotlib.text.Text at 0x12f9abb70>
In [28]:
plt.xlabel('$x$')
plt.ylabel('$\dot{x}$')
plt.plot(v[:,0],v[:,1], '-',
             mew=1, ms=8, mec='w')
Out[28]:
[<matplotlib.lines.Line2D at 0x7fa329621668>]
In [29]:
t = np.linspace(0., 300., 300)
A=2
omega=2
b=2
k=1
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(ff, v0, t, args=(b,k,A,omega))
    # We plot the particle's trajectory.
    plt.plot(v[:,0], v[:,1], '-',mew=1, ms=8, mec='w')
plt.legend()
plt.xlabel('$x$',size=16)
plt.ylabel('$\dot{x}$',size=16)
plt.title('$b$={} $k$={} $A$={} $\omega$={}'.format(b,k,A,omega))
plt.xlim(-2, 2)
plt.ylim(-2, 2)
/usr/local/lib/python3.4/dist-packages/matplotlib/axes/_axes.py:545: UserWarning: No labelled objects found. Use label='...' kwarg on individual plots.
  warnings.warn("No labelled objects found. "
Out[29]:
$$\left ( -2, \quad 2\right )$$
In [30]:
tmax=3000
steps=30000
t = np.linspace(0., tmax, steps) 
def maxinv(omega):
    v = spi.odeint(ff, v0, t, args=(0.1,4,1,omega))
    return [omega,max(v[-10000:,0]),min(-v[-10000:,0])] 
    
Amax=np.array([maxinv(omega) for omega in np.linspace(0.1, 5, num=100)])
In [35]:
plt.plot(Amax[:,0],Amax[:,1])
plt.plot(Amax[:,0],-Amax[:,2])
plt.show
Out[35]:
<function matplotlib.pyplot.show>
In [36]:
v = spi.odeint(ff, v0, t, args=(2,5,1,0.1))
plt.plot(t[-1000:],v[-1000:,0], '-',
             mew=1, ms=8, mec='w')
plt.ylabel('$x$')
plt.xlabel('$t$')
print(max(v[-10000:,0]))
print(min(v[-10000:,0]))
0.200240045859
-0.200240016589