Taller 7


In [1]:
%pylab inline


Populating the interactive namespace from numpy and matplotlib

Difusión

Ecuación a resolver $$ \frac{\partial u}{\partial t} = \nu\frac{\partial^2 u}{\partial x^2} $$

Discretizando

$$ \frac{\partial u}{\partial t} \approx \frac{u_{i}^{j+1} - u_{i}^{j}}{\Delta t} $$$$ \frac{\partial^2 u}{\partial x^2} \approx \frac{u_{i+1}^j - 2 u_{i}^{j} + u_{i-1}^{j}}{(\Delta x)^2} $$

En el esquema de diferencias finitas $$ \frac{u_{i}^{j+1} - u_{i}^{j}}{\Delta t} = \nu \frac{u_{i+1}^j - 2 u_{i}^{j} + u_{i-1}^{j}}{(\Delta x)^2} $$

$$ u_{i}^{j+1} = \nu\alpha u_{i+1}^{j} + (1 - 2\nu\alpha)u_{i}^{j} + \nu\alpha u_{i-1}^{j} $$

donde $\alpha = \Delta t/ (\Delta x)^2$.


In [25]:
n_x = 100
n_t = 100

nu = 0.3
sigma = 0.2 #parámetro que asegura que \alpha\nu < 0.5

x = linspace(0, 1.0, n_x)
dx = x[1]-x[0]


dt = sigma*dx**2/nu #dt se define utilizando sigma
alpha = dt/dx**2
print dt

u = zeros(n_x)

u[where((x<0.66) & (x>0.33))] = 1.0


6.80202700405e-05

In [26]:
plt.plot(x,u)
plt.title('Condicion inicial')
plt.xlabel('x')
plt.ylabel('u')


Out[26]:
<matplotlib.text.Text at 0x3c6a410>

In [27]:
for n in range(n_t):  # loop over time
    u_past = u.copy() 
    for i in range(1,n_x-1): #loop over space
        u[i] = nu * alpha * u_past[i+1]  + (1.0 - 2.0*nu*alpha)*u_past[i] + nu*alpha*u_past[i-1]

In [28]:
plt.plot(x,u)
plt.title('Condicion despes de $\Delta t=100$')
plt.xlabel('x')
plt.ylabel('u')


Out[28]:
<matplotlib.text.Text at 0x44ece10>

In [ ]: