3.利用 Crank-Nicolson 格式求解如下的抛物型方程(给出算法格式，附上代码，计算结果和图像)：\[ \begin{cases} \frac{\partial u}{\partial t} =u_{xx}+u_{yy},&0 \leq x,y \leq 1\\u\mid _{ \partial{G} } \ =0.\\u(x,y,0)=sin(\pi(x))sin(\pi{y}).\end{cases} \]

Crank-Nicolson 格式是一种经典的求解偏微分方程的数值方法，其可以通过求解一个线性方程组来计算下一个时间步的解。对于二维抛物型方程，Crank-Nicolson 格式的离散化形式为：

$$
\frac{U_{i,j}^{n+1}-U_{i,j}^n}{\Delta t}=\frac{1}{2}\left[\frac{U_{i+1,j}^n-2U_{i,j}^n+U_{i-1,j}^n}{\Delta x^2}+\frac{U_{i,j+1}^n-2U_{i,j}^n+U_{i,j-1}^n}{\Delta y^2}\right.\\
\left.+\frac{U_{i+1,j+1}^{n+1}-2U_{i,j}^{n+1}+U_{i-1,j-1}^{n+1}}{4\Delta x \Delta y}+\frac{U_{i+1,j-1}^{n+1}-2U_{i,j}^{n+1}+U_{i-1,j+1}^{n+1}}{4\Delta x \Delta y}\right]
$$

其中，$U_{i,j}^n$ 表示在时间 $t_n$ 和空间 $(i\Delta x,j\Delta y)$ 处的解，$\Delta x$ 和 $\Delta y$ 分别表示空间的步长，$\Delta t$ 表示时间的步长。

根据上式，我们可以得到关于 $U_{i,j}^{n+1}$ 的一个三对角线性方程组，可以使用追赶法（Thomas 算法）进行求解。求解该方程组的算法步骤如下：

1. 确定时间和空间步长 $\Delta t$、$\Delta x$、$\Delta y$，以及网格的大小 $N$。

2. 初始化解 $U_{i,j}^0=sin(\pi i\Delta x)sin(\pi j\Delta y)$。

3. 通过迭代求解三对角线性方程组，得到 $U_{i,j}^{n+1}$。

4. 重复步骤 3，直到达到所需的时间步数 $T/\Delta t$。

下面是 Python 代码实现：

import numpy as np
import matplotlib.pyplot as plt

# 空间和时间的步长
dx = 0.01
dy = 0.01
dt = 0.0001

# 网格大小
N = int(1/dx)
M = int(1/dy)

# 时间步数
T = 0.1
timesteps = int(T/dt)

# 初始化解
u = np.zeros((N+1,M+1))
for i in range(N+1):
    for j in range(M+1):
        u[i,j] = np.sin(np.pi*i*dx)*np.sin(np.pi*j*dy)

# 迭代求解
for n in range(timesteps):
    # 构造三对角线性方程组
    A = np.zeros((N+1,M+1))
    B = np.zeros((N+1,M+1))
    C = np.zeros((N+1,M+1))
    D = np.zeros((N+1,M+1))
    for i in range(1,N):
        for j in range(1,M):
            A[i,j] = -0.5*dt/dx**2
            B[i,j] = 1+dt/dx**2+dt/dy**2
            C[i,j] = -0.5*dt/dx**2
            D[i,j] = u[i,j] + 0.5*dt/dx**2*(u[i+1,j]-2*u[i,j]+u[i-1,j])\
                     + 0.5*dt/dy**2*(u[i,j+1]-2*u[i,j]+u[i,j-1])
    # 处理边界条件
    for i in range(N+1):
        A[i,0] = 0
        B[i,0] = 1
        C[i,0] = 0
        D[i,0] = 0
        A[i,M] = 0
        B[i,M] = 1
        C[i,M] = 0
        D[i,M] = 0
    for j in range(M+1):
        A[0,j] = 0
        B[0,j] = 1
        C[0,j] = 0
        D[0,j] = 0
        A[N,j] = 0
        B[N,j] = 1
        C[N,j] = 0
        D[N,j] = 0
    # 追赶法求解三对角线性方程组
    for j in range(1,M):
        alpha = np.zeros(N+1)
        beta = np.zeros(N+1)
        alpha[1] = A[1,j]/B[1,j]
        beta[1] = D[1,j]/B[1,j]
        for i in range(2,N):
            alpha[i] = A[i,j]/(B[i,j]-C[i,j]*alpha[i-1])
            beta[i] = (D[i,j]+C[i,j]*beta[i-1])/(B[i,j]-C[i,j]*alpha[i-1])
        u[N,j] = (D[N,j]+C[N,j]*beta[N-1])/(B[N,j]-C[N,j]*alpha[N-1])
        for i in range(N-1,0,-1):
            u[i,j] = alpha[i]*u[i+1,j]+beta[i]
    # 更新解
    u[:,0] = u[:,1]
    u[:,-1] = u[:,-2]
    u[0,:] = u[1,:]
    u[-1,:] = u[-2,:]

# 绘制解的图像
x = np.linspace(0,1,N+1)
y = np.linspace(0,1,M+1)
X,Y = np.meshgrid(x,y)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, u, cmap='winter')
plt.show()

最终得到的解的图像如下所示：