python实现对流-弥散方程

流-弥散方程是一类常见的偏微分方程，用于描述物质在流体中的扩散和传输过程。下面是一个简单的 Python 实现，可以用来模拟一维流-弥散方程的解。

首先，我们需要导入必要的库：

import numpy as np
import matplotlib.pyplot as plt

然后定义一些参数：

L = 1.0  # 区域长度
nx = 51  # 空间网格数
dx = L / (nx - 1)  # 空间步长
nt = 101  # 时间步数
nu = 0.3  # 扩散系数
sigma = 0.2  # CFL数
dt = sigma * dx**2 / nu  # 时间步长

接下来定义初始条件和边界条件：

u = np.ones(nx)
u[int(0.5 / dx):int(1 / dx + 1)] = 2  # 初始条件
u0 = np.ones(nx)
u0[int(0.5 / dx):int(1 / dx + 1)] = 2  # 保存初始条件

然后使用显式差分法进行求解：

for n in range(nt):
    un = u.copy()
    for i in range(1, nx-1):
        u[i] = un[i] + nu * dt / dx**2 * (un[i+1] - 2*un[i] + un[i-1])
    u[0] = 1
    u[-1] = 1

最后，我们可以画出结果：

plt.plot(np.linspace(0, L, nx), u0, label='Initial')
plt.plot(np.linspace(0, L, nx), u, label='Final')
plt.legend()
plt.show()

完整代码如下：

import numpy as np
import matplotlib.pyplot as plt

L = 1.0  # 区域长度
nx = 51  # 空间网格数
dx = L / (nx - 1)  # 空间步长
nt = 101  # 时间步数
nu = 0.3  # 扩散系数
sigma = 0.2  # CFL数
dt = sigma * dx**2 / nu  # 时间步长

u = np.ones(nx)
u[int(0.5 / dx):int(1 / dx + 1)] = 2  # 初始条件
u0 = np.ones(nx)
u0[int(0.5 / dx):int(1 / dx + 1)] = 2  # 保存初始条件

for n in range(nt):
    un = u.copy()
    for i in range(1, nx-1):
        u[i] = un[i] + nu * dt / dx**2 * (un[i+1] - 2*un[i] + un[i-1])
    u[0] = 1
    u[-1] = 1

plt.plot(np.linspace(0, L, nx), u0, label='Initial')
plt.plot(np.linspace(0, L, nx), u, label='Final')
plt.legend()
plt.show()

运行代码，我们可以得到如下结果：