对于这个方程∂t ∂u（x，t）=∂x² ∂²u（x，t）（0<x<1，0<t）它满足u（0，t）=u（1，t）=0和u（x，0）=4x（1-x）以这个方程为例，取x方向的步长h=0.2，取t方向的步长τ=1/150，计算j=1，2两层上的值，用古典显式差分方法写一个python程序

好的，根据您的要求，我编写了一个使用古典显式差分方法求解偏微分方程的Python程序。请注意，由于此方程的边界条件是周期性的，因此我们需要在左侧和右侧分别使用相同的边界条件，即u[0][j]=u[-1][j]=0。

# 导入必要的库
import numpy as np
import matplotlib.pyplot as plt

# 设置参数
L = 1  # 区间长度
T = 1  # 时间长度
n = 6  # x方向网格数
m = 150  # t方向网格数
h = L / (n - 1)  # x方向步长
tau = T / m  # t方向步长
r = tau / (h ** 2)  # 稳定性参数

# 创建初始值矩阵U
U = np.zeros((n, m))
U[:, 0] = 4 * np.arange(0, 1 + h, h) * (1 - np.arange(0, 1 + h, h))

# 对于每一层j，使用古典显式差分方法求解
for j in range(1, m):
    for i in range(1, n - 1):
        U[i][j] = U[i][j - 1] + r * (U[i - 1][j - 1] - 2 * U[i][j - 1] + U[i + 1][j - 1])
    U[0][j] = U[-1][j] = 0

# 输出结果
print("U[1][1] = ", U[1][1])
print("U[4][1] = ", U[4][1])
print("U[1][2] = ", U[1][2])
print("U[4][2] = ", U[4][2])

# 画图
x = np.linspace(0, L, n)
t = np.linspace(0, T, m)
X, T = np.meshgrid(x, t)
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X, T, U, cmap='coolwarm', linewidth=0, antialiased=False)
ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_zlabel('u')
plt.show()

程序运行后，输出了在第1层和第2层上x=0.2和x=0.8处的U值，以及使用matplotlib库绘制的3D图形。您可以根据需要调整n和m的值，以便获得更准确的结果和更高的精度。