二维热传导方程非齐次形式初边值问题程序

下面是一个求解二维热传导方程非齐次形式初边值问题的程序，使用有限差分方法：

import numpy as np
import matplotlib.pyplot as plt

# 定义初始条件
nx = 50
ny = 50
nt = 1000
dt = 0.01
dx = 0.1
dy = 0.1
x = np.linspace(0, dx*nx, nx)
y = np.linspace(0, dy*ny, ny)
X, Y = np.meshgrid(x, y)
T = np.zeros((ny, nx, nt))
T[:, :, 0] = np.sin(X)*np.sin(Y)

# 定义边界条件
T[:, 0, :] = 0
T[:, -1, :] = 0
T[0, :, :] = 0
T[-1, :, :] = 0

# 定义热传导系数和非齐次项
k = 1
q = np.zeros((ny, nx))
q[ny//3:2*ny//3, nx//3:2*nx//3] = 1

# 迭代求解
for n in range(1, nt):
    T[1:-1, 1:-1, n] = T[1:-1, 1:-1, n-1] + \
        k*dt/dx**2*(T[1:-1, 2:, n-1]-2*T[1:-1, 1:-1, n-1]+T[1:-1, 0:-2, n-1]) + \
        k*dt/dy**2*(T[2:, 1:-1, n-1]-2*T[1:-1, 1:-1, n-1]+T[0:-2, 1:-1, n-1]) + \
        dt*q[1:-1, 1:-1]

# 绘制结果
plt.imshow(T[:, :, -1], cmap='hot', origin='lower')
plt.colorbar()
plt.show()

其中，T是温度分布的三维数组，第一维表示y坐标，第二维表示x坐标，第三维表示时间。q是非齐次项的二维数组。程序使用了简单的欧拉显式差分方法求解热传导方程，通过迭代求解得到最终的温度分布，并使用matplotlib绘制结果。