Advection equation

$$\frac {\partial C}{\partial t} = - w \frac {\partial C}{\partial x}$$$$\frac {\partial C}{\partial x} \Bigg |_{x = 0, L} = 0$$

Discretise the equation

$$\frac {C^{t+1}_{x} - C^{t}_{x}}{\Delta t} = -\frac {w}{2} \Bigg( \frac {C^{t+1}_{x+1} - C^{t+1}_{x-1}}{2 \Delta x} + \frac {C^{t}_{x+1} - C^{t}_{x-1}}{2 \Delta x} \Bigg)$$$$C^{t+1}_{x} = C^{t}_{x} - \frac {w \Delta t}{4\Delta x} \Bigg(C^{t+1}_{x+1} - C^{t+1}_{x-1} + C^{t}_{x+1} - C^{t}_{x-1} \Bigg)$$

if $\frac {w \Delta t}{4\Delta x} = \sigma$

$$C^{t+1}_{x} = C^{t}_{x} - \sigma C^{t+1}_{x+1} + \sigma C^{t+1}_{x-1} - \sigma C^{t}_{x+1} + \sigma C^{t}_{x-1}$$

and rearranging,

$$-\sigma C^{t+1}_{x-1} + C^{t+1}_{x} + \sigma C^{t+1}_{x+1} = \sigma C^{t}_{x-1} + C^{t}_{x} - \sigma C^{t}_{x+1}$$

gives the discretised equation for $x > 1, x < L$.

Discretise the boundary conditions

$$\frac {C^t_{0} - C^t_{-1}}{\Delta x} = 0$$$$C^t_{0} = C^t_{-1}$$$$\frac {C^t_{L-1} - C^t_{L}}{\Delta x} = 0$$$$C^t_{L-1} = C^t_{L}$$

So, from the above equation, at $x = 0$ we have

$$-\sigma C^{t+1}_{-1} + C^{t+1}_{0} + \sigma C^{t+1}_{1} = \sigma C^{t}_{-1} + C^{t}_{0} - \sigma C^{t}_{1}$$

and substituting in $C^t_{0} = C^t_{-1}$ and $C^{t+1}_{0} = C^{t+1}_{-1}$

$$-\sigma C^{t+1}_{0} + C^{t+1}_{0} + \sigma C^{t+1}_{1} = \sigma C^{t}_{0} + C^{t}_{0} - \sigma C^{t}_{1}$$$$(1 - \sigma)C^{t+1}_{0} + \sigma C^{t+1}_{1} = (1 + \sigma) C^{t}_{0} - \sigma C^{t}_{1}$$

And at $x = L-1$ we have

$$-\sigma C^{t+1}_{L-2} + C^{t+1}_{L-1} + \sigma C^{t+1}_{L} = \sigma C^{t}_{L-2} + C^{t}_{L-1} - \sigma C^{t}_{L}$$

and substituting in $C^t_{L} = C^t_{L-1}$ and $C^{t+1}_{L} = C^{t+1}_{L-1}$

$$-\sigma C^{t+1}_{L-2} + C^{t+1}_{L-1} + \sigma C^{t+1}_{L-1} = \sigma C^{t}_{L-2} + C^{t}_{L-1} - \sigma C^{t}_{L-1}$$$$-\sigma C^{t+1}_{L-2} + (1 + \sigma)C^{t+1}_{L-1} = \sigma C^{t}_{L-2} + (1 - \sigma) C^{t}_{L-1}$$
$$ \begin{bmatrix} 1-\sigma & \sigma & 0 & 0 & \cdots & 0 & 0 & 0 & 0\\\\ -\sigma & 1 & \sigma & 0 & \cdots & 0 & 0 & 0 & 0 \\\\ 0 & -\sigma & 1 & \sigma & \cdots & 0 & 0 & 0 & 0 \\\\ 0 & 0 & \ddots & \ddots & \ddots & \ddots & \ddots & 0 & 0 \\\\ 0 & 0 & 0 & 0 & \cdots & 0 & -\sigma & 1 & \sigma \\\\ 0 & 0 & 0 & 0 & \cdots & 0 & 0 & -\sigma & 1+\sigma \end{bmatrix} \begin{bmatrix} C_0^{t+1} \\\\ C_1^{t+1} \\\\ C_2^{t+1} \\\\ \vdots \\\\ C_{L-2}^{t+1} \\\\ C_{L-1}^{t+1} \end{bmatrix} = \begin{bmatrix} 1+\sigma & -\sigma & 0 & 0 & \cdots & 0 & 0 & 0 & 0\\\\ \sigma & 1 & -\sigma & 0 & \cdots & 0 & 0 & 0 & 0 \\\\ 0 & \sigma & 1 & -\sigma & \cdots & 0 & 0 & 0 & 0 \\\\ 0 & 0 & \ddots & \ddots & \ddots & \ddots & 0 & 0 & 0 \\\\ 0 & 0 & 0 & 0 & 0 & 0 & \sigma & 1 & -\sigma \\\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \sigma & 1-\sigma \end{bmatrix} \begin{bmatrix} C_0^{t} \\\\ C_1^{t} \\\\ C_2^{t} \\\\ \vdots \\\\ C_{L-2}^{t} \\\\ C_{L-1}^{t} \end{bmatrix} $$

In [98]:
import numpy
from matplotlib import pyplot

In [99]:
numpy.set_printoptions(precision=3)

In [100]:
L = 1.
J = 100
dx = float(L)/float(J-1)
x_grid = numpy.array([j*dx for j in range(J)])

In [101]:
T = 200
N = 1000
dt = float(T)/float(N-1)
t_grid = numpy.array([n*dt for n in range(N)])

In [102]:
w = 0.1
sigma = float(w*dt)/float((4.*dx))

In [105]:
w*dt/4*dx
dt


Out[105]:
0.2002002002002002

In [106]:
#C =  numpy.zeros(J)
#C[10:20] = 5

gaussian = lambda z, height, position, hwhm: height * numpy.exp(-numpy.log(2) * ((z - position)/hwhm)**2)
C = gaussian(x_grid, 5, 0.5, 0.1)

In [107]:
print(C)


[  1.490e-07   2.980e-07   5.876e-07   1.142e-06   2.190e-06   4.138e-06
   7.710e-06   1.416e-05   2.566e-05   4.582e-05   8.068e-05   1.401e-04
   2.398e-04   4.046e-04   6.733e-04   1.105e-03   1.787e-03   2.850e-03
   4.481e-03   6.947e-03   1.062e-02   1.600e-02   2.378e-02   3.484e-02
   5.032e-02   7.167e-02   1.006e-01   1.393e-01   1.902e-01   2.560e-01
   3.397e-01   4.444e-01   5.733e-01   7.291e-01   9.142e-01   1.130e+00
   1.378e+00   1.656e+00   1.962e+00   2.293e+00   2.641e+00   3.000e+00
   3.359e+00   3.709e+00   4.037e+00   4.333e+00   4.585e+00   4.784e+00
   4.921e+00   4.991e+00   4.991e+00   4.921e+00   4.784e+00   4.585e+00
   4.333e+00   4.037e+00   3.709e+00   3.359e+00   3.000e+00   2.641e+00
   2.293e+00   1.962e+00   1.656e+00   1.378e+00   1.130e+00   9.142e-01
   7.291e-01   5.733e-01   4.444e-01   3.397e-01   2.560e-01   1.902e-01
   1.393e-01   1.006e-01   7.167e-02   5.032e-02   3.484e-02   2.378e-02
   1.600e-02   1.062e-02   6.947e-03   4.481e-03   2.850e-03   1.787e-03
   1.105e-03   6.733e-04   4.046e-04   2.398e-04   1.401e-04   8.068e-05
   4.582e-05   2.566e-05   1.416e-05   7.710e-06   4.138e-06   2.190e-06
   1.142e-06   5.876e-07   2.980e-07   1.490e-07]

In [108]:
A_C = numpy.diagflat([-sigma for i in range(J-1)], -1) +\
      numpy.diagflat([1.-sigma]+[1 for i in range(J-2)]+[1.+sigma]) +\
      numpy.diagflat([sigma for i in range(J-1)], 1)
        
B_C = numpy.diagflat([sigma for i in range(J-1)], -1) +\
      numpy.diagflat([1.+sigma]+[1. for i in range(J-2)]+[1.-sigma]) +\
      numpy.diagflat([-sigma for i in range(J-1)], 1)

In [109]:
print(A_C)


[[ 0.505  0.495  0.    ...,  0.     0.     0.   ]
 [-0.495  1.     0.495 ...,  0.     0.     0.   ]
 [ 0.    -0.495  1.    ...,  0.     0.     0.   ]
 ..., 
 [ 0.     0.     0.    ...,  1.     0.495  0.   ]
 [ 0.     0.     0.    ..., -0.495  1.     0.495]
 [ 0.     0.     0.    ...,  0.    -0.495  1.495]]

In [110]:
print(B_C)


[[ 1.495 -0.495  0.    ...,  0.     0.     0.   ]
 [ 0.495  1.    -0.495 ...,  0.     0.     0.   ]
 [ 0.     0.495  1.    ...,  0.     0.     0.   ]
 ..., 
 [ 0.     0.     0.    ...,  1.    -0.495  0.   ]
 [ 0.     0.     0.    ...,  0.495  1.    -0.495]
 [ 0.     0.     0.    ...,  0.     0.495  0.505]]

In [111]:
C_record = []

C_record.append(C)



print(C)


[  1.490e-07   2.980e-07   5.876e-07   1.142e-06   2.190e-06   4.138e-06
   7.710e-06   1.416e-05   2.566e-05   4.582e-05   8.068e-05   1.401e-04
   2.398e-04   4.046e-04   6.733e-04   1.105e-03   1.787e-03   2.850e-03
   4.481e-03   6.947e-03   1.062e-02   1.600e-02   2.378e-02   3.484e-02
   5.032e-02   7.167e-02   1.006e-01   1.393e-01   1.902e-01   2.560e-01
   3.397e-01   4.444e-01   5.733e-01   7.291e-01   9.142e-01   1.130e+00
   1.378e+00   1.656e+00   1.962e+00   2.293e+00   2.641e+00   3.000e+00
   3.359e+00   3.709e+00   4.037e+00   4.333e+00   4.585e+00   4.784e+00
   4.921e+00   4.991e+00   4.991e+00   4.921e+00   4.784e+00   4.585e+00
   4.333e+00   4.037e+00   3.709e+00   3.359e+00   3.000e+00   2.641e+00
   2.293e+00   1.962e+00   1.656e+00   1.378e+00   1.130e+00   9.142e-01
   7.291e-01   5.733e-01   4.444e-01   3.397e-01   2.560e-01   1.902e-01
   1.393e-01   1.006e-01   7.167e-02   5.032e-02   3.484e-02   2.378e-02
   1.600e-02   1.062e-02   6.947e-03   4.481e-03   2.850e-03   1.787e-03
   1.105e-03   6.733e-04   4.046e-04   2.398e-04   1.401e-04   8.068e-05
   4.582e-05   2.566e-05   1.416e-05   7.710e-06   4.138e-06   2.190e-06
   1.142e-06   5.876e-07   2.980e-07   1.490e-07]

In [112]:
for ti in range(1,2):
    C_new = numpy.linalg.solve(A_C, B_C.dot(C))
    
    C = C_new
    C
    C_record.append(C)

print(C)


[  8.107e-08   6.919e-08   1.043e-07   2.003e-07   4.035e-07   8.094e-07
   1.600e-06   3.113e-06   5.957e-06   1.121e-05   2.077e-05   3.787e-05
   6.796e-05   1.200e-04   2.088e-04   3.576e-04   6.030e-04   1.002e-03
   1.639e-03   2.640e-03   4.192e-03   6.555e-03   1.010e-02   1.533e-02
   2.293e-02   3.379e-02   4.907e-02   7.021e-02   9.901e-02   1.376e-01
   1.885e-01   2.544e-01   3.384e-01   4.437e-01   5.734e-01   7.304e-01
   9.170e-01   1.135e+00   1.384e+00   1.665e+00   1.973e+00   2.306e+00
   2.656e+00   3.016e+00   3.376e+00   3.726e+00   4.053e+00   4.348e+00
   4.598e+00   4.793e+00   4.927e+00   4.993e+00   4.989e+00   4.915e+00
   4.774e+00   4.572e+00   4.318e+00   4.021e+00   3.692e+00   3.342e+00
   2.984e+00   2.626e+00   2.280e+00   1.952e+00   1.648e+00   1.372e+00
   1.126e+00   9.116e-01   7.278e-01   5.731e-01   4.451e-01   3.409e-01
   2.576e-01   1.919e-01   1.410e-01   1.022e-01   7.310e-02   5.156e-02
   3.587e-02   2.462e-02   1.667e-02   1.113e-02   7.338e-03   4.771e-03
   3.060e-03   1.937e-03   1.210e-03   7.456e-04   4.535e-04   2.722e-04
   1.612e-04   9.428e-05   5.441e-05   3.100e-05   1.744e-05   9.679e-06
   5.319e-06   2.852e-06   1.593e-06   6.768e-07]

In [113]:
pyplot.ylim((0., 10.))
pyplot.xlabel('x')
pyplot.ylabel('concentration')
pyplot.plot(x_grid, C)
pyplot.show()

In [31]:
print(x_grid)


[ 0.     0.01   0.02   0.03   0.04   0.051  0.061  0.071  0.081  0.091
  0.101  0.111  0.121  0.131  0.141  0.152  0.162  0.172  0.182  0.192
  0.202  0.212  0.222  0.232  0.242  0.253  0.263  0.273  0.283  0.293
  0.303  0.313  0.323  0.333  0.343  0.354  0.364  0.374  0.384  0.394
  0.404  0.414  0.424  0.434  0.444  0.455  0.465  0.475  0.485  0.495
  0.505  0.515  0.525  0.535  0.545  0.556  0.566  0.576  0.586  0.596
  0.606  0.616  0.626  0.636  0.646  0.657  0.667  0.677  0.687  0.697
  0.707  0.717  0.727  0.737  0.747  0.758  0.768  0.778  0.788  0.798
  0.808  0.818  0.828  0.838  0.848  0.859  0.869  0.879  0.889  0.899
  0.909  0.919  0.929  0.939  0.949  0.96   0.97   0.98   0.99   1.   ]

In [ ]: