python求解二维热传导方程

二维热传导方程可以表示为：

$$\frac{\partial u}{\partial t} = \alpha\left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right)$$

其中，$u(x,y,t)$ 表示温度分布，$\alpha$ 表示热扩散系数。

我们可以使用有限差分法来求解这个方程。假设我们在一个矩形区域内进行求解，将这个矩形区域划分为 $N_x \times N_y$ 个网格，其中第 $i$ 行第 $j$ 列的网格的坐标为 $(x_i, y_j)$，时间步长为 $\Delta t$。我们用 $u_{i,j}^n$ 表示在 $(x_i, y_j)$ 处的温度在第 $n$ 个时间步的值。则有以下的差分方程：

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

其中，$\Delta x$ 和 $\Delta y$ 分别表示网格在 $x$ 和 $y$ 方向上的间距。

下面是 Python 代码实现：

import numpy as np
import matplotlib.pyplot as plt

# 矩形区域的边界
Lx = 1.0
Ly = 1.0

# 网格数
Nx = 41
Ny = 41

# 网格间距
dx = Lx / (Nx - 1)
dy = Ly / (Ny - 1)

# 时间步长和总时间
dt = 0.01
T = 1.0

# 热扩散系数
alpha = 0.01

# 初始化温度分布
u = np.zeros((Nx, Ny))
u[int(Nx / 4):int(3 * Nx / 4), int(Ny / 4):int(3 * Ny / 4)] = 1.0

# 迭代求解
nsteps = int(T / dt)
for n in range(nsteps):
    un = u.copy()
    u[1:-1, 1:-1] = un[1:-1, 1:-1] + alpha * dt / dx**2 * (un[2:, 1:-1] - 2*un[1:-1, 1:-1] + un[:-2, 1:-1]) + alpha * dt / dy**2 * (un[1:-1, 2:] - 2*un[1:-1, 1:-1] + un[1:-1, :-2])

# 绘制结果
x = np.linspace(0, Lx, Nx)
y = np.linspace(0, Ly, Ny)
X, Y = np.meshgrid(x, y)
plt.contourf(X, Y, u)
plt.show()

运行代码后，会得到一个矩形区域内的温度分布图像。