如何实现污染物分段释放并求解含衰减的一维扩散方程

我需要模拟一个在房间中心释放10分钟后关闭的污染源，采用含扩散与衰减的一维扩散方程是合理的：
$$\frac{dC}{dt} = D \cdot \frac{d^2C}{dx2} - k \cdot C$$
当初始条件设为C(50)=1（即初始时刻仅中心位置有浓度）时，求解这个方程很容易，但我不知道该如何设置分段的源项条件——也就是让污染源只在指定时间段内释放，之后停止。

初始版本代码（单次初始浓度） #

对应仅初始时刻中心有浓度的模拟：

import numpy as np

def point_source_pde(C, D, k, dx, dt):
    """求解含扩散与衰减的连续点源PDE方程
    
    参数:
    C (ndarray): t时刻各位置x的浓度分布
    D (float): 扩散系数
    k (float): 衰减速率
    dx (float): 空间离散步长
    dt (float): 时间离散步长
    
    返回:
    ndarray: t+dt时刻各位置x的浓度分布
    """
    # 网格点数量
    N = C.shape[0]
    
    # 初始化更新后的浓度数组
    C_new = np.zeros(N)
    
    # 遍历所有网格点
    for i in range(N):
        # 扩散项的有限差分计算
        C_diffusion = D * (C[(i+1)%N] - 2*C[i] + C[(i-1)%N]) / dx**2
        
        # 衰减项计算
        C_decay = -k * C[i]
        
        # 更新每个网格点的浓度
        C_new[i] = C[i] + dt * (C_diffusion + C_decay)
    
    return C_new


# 设置扩散系数与衰减速率
D = 0.1
k = 0.01

# 设置空间与时间离散步长
dx = 0.1
dt = 0.001

# 初始化浓度分布：仅中心位置（第50个网格）浓度为1
C = np.zeros(100)
C[50] = 1

# 迭代求解1000个时间步
for t in range(1000):
    C = point_source_pde(C, D, k, dx, dt)

# 绘制浓度分布
import matplotlib.pyplot as plt
plt.plot(C)
plt.xlabel('距离 [x]')
plt.ylabel('浓度 C')
plt.show()

运行上述代码得到的浓度分布：

尝试方案1的问题分析 #

你给出的尝试方案思路正确，但存在几个关键问题：

• 源项施加位置错误：代码把释放量加到了所有网格点，实际污染源仅在中心位置释放
• 时间概念混淆：用时间步的序号判断释放时段，未和实际释放时长（10分钟）对应
• 代码格式错误：函数定义的缩进不符合Python规范，会导致运行报错

修正后的实现方案 #

以下是调整后的代码，解决了上述问题，精准实现定时释放的源项条件：

import numpy as np
import matplotlib.pyplot as plt

def point_source_pde(C, D, k, dx, dt, current_time, release_start, release_end, source_pos):
    """求解含扩散与衰减的定时点源PDE方程，仅在指定时间段内释放污染源
    
    参数:
    C (ndarray): 当前时刻各位置的浓度分布
    D (float): 扩散系数
    k (float): 衰减速率
    dx (float): 空间离散步长
    dt (float): 时间离散步长
    current_time (float): 当前时刻的实际时间（单位：分钟）
    release_start (float): 污染源开始释放的实际时间（分钟）
    release_end (float): 污染源停止释放的实际时间（分钟）
    source_pos (int): 污染源所在的网格位置索引
    
    返回:
    ndarray: 下一时刻的浓度分布
    """
    N = C.shape[0]
    C_new = np.copy(C)
    
    # 计算所有网格点的扩散与衰减项
    for i in range(N):
        C_diffusion = D * (C[(i+1)%N] - 2*C[i] + C[(i-1)%N]) / dx**2
        C_decay = -k * C[i]
        C_new[i] += dt * (C_diffusion + C_decay)
    
    # 仅在释放时间段内，给中心位置添加源项
    if release_start <= current_time < release_end:
        # 源强度可根据需求调整，此处为单位时间释放量
        source_strength = 1.0
        C_new[source_pos] += dt * source_strength
    
    return C_new

# -------------------------- 参数设置 --------------------------
# 物理参数
D = 0.1                  # 扩散系数
k = 0.01                 # 衰减速率（1/分钟）
release_duration = 10    # 污染源释放时长（10分钟）
release_start = 0        # 从0时刻开始释放

# 离散参数
dx = 0.1                 # 空间步长
dt = 0.01                # 时间步长（0.01分钟=0.6秒）
total_time = 60          # 模拟总时长（60分钟）
num_time_steps = int(total_time / dt)

# 空间网格设置：共100个点，中心位置为第50个索引
num_grid = 100
source_pos = 50
C = np.zeros(num_grid)   # 初始浓度全为0

# -------------------------- 模拟求解 --------------------------
for step in range(num_time_steps):
    current_time = step * dt
    C = point_source_pde(C, D, k, dx, dt, current_time, release_start, release_start + release_duration, source_pos)

# -------------------------- 结果绘制 --------------------------
plt.plot(C)
plt.xlabel('距离 [x]')
plt.ylabel('浓度 C')
plt.title(f"污染源释放{release_duration}分钟后的浓度分布")
plt.show()

关键修正说明 #

1. 精准的源项施加：仅在中心网格点添加释放量，符合实际污染源位置
2. 实际时间对应：用实际时间（分钟）判断释放时段，直接匹配需求中的“10分钟释放”
3. 逻辑清晰的分段控制：通过current_time判断是否处于释放期，灵活控制源项的开启/关闭
4. 代码格式规范：修复缩进问题，保证代码可正常运行

内容的提问来源于stack exchange，提问作者HCAI