python 使用有限差分求解一维热传导方程

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

$$\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}$$

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

使用有限差分方法求解一维热传导方程可以分为以下几个步骤：

1. 离散化：将时间和空间分别离散化为 $t_i = i \Delta t$ 和 $x_j = j \Delta x$，其中 $\Delta t$ 和 $\Delta x$ 分别表示时间和空间的步长。

2. 近似：将偏导数用差分来近似，例如：

$$\frac{\partial u}{\partial t} \approx \frac{u_{j,i+1} - u_{j,i}}{\Delta t}$$

$$\frac{\partial^2 u}{\partial x^2} \approx \frac{u_{j+1,i} - 2u_{j,i} + u_{j-1,i}}{\Delta x^2}$$

3. 替换：将近似后的偏导数代入原方程，得到差分方程：

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

4. 边界条件：确定边界条件，例如固定边界条件 $u(0,t) = u(L,t) = 0$。

5. 时间推进：使用差分方程逐步推进时间，例如使用显式欧拉法：

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

6. 初始条件：确定初始条件，例如 $u(x,0) = f(x)$。

下面是一个 Python 实现的示例代码：

import numpy as np
import matplotlib.pyplot as plt

# 设置参数
L = 1.0         # 空间长度
T = 0.1         # 时间长度
alpha = 1.0     # 热扩散系数
N = 100         # 空间网格数
M = 100         # 时间网格数
dx = L / N      # 空间步长
dt = T / M      # 时间步长

# 初始化温度分布
u = np.zeros((N+1, M+1))
u[0,:] = u[N,:] = 0.0
u[:,0] = np.sin(np.pi*np.arange(N+1)/(N+1))
# u[:,0] = np.exp(-100*(np.arange(N+1)/(N+1) - 0.5)**2)  # 高斯分布

# 使用有限差分方法求解热传导方程
for i in range(M):
    for j in range(1, N):
        u[j,i+1] = u[j,i] + alpha*dt/dx**2 * (u[j+1,i] - 2*u[j,i] + u[j-1,i])

# 绘制温度分布曲线
x = np.linspace(0, L, N+1)
t = np.linspace(0, T, M+1)
X, T = np.meshgrid(x, t)
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(1,1,1,projection='3d')
ax.plot_surface(X, T, u.T, cmap='coolwarm')
ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_zlabel('u')
plt.show()

其中，初始化温度分布的代码可以根据实际情况进行修改。绘制温度分布曲线的代码使用了 `matplotlib` 库中的 `plot_surface` 函数。