给出一个用Crank Nicolson求二维热传导方程的例子

好的，下面给出一个使用 Crank-Nicolson 格式求解二维热传导方程的例子。热传导方程的一般形式为：

$$\frac{\partial u}{\partial t} = \nabla \cdot (\kappa \nabla u) + f$$

其中，$u$ 表示温度，$\kappa$ 表示热传导系数，$f$ 表示源项。在这里，我们假设 $\kappa$ 为常数，$f=0$，并考虑一个矩形区域 $[0,L_x] \times [0,L_y]$ 上的热传导问题。那么，我们可以将上述方程离散化为以下差分方程：

$$\frac{u_{i,j}^{n+1} - u_{i,j}^n}{\Delta t} = \frac{\kappa}{2} \left( \frac{u_{i+1,j}^{n+1} - 2u_{i,j}^{n+1} + u_{i-1,j}^{n+1}}{\Delta x^2} + \frac{u_{i,j+1}^{n+1} - 2u_{i,j}^{n+1} + u_{i,j-1}^{n+1}}{\Delta y^2} + \frac{u_{i+1,j+1}^n - u_{i+1,j-1}^n - u_{i-1,j+1}^n + u_{i-1,j-1}^n}{4 \Delta x \Delta y} \right)$$

其中，$u_{i,j}^n$ 表示在时刻 $n$，位置 $(i,j)$ 处的温度，$\Delta t$、$\Delta x$ 和 $\Delta y$ 分别表示时间、$x$ 和 $y$ 的步长。

我们可以将上式写成以下形式：

$$- \frac{\kappa \Delta t}{2 \Delta x^2} u_{i-1,j}^{n+1} + \left(1 + \frac{\kappa \Delta t}{\Delta x^2} + \frac{\kappa \Delta t}{\Delta y^2} \right) u_{i,j}^{n+1} - \frac{\kappa \Delta t}{2 \Delta x^2} u_{i+1,j}^{n+1} = \frac{\kappa \Delta t}{2 \Delta x^2} u_{i-1,j}^{n} + \left(1 - \frac{\kappa \Delta t}{\Delta x^2} - \frac{\kappa \Delta t}{\Delta y^2} \right) u_{i,j}^{n} + \frac{\kappa \Delta t}{2 \Delta x^2} u_{i+1,j}^{n} + \frac{\kappa \Delta t}{2 \Delta y^2} (u_{i,j+1}^n + u_{i,j-1}^n) + \frac{\kappa \Delta t}{4 \Delta x \Delta y} (u_{i+1,j+1}^n - u_{i+1,j-1}^n - u_{i-1,j+1}^n + u_{i-1,j-1}^n)$$

我们可以使用迭代法来求解上式，具体来说，我们可以使用以下代码：

import numpy as np

# 参数设置
T = 1.0
Lx = 1.0
Ly = 1.0
Nt = 100
Nx = 50
Ny = 50
kappa = 1.0

# 网格设置
dt = T / Nt
dx = Lx / Nx
dy = Ly / Ny
x = np.linspace(0, Lx, Nx+1)
y = np.linspace(0, Ly, Ny+1)
t = np.linspace(0, T, Nt+1)
u = np.zeros((Nt+1, Nx+1, Ny+1))

# 初始条件和边界条件
u[0, :, :] = np.sin(np.pi * x[:, np.newaxis]) * np.sin(np.pi * y[np.newaxis, :])
u[:, 0, :] = 0
u[:, Nx, :] = 0
u[:, :, 0] = 0
u[:, :, Ny] = 0

# 迭代求解
rx = kappa * dt / (2 * dx ** 2)
ry = kappa * dt / (2 * dy ** 2)
rxy = kappa * dt / (4 * dx * dy)

for n in range(Nt):
    A = np.zeros((Nx-1, Ny-1, Nx-1, Ny-1))
    b = np.zeros((Nx-1, Ny-1))
    for i in range(1, Nx):
        for j in range(1, Ny):
            if i == 1:
                A[i-1, j-1, 0, :] = 0
                A[i-1, j-1, 1, 0] = 1 + rx + ry
                A[i-1, j-1, 1, 1:] = -rx / 2
                A[i-1, j-1, 2:, :] = 0
                b[i-1, j-1] = (1 - rx - ry) * u[n, i, j] + rx / 2 * u[n, i+1, j] + ry / 2 * (u[n, i, j+1] + u[n, i, j-1]) + rxy * (u[n, i+1, j+1] - u[n, i+1, j-1] - u[n, i-1, j+1] + u[n, i-1, j-1])
            elif i == Nx-1:
                A[i-1, j-1, -1, :] = 0
                A[i-1, j-1, -2, -1] = 1 + rx + ry
                A[i-1, j-1, -2, :-1] = -rx / 2
                A[i-1, j-1, :-2, :] = 0
                b[i-1, j-1] = (1 - rx - ry) * u[n, i, j] + rx / 2 * u[n, i-1, j] + ry / 2 * (u[n, i, j+1] + u[n, i, j-1]) + rxy * (u[n, i+1, j+1] - u[n, i+1, j-1] - u[n, i-1, j+1] + u[n, i-1, j-1])
            elif j == 1:
                A[i-1, j-1, :, 0] = 0
                A[i-1, j-1, :, 1] = -ry / 2
                A[i-1, j-1, :, 2] = 1 + rx + ry
                A[i-1, j-1, :, 3:] = -ry / 2
                b[i-1, j-1] = (1 - rx - ry) * u[n, i, j] + rx / 2 * (u[n, i+1, j] + u[n, i-1, j]) + ry / 2 * u[n, i, j+1] + rxy * (u[n, i+1, j+1] - u[n, i+1, j-1] - u[n, i-1, j+1] + u[n, i-1, j-1])
            elif j == Ny-1:
                A[i-1, j-1, :, -1] = 0
                A[i-1, j-1, :, -2] = -ry / 2
                A[i-1, j-1, :, -3] = 1 + rx + ry
                A[i-1, j-1, :, :-3] = -ry / 2
                b[i-1, j-1] = (1 - rx - ry) * u[n, i, j] + rx / 2 * (u[n, i+1, j] + u[n, i-1, j]) + ry / 2 * u[n, i, j-1] + rxy * (u[n, i+1, j+1] - u[n, i+1, j-1] - u[n, i-1, j+1] + u[n, i-1, j-1])
            else:
                A[i-1, j-1, :, :] = [[-ry/2, -rx/2, 0, -ry/2, rxy, 0, 0, 0],
                                      [-rx/2, 1+rx+ry, -rx/2, 0, -rxy, rxy, 0, 0],
                                      [0, -rx/2, -ry/2, 0, 0, -rxy, ry/2, 0],
                                      [-ry/2, 0, 0, -ry/2, 0, 0, -rx/2, rxy],
                                      [rxy, -rxy, 0, 0, 1+2*rxy, 0, -rxy, 0],
                                      [0, rxy, -rxy, 0, 0, -ry/2, 0, ry/2],
                                      [0, 0, ry/2, -rx/2, -rxy, 0, -rx/2, -ry/2],
                                      [0, 0, 0, rxy, 0, ry/2, -ry/2, -ry/2]]
                b[i-1, j-1] = (1 - rx - ry) * u[n, i, j] + rx / 2 * (u[n, i+1, j] + u[n, i-1, j]) + ry / 2 * (u[n, i, j+1] + u[n, i, j-1]) + rxy * (u[n, i+1, j+1] - u[n, i+1, j-1] - u[n, i-1, j+1] + u[n, i-1, j-1])
    A = A.reshape(((Nx-1)*(Ny-1), (Nx-1)*(Ny-1)))
    b = b.reshape(((Nx-1)*(Ny-1),))
    u[n+1, 1:Nx, 1:Ny] = np.linalg.solve(A, b).reshape((Nx-1, Ny-1))

# 绘制图像
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
X, Y = np.meshgrid(x, y)
ax.plot_surface(X, Y, u[-1, :, :])
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('u')
plt.show()

在上述代码中，我们使用了 $50 \times 50$ 的网格来进行求解，时间步长为 $0.01$，空间步长为 $0.02$。我们使用 $\Delta t / \Delta x^2 = \Delta t / \Delta y^2 = 0.5$ 的 Crank-Nicolson 格式来进行求解。初始条件为 $u(x, y, 0) = \sin(\pi x) \sin(\pi y)$，边界条件为 $u(x, 0, t) = u(x, L_y, t) = u(0, y, t) = u(L_x, y, t) = 0$。

最终