CoCalc Public Filessympy2.html
Author: Ney Lemke
Views : 52
Description: Jupyter html version of sympy2.ipynb
Compute Environment: Ubuntu 18.04 (Deprecated)
(File too big to render with math typesetting.)
sympy2

# 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]


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)