建立河流-地下水系统中有机污染物 的对流、弥散及吸附作用的数学模型。考虑一维和二维，写出python代码的完整实现过程和建模思路并且给出代码解释

建立河流-地下水系统中有机污染物的对流、弥散及吸附作用的数学模型，可以使用对流-扩散-反应方程来描述。

对于一维情况，可以使用以下方程：

$$\frac{\partial C}{\partial t} = -v\frac{\partial C}{\partial x} + D\frac{\partial^2 C}{\partial x^2} - kC$$

其中，$C$表示有机污染物的浓度，$t$表示时间，$x$表示空间，$v$表示水流速度，$D$表示扩散系数，$k$表示吸附系数。

对于二维情况，可以使用以下方程：

$$\frac{\partial C}{\partial t} = -v_x\frac{\partial C}{\partial x} -v_y\frac{\partial C}{\partial y} + D\left(\frac{\partial^2 C}{\partial x^2} + \frac{\partial^2 C}{\partial y^2}\right) - kC$$

其中，$v_x$和$v_y$表示水流速度的$x$和$y$方向分量。

现在我们来考虑如何用 Python 实现这个模型。我们可以使用有限差分法来离散化这些方程，并使用显式欧拉方法来求解数值解。

对于一维情况，我们可以将空间离散化为$n$个网格点，时间离散化为$m$个时间步长。令$\Delta x$和$\Delta t$表示相邻网格点之间的距离，我们可以将方程离散化为：

$$\frac{C_i^{m+1}-C_i^m}{\Delta t} = -v\frac{C_i^m-C_{i-1}^m}{\Delta x} + D\frac{C_{i+1}^m-2C_i^m+C_{i-1}^m}{\Delta x^2} - kC_i^m$$

其中，$C_i^m$表示在时间步长$m$和空间点$i$处的浓度。我们可以通过重排这个方程来计算新的浓度$C_i^{m+1}$：

$$C_i^{m+1} = C_i^m + \Delta t\left(v\frac{C_i^m-C_{i-1}^m}{\Delta x} + D\frac{C_{i+1}^m-2C_i^m+C_{i-1}^m}{\Delta x^2} - kC_i^m\right)$$

这个式子可以用一个循环来计算所有网格点的浓度。我们还需要设置初始和边界条件来开始模拟。

对于二维情况，我们可以将空间离散化为$n_x\times n_y$个网格点，时间离散化为$m$个时间步长。令$\Delta x$和$\Delta t$表示相邻网格点之间的距离，我们可以将方程离散化为：

$$\frac{C_{i,j}^{m+1}-C_{i,j}^m}{\Delta t} = -v_x\frac{C_{i,j}^m-C_{i-1,j}^m}{\Delta x} -v_y\frac{C_{i,j}^m-C_{i,j-1}^m}{\Delta y} + D\left(\frac{C_{i+1,j}^m-2C_{i,j}^m+C_{i-1,j}^m}{\Delta x^2} + \frac{C_{i,j+1}^m-2C_{i,j}^m+C_{i,j-1}^m}{\Delta y^2}\right) - kC_{i,j}^m$$

我们可以通过重排这个方程来计算新的浓度$C_{i,j}^{m+1}$：

$$C_{i,j}^{m+1} = C_{i,j}^m + \Delta t\left(v_x\frac{C_{i-1,j}^m-C_{i,j}^m}{\Delta x} + v_y\frac{C_{i,j-1}^m-C_{i,j}^m}{\Delta y} + D\left(\frac{C_{i+1,j}^m-2C_{i,j}^m+C_{i-1,j}^m}{\Delta x^2} + \frac{C_{i,j+1}^m-2C_{i,j}^m+C_{i,j-1}^m}{\Delta y^2}\right) - kC_{i,j}^m\right)$$

这个式子可以用两个嵌套的循环来计算所有网格点的浓度。同样，我们还需要设置初始和边界条件来开始模拟。

以下是一个用 Python 实现一维情况的例子：

import numpy as np
import matplotlib.pyplot as plt

# 设置模拟参数
L = 1.0  # 河流长度
nx = 100  # 空间网格数
dx = L / (nx - 1)  # 空间步长
T = 100.0  # 模拟时间
nt = 1000  # 时间步长数
dt = T / nt  # 时间步长
v = 0.1  # 水流速度
D = 0.01  # 扩散系数
k = 0.1  # 吸附系数

# 设置初始和边界条件
C0 = np.zeros(nx)
C0[int(0.4*nx):int(0.6*nx)] = 1.0  # 污染物初始浓度
Cleft = 0.0  # 左边界条件
Cright = 0.0  # 右边界条件

# 初始化浓度矩阵
C = np.zeros((nt, nx))
C[0, :] = C0

# 循环计算浓度
for m in range(1, nt):
    for i in range(1, nx - 1):
        C[m, i] = C[m-1, i] + dt * (v*(C[m-1, i]-C[m-1, i-1])/dx + D*(C[m-1, i+1]-2*C[m-1, i]+C[m-1, i-1])/dx**2 - k*C[m-1, i])
    C[m, 0] = Cleft
    C[m, nx-1] = Cright

# 绘制浓度随时间变化的图像
plt.imshow(C.T, extent=[0, T, 0, L], cmap='plasma', aspect='auto')
plt.colorbar()
plt.xlabel('Time')
plt.ylabel('Position')
plt.show()

这个例子将河流长度设置为1.0，网格数设置为100，模拟时间设置为100.0，时间步长数设置为1000，水流速度设置为0.1，扩散系数设置为0.01，吸附系数设置为0.1。初始浓度设置为0，但在40%到60%的位置设置为1，左右边界条件都设置为0。在循环计算浓度时，使用了二阶中心差分来近似求解一阶和二阶导数。最后，将浓度随时间变化的图像绘制出来。

以下是一个用 Python 实现二维情况的例子：

import numpy as np
import matplotlib.pyplot as plt

# 设置模拟参数
Lx = 1.0  # 河流长度
Ly = 1.0  # 河流宽度
nx = 100  # x方向网格数
ny = 100  # y方向网格数
dx = Lx / (nx - 1)  # x方向步长
dy = Ly / (ny - 1)  # y方向步长
T = 100.0  # 模拟时间
nt = 1000  # 时间步长数
dt = T / nt  # 时间步长
vx = 0.1  # x方向水流速度
vy = 0.1  # y方向水流速度
D = 0.01  # 扩散系数
k = 0.1  # 吸附系数

# 设置初始和边界条件
C0 = np.zeros((nx, ny))
C0[int(0.4*nx):int(0.6*nx), int(0.4*ny):int(0.6*ny)] = 1.0  # 污染物初始浓度
Cleft = 0.0  # 左边界条件
Cright = 0.0  # 右边界条件
Ctop = 0.0  # 上边界条件
Cbottom = 0.0  # 下边界条件

# 初始化浓度矩阵
C = np.zeros((nt, nx, ny))
C[0, :, :] = C0

# 循环计算浓度
for m in range(1, nt):
    for i in range(1, nx - 1):
        for j in range(1, ny - 1):
            C[m, i, j] = C[m-1, i, j] + dt * (vx*(C[m-1, i-1, j]-C[m-1, i, j])/dx + vy*(C[m-1, i, j-1]-C[m-1, i, j])/dy + D*(C[m-1, i+1, j]-2*C[m-1, i, j]+C[m-1, i-1, j])/dx**2 + D*(C[m-1, i, j+1]-2*C[m-1, i, j]+C[m-1, i, j-1])/dy**2 - k*C[m-1, i, j])
    C[m, 0, :] = Cleft
    C[m, nx-1, :] = Cright
    C[m, :, 0] = Cbottom
    C[m, :, ny-1] = Ctop

# 绘制浓度随时间变化的图像
fig, ax = plt.subplots()
cax = ax.imshow(C[0, :, :].T, extent=[0, Lx, 0, Ly], cmap='plasma', aspect='auto')
ax.set_xlabel('X')
ax.set_ylabel('Y')
cbar = fig.colorbar(cax)
for m in range(1, nt):
    cax.set_data(C[m, :, :].T)
    plt.pause(0.01)
plt.show()

这个例子将河流长度和宽度都设置为1.0，x和y方向的网格数都设置为100，模拟时间设置为100.0，时间步长数设置为1000，x和y方向的水流速度都设置为0.1，扩散系数设置为0.01，吸附系数设置为0.1。初始浓度设置为0，但在40%到60%的位置和40%到60%的位置之间的区域设置为1，左右上下边界条件都设置为0。在循环计算浓度时，使用了二阶中心差分来近似求解一阶和二阶导数。最后，将浓度随时间变化的图像动态绘制出来。