Numpy

Biblioteca(modulo) numérico de Python

  • Arrays multimensionales (ndarray)
  • Más cercano al hardware(numéricamente más eficiente)
  • Designado para aplicaciones científicas
  • Utiliza memoria de manera más eficiente
  • Todos los modulos cientificos lo utilizan como tipo de dato básico(Opencv a partir de Opencv2¡¡¡¡)

Under the hood: la dispocisión de memoria en un numpy array es :

bloque de memoria + esquema de indexado + descriptor del tipo de dato

  • datos planos
  • como alocar elementos
  • como interpretar elementos


In [58]:
#Importamos el modulo numpy con el alias np
import numpy as np

In [59]:
#Creo un array 
a = np.array([1,0,0])
a


Out[59]:
array([1, 0, 0])

In [60]:
type(a)


Out[60]:
numpy.ndarray

Python posee por defecto un tipo de datos que se asemeja(listas), pero es numéricamente ineficiente


In [61]:
#Ejemplo creo una lista de Python de 0 a 1000 y calculo el cuadrado de cada elemento
L = range(1000)

In [62]:
%%timeit
[i**2 for i in L]


10000 loops, best of 3: 101 µs per loop

In [63]:
#Ahora hago lo mismo con Numpy
a = np.arange(1000)

In [64]:
%%timeit
a**2


The slowest run took 16.56 times longer than the fastest. This could mean that an intermediate result is being cached 
100000 loops, best of 3: 2.31 µs per loop

In [65]:
print 'Numpy es '
print  111/5.4
print 'veces mas rapido'


Numpy es 
20.5555555556
veces mas rapido

Caracteristicas y utilidades principales:


In [66]:
#Creando arrays

a = np.array([1,1,1]) # 1D
b = np.array([[1,1,1],[1,1,1]]) #2d (Matrix)
c = np.array([[[1,1,1],[1,1,1],[1,1,1]]]) #3D (Tensor...)

In [67]:
print a.shape
print b.shape
print c.shape


(3,)
(2, 3)
(1, 3, 3)

In [68]:
#Podemos crear arrays predeterminados con funciones muy utiles

a = np.arange(10) # un array de enteros 0 a 10(no lo incluye) 
a


Out[68]:
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

In [69]:
# Array de 1 a 9 con paso 2
b = np.arange(1, 9, 2)
b


Out[69]:
array([1, 3, 5, 7])

In [70]:
#Como en Matlab
a = np.ones((3,3))
a


Out[70]:
array([[ 1.,  1.,  1.],
       [ 1.,  1.,  1.],
       [ 1.,  1.,  1.]])

In [71]:
b = np.zeros((5,5))
b


Out[71]:
array([[ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.]])

In [72]:
c = np.eye(10)
c


Out[72]:
array([[ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.]])

In [73]:
d = np.diag(np.array([1, 2, 3, 4]))
d


Out[73]:
array([[1, 0, 0, 0],
       [0, 2, 0, 0],
       [0, 0, 3, 0],
       [0, 0, 0, 4]])

In [74]:
#Complex numbers
e = np.array([1+2j, 3+4j, 5+6*1j])
e


Out[74]:
array([ 1.+2.j,  3.+4.j,  5.+6.j])

In [75]:
#boolean
e = np.array([True, False, False, True])
e


Out[75]:
array([ True, False, False,  True], dtype=bool)

In [76]:
#String
f = np.array(['Bonjour', 'Hello', 'Hallo',])
f


Out[76]:
array(['Bonjour', 'Hello', 'Hallo'], 
      dtype='|S7')

Indexing and slicing

Los items de un array pueden accederse de una manera natural como en Matlab(hay que tener en cuenta que los indices comienzan en 0)


In [77]:
a = np.arange(10)
a


Out[77]:
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

In [78]:
#creo una tupla 
a[0],a[2],a[-1]


Out[78]:
(0, 2, 9)

In [79]:
# Slicing [comienzo:final:paso]
# Los tres no son necesarios explicitamente ya que por default comienzo=0, final=[-1] y el paso=1

a[2:8:2]


Out[79]:
array([2, 4, 6])

In [80]:
a[::4] # solo cambiamos el paso


Out[80]:
array([0, 4, 8])

Fancy indexing

Hay más maneras de indexar y asignar


In [81]:
np.random.seed(3)
a = np.random.random_integers(0, 20, 15)
a


Out[81]:
array([10,  3,  8,  0, 19, 10, 11,  9, 10,  6,  0, 20, 12,  7, 14])

In [82]:
(a % 3 == 0)


Out[82]:
array([False,  True, False,  True, False, False, False,  True, False,
        True,  True, False,  True, False, False], dtype=bool)

In [83]:
mask = (a % 3 == 0)
a_multiplos_3 = a[mask]
a_multiplos_3


Out[83]:
array([ 3,  0,  9,  6,  0, 12])

In [84]:
#Puedo indexar y asignar al mismo tiempo
a[a % 3 == 0] = -1
a


Out[84]:
array([10, -1,  8, -1, 19, 10, 11, -1, 10, -1, -1, 20, -1,  7, 14])

Elementwise operations

Operar con funciones en cada elemento de un array(NO USE FOR)


In [85]:
a = np.array([1, 2, 3, 4])
a


Out[85]:
array([1, 2, 3, 4])

In [86]:
#todas las operaciones aritmeticas funcionan
a+1


Out[86]:
array([2, 3, 4, 5])

In [87]:
j = np.arange(5)
2**(j + 1) - j


Out[87]:
array([ 2,  3,  6, 13, 28])

In [88]:
#Multiplicacion ES ELEMENTO A ELEMENTO
a * a


Out[88]:
array([ 1,  4,  9, 16])

In [89]:
#Para hacer multiplicaciones matriciales usamos dot

b = np.random.rand(3,3)
c = np.random.rand(3,3)

In [90]:
np.dot(b,c)


Out[90]:
array([[ 0.75172429,  0.76559311,  0.82038214],
       [ 0.5466532 ,  0.49578046,  0.62333632],
       [ 0.46871279,  0.54809788,  0.43332846]])

In [91]:
#cado objeto ndarray tiene muuchos metodos eje
c.sum()


Out[91]:
4.4680994807263517

Optimizaciones con Fortran y C


In [92]:
%%file hellofortran.f
C File  hellofortran.f
        subroutine hellofortran (n)
        integer n
       
        do 100 i=0, n
            print *, "Hola Soy Fortran tengo muuchos años"
100     continue
        end


Overwriting hellofortran.f

Generamos un modulo de Python con f2py


In [93]:
!f2py -c -m hellofortran hellofortran.f


running build
running config_cc
unifing config_cc, config, build_clib, build_ext, build commands --compiler options
running config_fc
unifing config_fc, config, build_clib, build_ext, build commands --fcompiler options
running build_src
build_src
building extension "hellofortran" sources
f2py options: []
f2py:> /tmp/tmpYoWbrk/src.linux-x86_64-2.7/hellofortranmodule.c
creating /tmp/tmpYoWbrk/src.linux-x86_64-2.7
Reading fortran codes...
	Reading file 'hellofortran.f' (format:fix,strict)
Post-processing...
	Block: hellofortran
			Block: hellofortran
Post-processing (stage 2)...
Building modules...
	Building module "hellofortran"...
		Constructing wrapper function "hellofortran"...
		  hellofortran(n)
	Wrote C/API module "hellofortran" to file "/tmp/tmpYoWbrk/src.linux-x86_64-2.7/hellofortranmodule.c"
  adding '/tmp/tmpYoWbrk/src.linux-x86_64-2.7/fortranobject.c' to sources.
  adding '/tmp/tmpYoWbrk/src.linux-x86_64-2.7' to include_dirs.
copying /home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/f2py/src/fortranobject.c -> /tmp/tmpYoWbrk/src.linux-x86_64-2.7
copying /home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/f2py/src/fortranobject.h -> /tmp/tmpYoWbrk/src.linux-x86_64-2.7
build_src: building npy-pkg config files
running build_ext
customize UnixCCompiler
customize UnixCCompiler using build_ext
customize Gnu95FCompiler
Found executable /usr/bin/gfortran
customize Gnu95FCompiler
customize Gnu95FCompiler using build_ext
building 'hellofortran' extension
compiling C sources
C compiler: gcc -pthread -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC

creating /tmp/tmpYoWbrk/tmp
creating /tmp/tmpYoWbrk/tmp/tmpYoWbrk
creating /tmp/tmpYoWbrk/tmp/tmpYoWbrk/src.linux-x86_64-2.7
compile options: '-I/tmp/tmpYoWbrk/src.linux-x86_64-2.7 -I/home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/core/include -I/home/elsuizo/anaconda/include/python2.7 -c'
gcc: /tmp/tmpYoWbrk/src.linux-x86_64-2.7/hellofortranmodule.c
In file included from /home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/ndarraytypes.h:1804:0,
                 from /home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/ndarrayobject.h:17,
                 from /home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/arrayobject.h:4,
                 from /tmp/tmpYoWbrk/src.linux-x86_64-2.7/fortranobject.h:13,
                 from /tmp/tmpYoWbrk/src.linux-x86_64-2.7/hellofortranmodule.c:17:
/home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: #warning "Using deprecated NumPy API, disable it by " "#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp]
 #warning "Using deprecated NumPy API, disable it by " \
  ^
gcc: /tmp/tmpYoWbrk/src.linux-x86_64-2.7/fortranobject.c
In file included from /home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/ndarraytypes.h:1804:0,
                 from /home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/ndarrayobject.h:17,
                 from /home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/arrayobject.h:4,
                 from /tmp/tmpYoWbrk/src.linux-x86_64-2.7/fortranobject.h:13,
                 from /tmp/tmpYoWbrk/src.linux-x86_64-2.7/fortranobject.c:2:
/home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: #warning "Using deprecated NumPy API, disable it by " "#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp]
 #warning "Using deprecated NumPy API, disable it by " \
  ^
compiling Fortran sources
Fortran f77 compiler: /usr/bin/gfortran -Wall -g -ffixed-form -fno-second-underscore -fPIC -O3 -funroll-loops
Fortran f90 compiler: /usr/bin/gfortran -Wall -g -fno-second-underscore -fPIC -O3 -funroll-loops
Fortran fix compiler: /usr/bin/gfortran -Wall -g -ffixed-form -fno-second-underscore -Wall -g -fno-second-underscore -fPIC -O3 -funroll-loops
compile options: '-I/tmp/tmpYoWbrk/src.linux-x86_64-2.7 -I/home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/core/include -I/home/elsuizo/anaconda/include/python2.7 -c'
gfortran:f77: hellofortran.f
/usr/bin/gfortran -Wall -g -Wall -g -shared /tmp/tmpYoWbrk/tmp/tmpYoWbrk/src.linux-x86_64-2.7/hellofortranmodule.o /tmp/tmpYoWbrk/tmp/tmpYoWbrk/src.linux-x86_64-2.7/fortranobject.o /tmp/tmpYoWbrk/hellofortran.o -L/home/elsuizo/anaconda/lib -lpython2.7 -lgfortran -o ./hellofortran.so
Removing build directory /tmp/tmpYoWbrk

Importamos el modulo que generamos y lo utilizamos


In [94]:
%%file hello.py
import hellofortran

hellofortran.hellofortran(5)


Overwriting hello.py

In [95]:
# corremos el script
!python hello.py


 Hola Soy Fortran tengo muuchos años
 Hola Soy Fortran tengo muuchos años
 Hola Soy Fortran tengo muuchos años
 Hola Soy Fortran tengo muuchos años
 Hola Soy Fortran tengo muuchos años
 Hola Soy Fortran tengo muuchos años

Ejemplo 2: suma acumulativa, entrada vector, salida vector

Primero hacemos la implementacion en Python puro, este tipo de suma es particularmente costoso ya que se necesitan loops for


In [96]:
# Esta no es la mejor implementacion
#Porque el loop esta implementado en Python
def py_dcumsum(a):
    b = np.empty_like(a)
    b[0] = a[0]
    for n in range(1,len(a)):
        b[n] = b[n-1]+a[n]
    return b

Ahora hacemos la implementacion en Fortran


In [97]:
%%file dcumsum.f
c File dcumsum.f
       subroutine dcumsum(a, b, n)
       double precision a(n)
       double precision b(n)
       integer n
cf2py  intent(in) :: a
cf2py  intent(out) :: b
cf2py  intent(hide) :: n

       b(1) = a(1)
       do 100 i=2, n
           b(i) = b(i-1) + a(i)
100    continue
       end


Overwriting dcumsum.f

Compilamos directamente a un modulo de python


In [98]:
!f2py -c dcumsum.f -m dcumsum


running build
running config_cc
unifing config_cc, config, build_clib, build_ext, build commands --compiler options
running config_fc
unifing config_fc, config, build_clib, build_ext, build commands --fcompiler options
running build_src
build_src
building extension "dcumsum" sources
f2py options: []
f2py:> /tmp/tmpztScmZ/src.linux-x86_64-2.7/dcumsummodule.c
creating /tmp/tmpztScmZ/src.linux-x86_64-2.7
Reading fortran codes...
	Reading file 'dcumsum.f' (format:fix,strict)
Post-processing...
	Block: dcumsum
			Block: dcumsum
Post-processing (stage 2)...
Building modules...
	Building module "dcumsum"...
		Constructing wrapper function "dcumsum"...
		  b = dcumsum(a)
	Wrote C/API module "dcumsum" to file "/tmp/tmpztScmZ/src.linux-x86_64-2.7/dcumsummodule.c"
  adding '/tmp/tmpztScmZ/src.linux-x86_64-2.7/fortranobject.c' to sources.
  adding '/tmp/tmpztScmZ/src.linux-x86_64-2.7' to include_dirs.
copying /home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/f2py/src/fortranobject.c -> /tmp/tmpztScmZ/src.linux-x86_64-2.7
copying /home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/f2py/src/fortranobject.h -> /tmp/tmpztScmZ/src.linux-x86_64-2.7
build_src: building npy-pkg config files
running build_ext
customize UnixCCompiler
customize UnixCCompiler using build_ext
customize Gnu95FCompiler
Found executable /usr/bin/gfortran
customize Gnu95FCompiler
customize Gnu95FCompiler using build_ext
building 'dcumsum' extension
compiling C sources
C compiler: gcc -pthread -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC

creating /tmp/tmpztScmZ/tmp
creating /tmp/tmpztScmZ/tmp/tmpztScmZ
creating /tmp/tmpztScmZ/tmp/tmpztScmZ/src.linux-x86_64-2.7
compile options: '-I/tmp/tmpztScmZ/src.linux-x86_64-2.7 -I/home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/core/include -I/home/elsuizo/anaconda/include/python2.7 -c'
gcc: /tmp/tmpztScmZ/src.linux-x86_64-2.7/dcumsummodule.c
In file included from /home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/ndarraytypes.h:1804:0,
                 from /home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/ndarrayobject.h:17,
                 from /home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/arrayobject.h:4,
                 from /tmp/tmpztScmZ/src.linux-x86_64-2.7/fortranobject.h:13,
                 from /tmp/tmpztScmZ/src.linux-x86_64-2.7/dcumsummodule.c:18:
/home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: #warning "Using deprecated NumPy API, disable it by " "#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp]
 #warning "Using deprecated NumPy API, disable it by " \
  ^
/tmp/tmpztScmZ/src.linux-x86_64-2.7/dcumsummodule.c:111:12: warning: ‘f2py_size’ defined but not used [-Wunused-function]
 static int f2py_size(PyArrayObject* var, ...)
            ^
gcc: /tmp/tmpztScmZ/src.linux-x86_64-2.7/fortranobject.c
In file included from /home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/ndarraytypes.h:1804:0,
                 from /home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/ndarrayobject.h:17,
                 from /home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/arrayobject.h:4,
                 from /tmp/tmpztScmZ/src.linux-x86_64-2.7/fortranobject.h:13,
                 from /tmp/tmpztScmZ/src.linux-x86_64-2.7/fortranobject.c:2:
/home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: #warning "Using deprecated NumPy API, disable it by " "#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp]
 #warning "Using deprecated NumPy API, disable it by " \
  ^
compiling Fortran sources
Fortran f77 compiler: /usr/bin/gfortran -Wall -g -ffixed-form -fno-second-underscore -fPIC -O3 -funroll-loops
Fortran f90 compiler: /usr/bin/gfortran -Wall -g -fno-second-underscore -fPIC -O3 -funroll-loops
Fortran fix compiler: /usr/bin/gfortran -Wall -g -ffixed-form -fno-second-underscore -Wall -g -fno-second-underscore -fPIC -O3 -funroll-loops
compile options: '-I/tmp/tmpztScmZ/src.linux-x86_64-2.7 -I/home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/core/include -I/home/elsuizo/anaconda/include/python2.7 -c'
gfortran:f77: dcumsum.f
/usr/bin/gfortran -Wall -g -Wall -g -shared /tmp/tmpztScmZ/tmp/tmpztScmZ/src.linux-x86_64-2.7/dcumsummodule.o /tmp/tmpztScmZ/tmp/tmpztScmZ/src.linux-x86_64-2.7/fortranobject.o /tmp/tmpztScmZ/dcumsum.o -L/home/elsuizo/anaconda/lib -lpython2.7 -lgfortran -o ./dcumsum.so
Removing build directory /tmp/tmpztScmZ

In [99]:
#importamos el modulo recien creado

import dcumsum

In [100]:
a = np.array([1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0])

In [101]:
py_dcumsum(a)


Out[101]:
array([  1.,   3.,   6.,  10.,  15.,  21.,  28.,  36.])

In [102]:
dcumsum.dcumsum(a)


Out[102]:
array([  1.,   3.,   6.,  10.,  15.,  21.,  28.,  36.])

Ahora los ponemos a prueba


In [103]:
a = np.random.rand(10000)

In [104]:
%%timeit 
py_dcumsum(a)


100 loops, best of 3: 4.38 ms per loop

In [105]:
%%timeit
dcumsum.dcumsum(a)


The slowest run took 16.02 times longer than the fastest. This could mean that an intermediate result is being cached 
100000 loops, best of 3: 17.3 µs per loop

In [106]:
%%timeit
a.cumsum()


10000 loops, best of 3: 45.7 µs per loop

Guauuuu¡¡¡¡

Ejemplo de vectorizacion


In [107]:
%run srinivasan_pruebas.py


srinivasan_pruebas.py:58: RuntimeWarning: overflow encountered in square
  a_2_2 = np.sum((imagen5 - imagen4) ** 2 * mask)
srinivasan_pruebas.py:58: RuntimeWarning: invalid value encountered in multiply
  a_2_2 = np.sum((imagen5 - imagen4) ** 2 * mask)
---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
/home/elsuizo/Charlas_presentaciones/Python_scientific/srinivasan_pruebas.py in <module>()
     71 
     72 #Resolvemos el sistema
---> 73 x = np.linalg.solve(A,B)
     74 
     75 #tiempo de procesamiento

/home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/linalg/linalg.pyc in solve(a, b)
    379     signature = 'DD->D' if isComplexType(t) else 'dd->d'
    380     extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
--> 381     r = gufunc(a, b, signature=signature, extobj=extobj)
    382 
    383     return wrap(r.astype(result_t))

/home/elsuizo/anaconda/lib/python2.7/site-packages/numpy/linalg/linalg.pyc in _raise_linalgerror_singular(err, flag)
     88 
     89 def _raise_linalgerror_singular(err, flag):
---> 90     raise LinAlgError("Singular matrix")
     91 
     92 def _raise_linalgerror_nonposdef(err, flag):

LinAlgError: Singular matrix

In [108]:
%run srinivasan_pruebas_vec.py


Tiempo de procesamiento: 0.068244934082
[-0.5 -0.5]

Matplotlib

Es el modulo para realizar y visualizar graficos

Pagina oficial


In [109]:
#para que los graficos queden empotrados
%pylab inline


Populating the interactive namespace from numpy and matplotlib
WARNING: pylab import has clobbered these variables: ['f', 'arrow', 'e']
`%matplotlib` prevents importing * from pylab and numpy

In [110]:
X = np.linspace(-np.pi, np.pi, 256, endpoint=True)
C, S = np.cos(X), np.sin(X)

In [111]:
plot(X, C)
plot(X, S)


Out[111]:
[<matplotlib.lines.Line2D at 0x7f94fd2a21d0>]

In [112]:
t = 2 * np.pi / 3

plot(X, C, color="blue", linewidth=2.5, linestyle="-", label="cosine")
plot(X, S, color="red",  linewidth=2.5, linestyle="-", label="sine")
plot([t, t], [0, np.cos(t)], color='blue', linewidth=2.5, linestyle="--")
scatter([t, ], [np.cos(t), ], 50, color='blue')

annotate(r'$sin(\frac{2\pi}{3})=\frac{\sqrt{3}}{2}$',
            xy=(t, np.sin(t)), xycoords='data',
            xytext=(+10, +30), textcoords='offset points', fontsize=16,
            arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))

plot([t, t],[0, np.sin(t)], color='red', linewidth=2.5, linestyle="--")
scatter([t, ],[np.sin(t), ], 50, color='red')

annotate(r'$cos(\frac{2\pi}{3})=-\frac{1}{2}$',
            xy=(t, np.cos(t)), xycoords='data',
            xytext=(-90, -50), textcoords='offset points', fontsize=16,
            arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))


Out[112]:
<matplotlib.text.Annotation at 0x7f94f4b6ccd0>

Ejemplos

Ejemplo 1 Ecuaciones de Navier-Stokes

con convección no-lineal en 2D

Now we solve 2D Convection, represented by the pair of coupled partial differential equations below:

$$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = 0$$$$\frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} = 0$$

Discretizing these equations using the methods we've applied previously yields:

$$\frac{u_{i,j}^{n+1}-u_{i,j}^n}{\Delta t} + u_{i,j}^n \frac{u_{i,j}^n-u_{i-1,j}^n}{\Delta x} + v_{i,j}^n \frac{u_{i,j}^n-u_{i,j-1}^n}{\Delta y} = 0$$$$\frac{v_{i,j}^{n+1}-v_{i,j}^n}{\Delta t} + u_{i,j}^n \frac{v_{i,j}^n-v_{i-1,j}^n}{\Delta x} + v_{i,j}^n \frac{v_{i,j}^n-v_{i,j-1}^n}{\Delta y} = 0$$

Initial Conditions

The initial conditions are the same that we used for 1D convection, applied in both the x and y directions.

$$u,\ v\ = \begin{cases}\begin{matrix} 2 & \text{for } x,y \in (0.5, 1)\times(0.5,1) \cr 1 & \text{everywhere else} \end{matrix}\end{cases}$$

Boundary Conditions

The boundary conditions hold u and v equal to 1 along the boundaries of the grid .

$$u = 1,\ v = 1 \text{ for } \begin{cases} \begin{matrix}x=0,2\cr y=0,2 \end{matrix}\end{cases}$$

In [113]:
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np

###variable declarations
nx = 101
ny = 101
nt = 80
c = 1
dx = 2.0/(nx-1)
dy = 2.0/(ny-1)
sigma = .2
dt = sigma*dx

x = np.linspace(0,2,nx)
y = np.linspace(0,2,ny)

u = np.ones((ny,nx)) ##create a 1xn vector of 1's
v = np.ones((ny,nx))
un = np.ones((ny,nx))
vn = np.ones((ny,nx))

###Assign initial conditions

u[.5/dy:1/dy+1,.5/dx:1/dx+1]=2 ##set hat function I.C. : u(.5<=x<=1 && .5<=y<=1 ) is 2
v[.5/dy:1/dy+1,.5/dx:1/dx+1]=2 ##set hat function I.C. : u(.5<=x<=1 && .5<=y<=1 ) is 2

for n in range(nt+1): ##loop across number of time steps
    un[:] = u[:]
    vn[:] = v[:]

    u[1:,1:]=un[1:,1:]-(un[1:,1:]*dt/dx*(un[1:,1:]-un[0:-1,1:]))-vn[1:,1:]*dt/dy*(un[1:,1:]-un[1:,0:-1]) 
    v[1:,1:]=vn[1:,1:]-(un[1:,1:]*dt/dx*(vn[1:,1:]-vn[0:-1,1:]))-vn[1:,1:]*dt/dy*(vn[1:,1:]-vn[1:,0:-1])
    
    u[0,:] = 1
    u[-1,:] = 1
    u[:,0] = 1
    u[:,-1] = 1
    
    v[0,:] = 1
    v[-1,:] = 1
    v[:,0] = 1
    v[:,-1] = 1

In [114]:
from matplotlib import cm ##cm = "colormap" for changing the 3d plot color palette
fig = plt.figure(figsize=(11,7), dpi=100)
ax = fig.gca(projection='3d')
X,Y = np.meshgrid(x,y)

ax.plot_surface(X,Y,u, cmap=cm.coolwarm)


Out[114]:
<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x7f94fd2a0790>

In [ ]: