CoCalc Public Filessympy2.ipynbOpen with one click!
Author: Ney Lemke
Views : 48
Description: Jupyter notebook sympy2.ipynb
Compute Environment: Ubuntu 18.04 (Deprecated)

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

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$')
<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) plt.suptitle('Osciladores Amortecidos',size=24) plt.show()