Difusión en una dimensión

In [1]:

    

%matplotlib inline


    

In [2]:

    

import scipy as np
from matplotlib import pyplot as plt


    

In [18]:

    

L  = 1.0           # longitud del sistema 1D
nx = 42            # nodos espaciales
dx = L / (nx-2)    # sí, quitamos dos nodos ...
x = np.linspace( 0 , L , num=nx )

T= 0.02            # tiempo total
nt = 400            # pasos temporales
dt = T / nt

c = 1              # velocidad de la onda


    

In [22]:

    

Co = c * dt / dx**2
Co


    

    
Out[22]:





0.07999999999999999

In [23]:

    

u0 = 1 * np.ones(nx)                    # todo uno
x1 = L/4 ; n1 = int(x1 / dx)
x2 = L/2 ; n2 = int(x2 / dx)
u0[ n1 : n2 ] = 2


    

In [24]:

    

plt.plot( x , u0 )


    

    
Out[24]:





[<matplotlib.lines.Line2D at 0x7f5b70f4c9d0>]






    

In [25]:

    

u = u0.copy()


    

In [26]:

    

un = u.copy()         # distribución actual

for i in range( 1 , nx - 1 ):    # Ahora queda claro por qué hemos quitado los extremos !!
   u[i] = un[i] + Co * (un[i+1] + un[i-1] - 2 * un[i] )


    

In [27]:

    

plt.plot(x,u)


    

    
Out[27]:





[<matplotlib.lines.Line2D at 0x7f5b70e83610>]






    

In [28]:

    

u = u0.copy()


    

In [29]:

    

for n in range(nt):
    un = u.copy()
    for i in range( 1 , nx - 1 ): 
       u[i] = un[i] + Co * (un[i+1] + un[i-1] - 2 * un[i] )


    

In [30]:

    

plt.plot(x , u , x , u0 , 'r')


    

    
Out[30]:





[<matplotlib.lines.Line2D at 0x7f5b70da8f90>,
 <matplotlib.lines.Line2D at 0x7f5b70db61d0>]