我想通过有限差分法确定时间导数和拉普拉斯算子,以解决一维情况下的热方程。我的目标是得到源项,这就是为什么我需要得到方程的所有左项。
当我将源项绘制为时间的函数时,结果令人惊讶……与预期相去甚远。这是 np.reshape 的问题吗?感谢帮助 ;)
#profil1D_theta.shape[0] : time
#profil1D_theta.shape[1] : space
#profil1D_theta_reshape.shape
#Out[986]: (328, 411)
T_tempo = 1./115. #time period
#Time derivative
dt_theta_list = []
for i in range(profil1D_theta.shape[0]-1):
for j in range(profil1D_theta.shape[1]-1):
dt_theta = (profil1D_theta[i+1,j] - profil1D_theta[i,j])/T_tempo
dt_theta_list.append(dt_theta)
dt_theta = np.asarray(dt_theta_list)
dt_theta_reshape = dt_theta.reshape(328,411)
#Laplacian
T_spatial = 0.000021661 #spatial period
dx2test_theta_list =[]
for i in range(0,profil1D_theta.shape[0]-1):
for j in range(0,profil1D_theta.shape[1]-1):
dxtest_theta = (profil1D_theta[i,j+1] - 2*profil1D_theta[i,j] + profil1D_theta[i,j-1])/T_spatial**2
dx2test_list.append(dxtest_theta)
dx2test_theta = np.asarray(dx2test_theta_list)
dx2_theta_reshape = dx2test_theta.reshape(328,411)
rho = 7850.
C = 500.
k = 42.
D = k/(rho*C)
e = 3.*0.001
b = 6.*0.001
Boltz = 5.679373*10**(-8)
Emissivite_echant = 0.97 #
T0 = 273.15 + 26.
h = 5.
tau_th = rho*C*e/(2*h+8*Boltz*Emissivite_echant*T0**3)
#Sources
profil1D_theta_reshape = profil1D_theta[0:328,0:411]
w1_chaleur_list=[]
for i in range(0,profil1D_theta_reshape.shape[1]):
for j in range(0,profil1D_theta_reshape.shape[0]):
w_chaleur1D = rho*C*(dt_theta_reshape[j,i]+1./(tau_th)*profil1D_theta_reshape[j,i] - D*dx_theta_reshape[j,i])
w1_chaleur_list.append(w_chaleur1D)
d,e = dt_theta_reshape.shape
w1_chaleur = np.asarray(w1_chaleur_list)
w1_chaleur_reshape = np.reshape(w1_chaleur,(e,d))
clim_min = 0
clim_max = 1*10**6
plt.imshow(w1_chaleur_reshape);plt.clim([clim_min,clim_max]);plt.colorbar();
plt.ylabel(r'Axe longitudinal')
plt.xlabel(r'Temps (images)')
plt.show()