In [2]:
#Importamos el modulo numpy con el alias np
import numpy as np
In [3]:
#Creo un array
a = np.array([1,0,0])
a
Out[3]:
In [4]:
type(a)
Out[4]:
In [5]:
#Ejemplo creo una lista de Python de 0 a 1000 y calculo el cuadrado de cada elemento
L = range(1000)
In [5]:
%%timeit
[i**2 for i in L]
In [6]:
#Ahora hago lo mismo con Numpy
a = np.arange(1000)
In [7]:
%%timeit
a**2
In [8]:
print 'Numpy es '
print 111/5.4
print 'veces mas rapido'
In [9]:
#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 [10]:
print a.shape
print b.shape
print c.shape
In [11]:
#Podemos crear arrays predeterminados con funciones muy utiles
a = np.arange(10) # un array de enteros 0 a 10(no lo incluye)
a
Out[11]:
In [12]:
# Array de 1 a 9 con paso 2
b = np.arange(1, 9, 2)
b
Out[12]:
In [13]:
#Como en Matlab
a = np.ones((3,3))
a
Out[13]:
In [14]:
b = np.zeros((5,5))
b
Out[14]:
In [15]:
c = np.eye(10)
c
Out[15]:
In [16]:
d = np.diag(np.array([1, 2, 3, 4]))
d
Out[16]:
In [17]:
#Complex numbers
e = np.array([1+2j, 3+4j, 5+6*1j])
e
Out[17]:
In [41]:
#boolean
e = np.array([True, False, False, True])
e
Out[41]:
In [42]:
#String
f = np.array(['Bonjour', 'Hello', 'Hallo',])
f
Out[42]:
In [44]:
a = np.arange(10)
a
Out[44]:
In [45]:
#creo una tupla
a[0],a[2],a[-1]
Out[45]:
In [48]:
# 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[48]:
In [51]:
a[::4] # solo cambiamos el paso
Out[51]:
In [53]:
np.random.seed(3)
a = np.random.random_integers(0, 20, 15)
a
Out[53]:
In [54]:
(a % 3 == 0)
Out[54]:
In [58]:
mask = (a % 3 == 0)
a_multiplos_3 = a[mask]
a_multiplos_3
Out[58]:
In [59]:
#Puedo indexar y asignar al mismo tiempo
a[a % 3 == 0] = -1
a
Out[59]:
In [60]:
a = np.array([1, 2, 3, 4])
a
Out[60]:
In [62]:
#todas las operaciones aritmeticas funcionan
a+1
Out[62]:
In [63]:
j = np.arange(5)
2**(j + 1) - j
Out[63]:
In [65]:
#Multiplicacion ES ELEMENTO A ELEMENTO
a * a
Out[65]:
In [68]:
#Para hacer multiplicaciones matriciales usamos dot
b = np.random.rand(3,3)
c = np.random.rand(3,3)
In [69]:
np.dot(b,c)
Out[69]:
In [43]:
#cado objeto ndarray tiene muuchos metodos eje
c.sum()
Out[43]:
In [22]:
%%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 [23]:
!f2py -c -m hellofortran hellofortran.f
In [24]:
%%file hello.py
import hellofortran
hellofortran.hellofortran(5)
In [25]:
# 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 [34]:
# 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 [27]:
%%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 [29]:
!f2py -c dcumsum.f -m dcumsum
In [30]:
#importamos el modulo recien creado
import dcumsum
In [32]:
a = np.array([1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0])
In [35]:
py_dcumsum(a)
Out[35]:
In [36]:
dcumsum.dcumsum(a)
Out[36]:
In [39]:
a = np.random.rand(10000)
In [40]:
%%timeit
py_dcumsum(a)
In [41]:
%%timeit
dcumsum.dcumsum(a)
In [42]:
%%timeit
a.cumsum()
In [6]:
%run srinivasan_pruebas.py
In [7]:
%run srinivasan_pruebas_vec.py
In [92]:
#para que los graficos queden empotrados
%pylab inline
In [93]:
X = np.linspace(-np.pi, np.pi, 256, endpoint=True)
C, S = np.cos(X), np.sin(X)
In [94]:
plot(X, C)
plot(X, S)
Out[94]:
In [95]:
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[95]:
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 [83]:
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 [84]:
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[84]:
In [8]:
%run MagicCube/code/cube_interactive.py 5
In [ ]: