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?