代码存在报错 D[:] = u[i-1, 1:-1] + alpha*(u[i-1, :-2] - 2*u[i-1, 1:-1] + u[i-1, 2:]) ValueError: could not broadcast input array from shape (99,) into shape (100,)

非常抱歉，这是我的疏忽，代码中的 `D` 应该是包含 $N-1$ 个元素的数组，而不是包含 $N$ 个元素的数组，因此在计算 `D` 时出现了维度不匹配的错误。以下是修改后的代码：

import numpy as np
import matplotlib.pyplot as plt

# 定义常数
L = 1  # 计算区间长度
N = 101  # 离散点数
dx = L / (N - 1)  # 离散点间距
dt = 0.001  # 时间步长
T = 0.1  # 计算时长
nt = int(T/dt) + 1  # 时间步数
nu = 0.1  # 粘度系数
Re = 5000  # 雷诺数
u0 = np.zeros(N)  # 初始速度

# 设置边界条件
u0[0] = 0  # 下边界条件
u0[-1] = 1  # 上边界条件

# 定义函数
def thomas(A, B, C, D):
    n = len(B)
    c_ = np.zeros(n-1)
    d_ = np.zeros(n)
    u = np.zeros(n)
    c_[0] = C[0] / B[0]
    d_[0] = D[0] / B[0]
    for i in range(1, n-1):
        c_[i] = C[i] / (B[i] - A[i]*c_[i-1])
    for i in range(1, n):
        d_[i] = (D[i] - A[i-1]*d_[i-1]) / (B[i] - A[i-1]*c_[i-1])
    u[-1] = d_[-1]
    for i in range(n-2, -1, -1):
        u[i] = d_[i] - c_[i]*u[i+1]
    return u

def CrankNicolson(u0, nu, dx, dt, nt):
    u = np.zeros((nt, N))
    u[0] = u0
    alpha = nu*dt/(2*dx**2)
    A = np.zeros(N-2)
    B = np.zeros(N-1)
    C = np.zeros(N-2)
    D = np.zeros(N-1)
    for i in range(1, nt):
        A[:] = -alpha
        B[:] = 1 + 2*alpha
        C[:] = -alpha
        D[:-1] = u[i-1, 1:-1] + alpha*(u[i-1, :-2] - 2*u[i-1, 1:-1] + u[i-1, 2:])
        D[-1] = u[i-1, -1]  # 上边界条件
        u[i, 1:-1] = thomas(A, B, C, D)
        u[i, 0] = 0  # 下边界条件
        u[i, -1] = 1  # 上边界条件
    return u

# 计算数值解
u = CrankNicolson(u0, nu, dx, dt, nt)

# 绘制曲线
plt.plot(np.linspace(0, L, N), u[-1])
plt.xlabel('位置')
plt.ylabel('速度')
plt.show()

在计算 `D` 时，我们首先计算前 $N-1$ 个元素，然后将最后一个元素设置为上边界条件。这样就可以避免维度不匹配的错误了。