求助：修正基于有限差分法的扩散方程数值模拟绘图结果

我有一套求解1D扩散方程（$\partial^2 C/\partial x^2 - \partial C/\partial t = 0$）的代码，离散格式为：
$$C[n+1,j] = C[n,j] + \frac{dt}{dx^2} \times (C[n,j+1] - 2C[n,j] + C[n,j-1])$$

目标是生成随时间演化的浓度分布曲线，但当前绘图结果与预期不符，请求修正绘图逻辑。以下是原代码：

import numpy as np
import matplotlib.pyplot as plt

dt = 0.001 # grid size for time (s)
dx = 0.05 # grid size for space (m)
x_max = 1 # in m
t_max = 1 # total time in s
C0 = 1 # concentration
 
# function to calculate concentration profiles based on a 
# finite difference approximation to the 1D diffusion 
# equation:
def diffusion(dt,dx,t_max,x_max,C0):
    # diffusion number:
    s = dt/dx**2
    x = np.arange(0,x_max+dx,dx) 
    t = np.arange(0,t_max+dt,dt)
    r = len(t)
    a = len(x)
    C = np.zeros([r,a]) # initial condition
    C[:,0] = C0 # boundary condition on left side
    C[:,-1] = 0 # boundary condition on right side
    for n in range(0,r-1): # time
        for j in range(1,a-1): # space
            C[n+1,j] = C[n,j] + s*(C[n,j-1] -
                2*C[n,j] + C[n,j+1]) 
    return x,C,r,a
# note that this can be written without the for-loop 
# in space, but it is easier to read it this way
 
x,C,r,a = diffusion(dt,dx,t_max,x_max,C0)
 
# plotting:
plt.figure()
plt.xlim([0,1])
plt.ylim([0,1])
plot_times = np.arange(0,1,0.02)
for t in plot_times:
    plt.plot(x,C[int(t/dt),:],'Gray',label='numerical')
plt.xlabel('Membrane position x',fontsize=12)
plt.ylabel('Concentration',fontsize=12)

问题分析与修正方案 #

原代码的数值计算部分是稳定的（扩散数$s=0.4<0.5$，满足显式格式稳定性条件），问题出在绘图逻辑上：

1. 时间采样过于密集（每隔0.02秒绘制一条曲线），50条灰色曲线重叠后无法区分演化趋势
2. 重复添加相同图例，导致图例无效
3. 可视化样式缺乏区分度，无法直观观察时间演化

以下是修正后的完整代码：

import numpy as np
import matplotlib.pyplot as plt

dt = 0.001  # 时间步长 (s)
dx = 0.05   # 空间步长 (m)
x_max = 1   # 空间范围 (m)
t_max = 1   # 总时间 (s)
C0 = 1      # 左边界浓度

def diffusion(dt, dx, t_max, x_max, C0):
    s = dt / dx**2  # 扩散数
    x = np.arange(0, x_max + dx, dx)
    t = np.arange(0, t_max + dt, dt)
    r, a = len(t), len(x)
    C = np.zeros([r, a])
    C[:, 0] = C0       # 左边界：x=0处浓度保持C0
    C[:, -1] = 0       # 右边界：x=x_max处浓度保持0
    
    for n in range(r-1):
        for j in range(1, a-1):
            C[n+1, j] = C[n, j] + s * (C[n, j-1] - 2*C[n, j] + C[n, j+1])
    return x, C, t

x, C, t = diffusion(dt, dx, t_max, x_max, C0)

# 优化绘图逻辑
plt.figure(figsize=(8, 6))
plt.xlim([0, x_max])
plt.ylim([0, C0])

# 每隔0.1秒采样一次，减少曲线数量
sample_interval = 0.1
plot_indices = np.arange(0, len(t), int(sample_interval / dt))
# 用渐变颜色区分不同时间的曲线
colors = plt.cm.viridis(np.linspace(0, 1, len(plot_indices)))

for idx, color in zip(plot_indices, colors):
    current_time = t[idx]
    plt.plot(x, C[idx, :], color=color, label=f't={current_time:.1f}s')

plt.xlabel('膜位置 x (m)', fontsize=12)
plt.ylabel('浓度 C', fontsize=12)
plt.title('1D扩散方程数值解：浓度随时间演化', fontsize=14)
# 将图例放在图外，避免遮挡曲线
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

修正说明 #

• 增大时间采样间隔至0.1秒，减少曲线数量，清晰展示演化趋势
• 使用颜色渐变区分不同时间点的曲线，直观体现时间推进
• 为每条曲线标注对应时间，图例移至图外避免遮挡
• 添加网格和标题，优化整体可视化效果

内容的提问来源于stack exchange，提问作者geo.freitas