Convección en una dimensión, con álgebra lineal

... igual que el otro día



In [1]:

    
%matplotlib inline



In [2]:

    
import scipy as np
from matplotlib import pyplot as plt

Parámetros



In [3]:

    
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.1            # tiempo total
nt = 100            # pasos temporales
dt = T / nt

c = 1              # velocidad de la onda

¡Número de Courant !



In [4]:

    
Co = c * dt / dx
Co









    Out[4]:





0.04

Condiciones iniciales



In [5]:

    
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 [6]:

    
plt.plot( x , u0 )









    Out[6]:





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

Un paso en el tiempo

Recordemos que queremos implementar $u_i^{n+1} = u_i^n - \mathrm{Co}/2 (u_{i+1}^n-u_{i-1}^n)$



In [7]:

    
u = u0.copy()



In [8]:

    
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 / 2.0) * (un[i+1] - un[i-1])

¡Empiezan los cambios!

... pero eso se parece a una operación de álgebra lineal



In [9]:

    
D=np.eye( nx )



In [10]:

    
plt.matshow( D )









    Out[10]:





<matplotlib.image.AxesImage at 0x7efecb7d5f50>



In [11]:

    
D = D + np.diag( -(Co/2.0)*np.ones( nx - 1 ) , k=1 ) + np.diag( (Co/2.0)*np.ones( nx - 1 ) , k=-1 )



In [14]:

    
D









    Out[14]:





array([[ 1.  , -0.02,  0.  , ...,  0.  ,  0.  ,  0.  ],
       [ 0.02,  1.  , -0.02, ...,  0.  ,  0.  ,  0.  ],
       [ 0.  ,  0.02,  1.  , ...,  0.  ,  0.  ,  0.  ],
       ..., 
       [ 0.  ,  0.  ,  0.  , ...,  1.  , -0.02,  0.  ],
       [ 0.  ,  0.  ,  0.  , ...,  0.02,  1.  , -0.02],
       [ 0.  ,  0.  ,  0.  , ...,  0.  ,  0.02,  1.  ]])



In [17]:

    
plt.matshow( D )









    Out[17]:





<matplotlib.image.AxesImage at 0x7efecb501d10>



In [22]:

    
plt.matshow( D[ 1:10, 1:10 ] )









    Out[22]:





<matplotlib.image.AxesImage at 0x7efecb10e450>



In [23]:

    
u = u0.copy()

cada paso temporal puede verse como la aplicación de la matriz $D$ al vector $u$, o sea, $D \cdot u$



In [24]:

    
u = np.dot( D , u)



In [25]:

    
plt.plot(x,u)









    Out[25]:





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

Tiempo completo



In [26]:

    
u = u0.copy()



In [27]:

    
for n in range(nt):
    u = np.dot( D , u)



In [28]:

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









    Out[28]:





[<matplotlib.lines.Line2D at 0x7efecafbfed0>,
 <matplotlib.lines.Line2D at 0x7efecafcd110>]

Formulaciones alternativas

Probemos ahora el algoritmo upwind: $u_i^{n+1} = u_i^n - \mathrm{Co} (u_{i}^n-u_{i-1}^n)$



In [33]:

    
u = u0.copy()



In [35]:

    
D=  ( 1 - Co ) * np.eye( nx )
D = D + np.diag( Co * np.ones( nx - 1 ) , k=-1 )
D









    Out[35]:





array([[ 0.96,  0.  ,  0.  , ...,  0.  ,  0.  ,  0.  ],
       [ 0.04,  0.96,  0.  , ...,  0.  ,  0.  ,  0.  ],
       [ 0.  ,  0.04,  0.96, ...,  0.  ,  0.  ,  0.  ],
       ..., 
       [ 0.  ,  0.  ,  0.  , ...,  0.96,  0.  ,  0.  ],
       [ 0.  ,  0.  ,  0.  , ...,  0.04,  0.96,  0.  ],
       [ 0.  ,  0.  ,  0.  , ...,  0.  ,  0.04,  0.96]])



In [36]:

    
plt.matshow( D[ 1:10, 1:10 ] )









    Out[36]:





<matplotlib.image.AxesImage at 0x7efecae298d0>



In [37]:

    
for n in range(nt):
    u = np.dot( D , u)



In [38]:

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









    Out[38]:





[<matplotlib.lines.Line2D at 0x7efecadd3dd0>,
 <matplotlib.lines.Line2D at 0x7efecadd3fd0>]

Métodos implícitos

... Esto no parece servir de mucho, pero ¿qué sucede si consideramos un método implícito ?

$u_i^{n+1} = u_i^n - \mathrm{Co} (u_{i}^{n+1}-u_{i-1}^{n+1})$

$u_i^{n+1} (1 + \mathrm{Co} ) - \mathrm{Co} \, u_{i-1}^{n+1} = u_i^n $

Este es un problema de álgebra lineal, $ D u' = u$ (otra matriz $D$, pero mantenemos el símbolo)



In [40]:

    
D=  (1+Co) * np.eye( nx )
D = D - np.diag( Co * np.ones( nx - 1 ) , k=-1 )
D









    Out[40]:





array([[ 1.04,  0.  ,  0.  , ...,  0.  ,  0.  ,  0.  ],
       [-0.04,  1.04,  0.  , ...,  0.  ,  0.  ,  0.  ],
       [ 0.  , -0.04,  1.04, ...,  0.  ,  0.  ,  0.  ],
       ..., 
       [ 0.  ,  0.  ,  0.  , ...,  1.04,  0.  ,  0.  ],
       [ 0.  ,  0.  ,  0.  , ..., -0.04,  1.04,  0.  ],
       [ 0.  ,  0.  ,  0.  , ...,  0.  , -0.04,  1.04]])



In [41]:

    
plt.matshow( D[ 1:10, 1:10 ] )









    Out[41]:





<matplotlib.image.AxesImage at 0x7efecacfdf50>



In [42]:

    
u = u0.copy()



In [43]:

    
from scipy import linalg



In [44]:

    
u = np.linalg.solve( D , u )



In [45]:

    
plt.plot(x , u )









    Out[45]:





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



In [46]:

    
u = u0.copy()



In [47]:

    
for n in range(nt):
    u = np.linalg.solve( D , u )



In [48]:

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









    Out[48]:





[<matplotlib.lines.Line2D at 0x7efec882f250>,
 <matplotlib.lines.Line2D at 0x7efec882f450>]

Reflexiones

De nuevo, ¿qué hemos hecho en los bordes? (condiciones de contorno)
¿Por qué puede ser interesante un método explícito?