二维热传导方程初边值问题数值求解，给出Crank-Nicholson离散格式和python代码

二维热传导方程描述了物体内部温度随时间和空间的变化过程，其基本形式是时间偏微分方程。初边值问题是指在给定初始条件和边界条件下，求解这个方程。Crank-Nicolson方法是一种半隐式有限差分法，它结合了向前、向后的显式和隐式步骤，用于求解稳定性好且收敛速度快。

Crank-Nicholson公式的具体离散格式可以这样表达：
假设我们要计算网格点 (i, j) 的下一个时间步 (n+1) 的温度 \( T_{i,j}^{(n+1)} \)，则有：
\[ \frac{T_{i,j}^{(n+1)} - T_{i,j}^{(n)}}{\Delta t} = \alpha \left(\frac{T_{i+1,j}^{(n+\frac{1}{2})} + T_{i-1,j}^{(n+\frac{1}{2})} - 2T_{i,j}^{(n+\frac{1}{2})}}{\Delta x^2} + \frac{T_{i,j+1}^{(n+\frac{1}{2})} + T_{i,j-1}^{(n+\frac{1}{2})} - 2T_{i,j}^{(n+\frac{1}{2})}}{\Delta y^2}\right) \]

这里，\( \Delta t \) 是时间步长，\( \Delta x \) 和 \( \Delta y \) 是空间步长，\( \alpha \) 是热扩散系数。\( T_{i+1,j}^{(n+\frac{1}{2})} \) 等是邻近节点的时间平滑值，通常通过平均当前时间步和下一时间步的值得到。

Python代码示例（简化版）如下：

import numpy as np

def crank_nicolson(dt, dx, dy, alpha, T_n, T_left, T_right, T_top, T_bottom):
    T_half = 0.5 * (T_n[1:-1, 1:-1] + T_n[:-2, 1:-1] + T_n[2:, 1:-1] + T_n[1:-1, :-2] + T_n[1:-1, 2:])
    
    # 更新中间网格点
    for i in range(1, T_n.shape[0]-1):
        for j in range(1, T_n.shape[1]-1):
            T_next[i, j] = T_n[i, j] + dt * alpha * (dx**2 * (T_half[i+1, j] - 2*T_half[i, j] + T_half[i-1, j]) + dy**2 * (T_half[i, j+1] - 2*T_half[i, j] + T_half[i, j-1]))

    # 应用边界条件
    T_next[0, :] = T_left
    T_next[-1, :] = T_right
    T_next[:, 0] = T_bottom
    T_next[:, -1] = T_top
    
    return T_next

# 假设已有的变量初始化和网格设置...