python 求解一维热传导方程

一维热传导方程可以用偏微分方程来描述：

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

其中，$u(x,t)$ 是温度分布，$\alpha$ 是热扩散系数。为了求解这个方程，我们需要设定初始条件和边界条件。假设初始温度分布为 $u(x,0) = f(x)$，边界条件为 $u(0,t) = u(L,t) = 0$。

我们可以使用有限差分法来求解这个方程。假设我们把时间分成 $N$ 个时间步长，把空间分成 $M$ 个网格点，那么我们可以把偏微分方程转化成一个差分方程：

$$\frac{u_i^{n+1} - u_i^n}{\Delta t} = \alpha \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\Delta x^2}$$

其中，$u_i^n$ 表示在时间步 $n$ 和空间点 $i$ 处的温度。

使用这个差分方程，我们可以递推求解出整个系统在时间 $t$ 处的温度分布。下面是一个用 Python 实现的一维热传导方程的求解代码：

import numpy as np
import matplotlib.pyplot as plt

# 设置参数和初始条件
L = 1
alpha = 1
T = 1
M = 100
N = 1000
dx = L / M
dt = T / N
x = np.linspace(0, L, M+1)
u = np.zeros((M+1, N+1))
u[:,0] = np.sin(np.pi*x)

# 使用差分方程递推求解
for n in range(N):
    for i in range(1, M):
        u[i,n+1] = u[i,n] + alpha * dt / dx**2 * (u[i+1,n] - 2*u[i,n] + u[i-1,n])
    u[0,n+1] = 0
    u[M,n+1] = 0

# 绘制结果
plt.plot(x, u[:,0], label='t=0')
plt.plot(x, u[:,-1], label='t=1')
plt.legend()
plt.show()

运行上面的代码，可以得到一个图像，显示初始温度分布和 $t=1$ 时的温度分布。可以看到，随着时间的推移，温度分布逐渐趋于平稳状态。