Deberías de haber completado tu propio código hasta el Paso 5 antes de continuar con esta lección. Al igual que con los pasos 1 a 4, vamos a avanzar gradualmente, por lo que es importante completar el paso anterior!
Seguimos...
Ahora resolveremos un problema de convección 2D, representado por el par de ecuaciones diferenciales parciales acopladas que se presentan a continuación:
$$\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$$Discretizando las ecuaciones con los métodos que hemos aplicado anteriormente se obtiene:
Reorganización de las dos ecuaciones, se resuelve por $u_{i,j}^{n+1}$ and $v_{i,j}^{n+1}$, respectivamente. Ten en cuenta que estas ecuaciones también están acopladas.
Las condiciones iniciales son las mismos que se utilizaron para la convección en 1D, aplicadas tanto en la dirección x como en la dirección y.
$$u,\ v\ = \begin{cases}\begin{matrix} 2 & \text{para } x,y \in (0.5, 1)\times(0.5,1) \cr 1 & \text{en el resto de puntos} \end{matrix}\end{cases}$$Las condiciones de contorno son u y v igual a 1 a lo largo de los límites de la malla (cuadrícula).
$$u = 1,\ v = 1 \text{ para } \begin{cases} \begin{matrix}x=0,2\cr y=0,2 \end{matrix}\end{cases}$$%pylab inline
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
### Declaración de variables
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)) ## crea un vector 1xn de unos
v = np.ones((ny,nx))
un = np.ones((ny,nx))
vn = np.ones((ny,nx))
### Asignar variables inciales
u[.5/dy:1/dy+1,.5/dx:1/dx+1]=2 ## establece función de sobrero como C.I.: u(.5<=x<=1 && .5<=y<=1 ) is 2
v[.5/dy:1/dy+1,.5/dx:1/dx+1]=2 ##establece función de sobrero como C.I: u(.5<=x<=1 && .5<=y<=1 ) is 2
for n in range(nt+1): ## Bucle para los incrementos de tiempo
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
Populating the interactive namespace from numpy and matplotlib
from matplotlib import cm ##cm = "colormap" para cambiar la paleta de color de la figura 3d
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)
<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x7f41c1cac710>
La lección del vídeo que te lleva a través de los detalles de los pasos 5 al 8 es la Video Lesson 6 en YouTube:
from IPython.display import YouTubeVideo
YouTubeVideo('tUg_dE3NXoY')
from IPython.core.display import HTML
def css_styling():
styles = open("../styles/custom.css", "r").read()
return HTML(styles)
css_styling()
(La celda de arriba establece el formato de este notebook)