Advection equation
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}$$
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]:
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)
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)
In [110]:
print(B_C)
In [111]:
C_record = []
C_record.append(C)
print(C)
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)
In [113]:
pyplot.ylim((0., 10.))
pyplot.xlabel('x')
pyplot.ylabel('concentration')
pyplot.plot(x_grid, C)
pyplot.show()
In [31]:
print(x_grid)
In [ ]: