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]:
In [60]:
type(a)
Out[60]:
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]
In [63]:
#Ahora hago lo mismo con Numpy
a = np.arange(1000)
In [64]:
%%timeit
a**2
In [65]:
print 'Numpy es '
print 111/5.4
print 'veces mas rapido'
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
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]:
In [69]:
# Array de 1 a 9 con paso 2
b = np.arange(1, 9, 2)
b
Out[69]:
In [70]:
#Como en Matlab
a = np.ones((3,3))
a
Out[70]:
In [71]:
b = np.zeros((5,5))
b
Out[71]:
In [72]:
c = np.eye(10)
c
Out[72]:
In [73]:
d = np.diag(np.array([1, 2, 3, 4]))
d
Out[73]:
In [74]:
#Complex numbers
e = np.array([1+2j, 3+4j, 5+6*1j])
e
Out[74]:
In [75]:
#boolean
e = np.array([True, False, False, True])
e
Out[75]:
In [76]:
#String
f = np.array(['Bonjour', 'Hello', 'Hallo',])
f
Out[76]:
In [77]:
a = np.arange(10)
a
Out[77]:
In [78]:
#creo una tupla
a[0],a[2],a[-1]
Out[78]:
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]:
In [80]:
a[::4] # solo cambiamos el paso
Out[80]:
In [81]:
np.random.seed(3)
a = np.random.random_integers(0, 20, 15)
a
Out[81]:
In [82]:
(a % 3 == 0)
Out[82]:
In [83]:
mask = (a % 3 == 0)
a_multiplos_3 = a[mask]
a_multiplos_3
Out[83]:
In [84]:
#Puedo indexar y asignar al mismo tiempo
a[a % 3 == 0] = -1
a
Out[84]:
In [85]:
a = np.array([1, 2, 3, 4])
a
Out[85]:
In [86]:
#todas las operaciones aritmeticas funcionan
a+1
Out[86]:
In [87]:
j = np.arange(5)
2**(j + 1) - j
Out[87]:
In [88]:
#Multiplicacion ES ELEMENTO A ELEMENTO
a * a
Out[88]:
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]:
In [91]:
#cado objeto ndarray tiene muuchos metodos eje
c.sum()
Out[91]:
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
In [93]:
!f2py -c -m hellofortran hellofortran.f
In [94]:
%%file hello.py
import hellofortran
hellofortran.hellofortran(5)
In [95]:
# corremos el script
!python hello.py
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
Compilamos directamente a un modulo de python
In [98]:
!f2py -c dcumsum.f -m dcumsum
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]:
In [102]:
dcumsum.dcumsum(a)
Out[102]:
In [103]:
a = np.random.rand(10000)
In [104]:
%%timeit
py_dcumsum(a)
In [105]:
%%timeit
dcumsum.dcumsum(a)
In [106]:
%%timeit
a.cumsum()
In [107]:
%run srinivasan_pruebas.py
In [108]:
%run srinivasan_pruebas_vec.py
In [109]:
#para que los graficos queden empotrados
%pylab inline
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]:
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]:
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:
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]:
In [ ]: