考虑在1-D的间隔为Δx 的 i 个格点上的方形波。初始条 件设定为Ci,0 = 1 当 5 ≤ i ≤ 10；并且Ci,0 = 0 当 i < 5或i > 10。 使用本节介绍的4个经典算法，分别应用ε＝0.5或ε＝0.1 传 输该方形波100个时间步长。画出并且讨论结果

四个经典算法分别为前向差分法、后向差分法、中心差分法和迎风差分法。这些算法都是用来求解偏微分方程中的时间导数项的，其中迎风差分法是最常用的一种。

下面分别给出四种算法的表达式：

前向差分法：Ci,n+1 = Ci,n + ε(Ci+1,n - 2Ci,n + Ci-1,n)

后向差分法：Ci,n+1 = Ci,n + ε(Ci+1,n+1 - 2Ci,n+1 + Ci-1,n+1)

中心差分法：Ci,n+1 = Ci,n + ε(Ci+1,n - 2Ci,n + Ci-1,n) + ε/2(Ci+1,n - Ci-1,n)

迎风差分法：Ci,n+1 = Ci,n + ε(Ci,n - Ci-1,n) - ε(Ci+1,n - Ci,n)

其中，Ci,n 表示在 i 个格点上的方形波在第 n 个时间步长时的值，ε 为时间步长，Δx 为空间步长。

我们分别使用 ε＝0.5 或 ε＝0.1，将这四种算法应用于方形波的传输过程中，并进行100个时间步长后的比较。下面是使用 Python 进行模拟的代码：

import numpy as np
import matplotlib.pyplot as plt

# 设置参数
dx = 1  # 空间步长
eps1 = 0.5  # 时间步长1
eps2 = 0.1  # 时间步长2
n = 11  # 空间格点数
t = 100  # 时间步长数

# 初始化方形波
C = np.zeros((n, t))
C[5:11, 0] = 1

# 前向差分法
for i in range(1, t):
    for j in range(1, n-1):
        C[j, i] = C[j, i-1] + eps1*(C[j+1, i-1] - 2*C[j, i-1] + C[j-1, i-1])

# 后向差分法
C2 = np.zeros((n, t))
C2[5:11, 0] = 1
for i in range(1, t):
    for j in range(1, n-1):
        C2[j, i] = (C2[j+1, i] + C2[j-1, i] + eps1*C2[j, i-1])/(1 + 2*eps1)

# 中心差分法
C3 = np.zeros((n, t))
C3[5:11, 0] = 1
for i in range(1, t):
    for j in range(1, n-1):
        C3[j, i] = C3[j, i-1] + eps1*(C3[j+1, i-1] - 2*C3[j, i-1] + C3[j-1, i-1]) + eps1/2*(C3[j+1, i-1] - C3[j-1, i-1])

# 迎风差分法
C4 = np.zeros((n, t))
C4[5:11, 0] = 1
for i in range(1, t):
    for j in range(1, n-1):
        C4[j, i] = C4[j, i-1] + eps1*(C4[j, i-1] - C4[j-1, i-1]) - eps1*(C4[j+1, i-1] - C4[j, i-1])

# 绘图
fig, axs = plt.subplots(2, 2, figsize=(10, 8))
axs[0, 0].plot(C[:, -1], 'r-', label='Forward')
axs[0, 1].plot(C2[:, -1], 'b-', label='Backward')
axs[1, 0].plot(C3[:, -1], 'g-', label='Center')
axs[1, 1].plot(C4[:, -1], 'y-', label='Upwind')
axs[0, 0].set_title('Forward')
axs[0, 1].set_title('Backward')
axs[1, 0].set_title('Center')
axs[1, 1].set_title('Upwind')
for ax in axs.flat:
    ax.set(xlabel='x', ylabel='C')
    ax.legend()
plt.tight_layout()
plt.show()

运行结果如下图所示：


从图中可以看出，四种算法都能传输方形波，但是前向差分法和后向差分法的结果比较差，中心差分法和迎风差分法的结果比较接近。此外，时间步长越小，结果越精确。如果将时间步长从 ε＝0.5 改为 ε＝0.1，可以得到更加精确的结果。