Sympy & Scipy

%matplotlib inline

from sympy import *
init_printing()

from sympy import *
x, t, z, nu = symbols('x t z nu')
diff(sin(x)*exp(x), x)

from sympy import symbols
from sympy.plotting import plot,plot3d
from sympy.abc import *

x = symbols('x')
plot(x,(x,-5,5))

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

x+1

 plot3d((x**2 + y**2, (x, -5, 5), (y, -5, 5)),(x*y, (x, -3, 3), (y, -3, 3)))

plot3d(sin(x*10)*cos(y*5) - x*y,(x,-10,10),(y,-10,10))

import numpy as np
import scipy.integrate as spi
import matplotlib.pyplot as plt

m = 1.  # particle's mass
k = 1.  # drag coefficient
g = 9.81  # gravity acceleration

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

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

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]

m = 1.  # particle's mass
b = 1.  # drag coefficient
k= 1. #spring constant

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

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

t = np.linspace(0., 10., 300) 
v = spi.odeint(f, v0, t, args=(4,2))

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

plt.xlabel('$x$')
plt.ylabel('$\dot{x}$')
plt.plot(v[:,0],v[:,1], '-',
             mew=1, ms=8, mec='w')

from random import random
from math import cos,sin

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)

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()

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

t = np.linspace(0., 30., 300) 
v = spi.odeint(ff, v0, t, args=(4,2,1,1))

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

plt.xlabel('$x$')
plt.ylabel('$\dot{x}$')
plt.plot(v[:,0],v[:,1], '-',
             mew=1, ms=8, mec='w')

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)

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

plt.plot(Amax[:,0],Amax[:,1])
plt.plot(Amax[:,0],-Amax[:,2])
plt.show

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