我对使用 numpy 和 sympy 库还是很陌生。抱歉，如果我在那里打印了很多照片，我只是想检查代码是否正常工作。我正在尝试解决这个 2D 热方程问题，并且有点难以理解我如何添加初始条件（30 度的温度）并添加均匀的狄利克雷边界条件：温度到板的每一侧（即顶部，底部和两侧）。例如，我尝试将板顶部的温度添加为 1000 度，但出现错误“IndexError: index 18 is out of bounds for axis 0 with size 18”。我已经设法将函数添加到代码中，它似乎也可以输出，只是在添加条件时有点迷失......我试过在网上寻找类似的问题，但不能
这是原始二维热方程的链接：https ://drive.google.com/file/d/1q3KyfaHD7hZTbdIYc9_e5otusELqYJaD/view?usp=sharing
这是显式有限差分的最终方程的链接：https ://drive.google.com/file/d/1unJHXO3b3xig8RhBen5k0eOg40YUfjyv/view?usp=sharing
i;j 是位置（节点编号）和 k= 时间（时间步长编号）
我也在试图找出温度达到稳定状态的时间
我已经添加了到目前为止我所做的整个代码
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import sympy as sp
from sympy.solvers import solve
from sympy import Function, Symbol
#Length of workpiece
Lx= 0.05 #Length of work piece in the x direction, measured in meters (m)
Nx= 18 #Number of mesh points along the x-direction
Ly= 0.1 #Length of work piece in the y direction, measured in meters (m)
Ny= 9 #Number of mesh points along the y-direction
dy=Lx/Nx #distance between nodes in x
dx=Ly/Ny #distance between nodes in x
#Defining material properties
conduc= 16.26 #Conductivity of the material in (W/m (degrees)C)
rho= 8020 #Density of the material in (kg/m^3)
c=502 #Specific heat of the material (J/kg(degrees)C)
#Calculating the alpha (thermal diffusivity) for the plate, in (m^2/s)
a=((conduc)/(rho*c))
print (a)
t=215;
nt= 215; #number of steps wrt time
dt= t/nt; #timestep
#Define the mesh in space
x= np.linspace(1,Lx, Nx)[np.newaxis] #vector to create mesh
y= np.linspace(1,Ly,Ny)[np.newaxis]
x,y=np.meshgrid(x,y) #mesh grid function return
print (x,y)
T = np.ones(([Nx,Ny]))
#Adding the other conditions
T[Nx,:]= 1000 #Assining the temp for the top
sym= sp.symbols('i, j, k')
function = sp.Function('T(i,j,k)+(dt*a(((T(i+1,j,k)-2*T(i,j,k)+T(i-1,j,k))/(dx**2))+(a*dt(((T(i,j+1,k)-2*T(i,k,k)+T(i,j-1,k))/dy**2)')
print(function)
更新：我已将此添加到代码中，这些似乎有效。但是，我注意到我
T = np.ones(([Nx,Ny]))的尺寸是 18 x 9，但是水平边的盘子更长。这是否意味着我只需要翻转矩阵？
#Temperatures for the plate
Ttop=1000 # temp for the top side of the plate
Tbottom = 500 # temnp for the bottom side of the plate
Tright = 500 #temp for the right side of the plate
Tleft = 1000 #temp for the left side of the plate
#Adding the boundaries on the system
T[:, 0] = Ttop #temp for the top
print (T[:,0])
T[:, -1] = Tbottom #temp for bottom
print (T[:, -1])
T[0, :] = Tright #temp for the right
print (T[0, :])
T[-1, :] = Tleft #temp for the left
print (T[-1, :])