我见过其他人问类似的事情,但我无论如何都无法弄清楚问题所在。 我正在尝试将 Matlab 代码转换为 Python,但在 FOR 循环中的以下行之后遇到了问题:
dx = abs(np.diff(g_coord(num
(((下面是该循环的代码。任何帮助将不胜感激。我真的试图自己修复它,但没有成功。对不起,如果这是一个愚蠢的错误。MATLAB行保留为Python注释,以防万一。
import numpy as np
from scipy.sparse import lil_matrix
# physical parameters
seconds_per_yr = 60*60*24*365; # number of seconds in one year
lx = 10000 ; #length of spatial domain (m)
Cp = 1e3 ; # rock heat capacity (J/kg/K)
rho = 2700 ; # rock density (kg/mˆ3)
K = 3.3 ; # bulk thermal conductivity (W/m/K)
kappa = K/(Cp*rho); # thermal diffusivity (mˆ2/s)
Tb = 0 ; # temperatures at boundaries (o C)
A = 2.6e-6 ; # heat production (W/mˆ3)
H = A/(rho*Cp); # heat source term (o K/s) % numerical parameters
dt = 1000*seconds_per_yr ; # time step (s)
ntime = 5000 ; # number of time steps
nels = 40 ; # total number of elements
nod = 2 ; # number of nodes per element
nn = nels+1 # total number of nodes
dx = lx/nels ; # element size
g_coord = np.arange(0, lx+1, dx)#[0:dx:lx]
bcdof = np.array([1, nn]); #[ 1 nn ] ; boundary nodes
bcval = np.array([Tb, Tb]); #[ Tb Tb ] ; # boudary values
g_num = np.zeros((nod, nels), float); #zeros(nod,nels) ;
g_num[0,:]=np.arange(1, nn); #g_num(1,:) = [1:nn-1] ;
g_num[1,:]=np.arange(2, nn+1); #g_num(2,:) = [2:nn] ;
# initialise matrices and vectors
ff = np.zeros((nn,1), float); # system load vector
b = np.zeros((nn,1), float); # system rhs vector
lhs=lil_matrix((nn, nn)) #lhs = sparse(nn,nn); system lhs matrix
rhs=lil_matrix((nn, nn)) #rhs = sparse(nn,nn); system rhs matrix
displ = np.zeros((nn,1), float); # initial temperature (o C)
#-----------------------------------------------------
# matrix assembly
#-----------------------------------------------------
# Matlab version of the loop
#-----------------------------------------------------
#for iel=1:nels # loop over all elements
# num = g_num(:,iel) ; # retrieve equation number
# dx = abs(diff(g_coord(num))) ; # length of element
# MM = dx*[1/3 1/6 ; 1/6 1/3 ] ;# mass matrix
# KM = [kappa/dx -kappa/dx ; -kappa/dx kappa/dx ]; #diffn matrix
# F = dx*H*[1/2 ; 1/2] ; # load vector
# lhs(num,num) = lhs(num,num) + MM/dt + KM ; # assemble lhs
# rhs(num,num) = rhs(num,num) + MM/dt ; # assemble rhs
# ff(num) = ff(num) + F ; # assemble load
#end # end of element loop
#Python version of the loop
#-----------------------------------------------------
for iel in range(0, nels): # loop over all elements
num = g_num[:,iel] # retrieve equation number
#print(num)
dx = abs(np.diff(g_coord[num])) # length of element
MM = dx*(np.array([[1/3, 1/6],[1/6, 1/3]])) # mass matrix
KM = np.array([[kappa/dx, -kappa/dx],[-kappa/dx, kappa/dx]])
F = dx*H*(np.array([1/2, 1/2])).reshape(-1,1) # load vector
lhs[num,num] = lhs[num,num] + MM/dt + KM # assemble lhs
rhs[num,num] = rhs[num,num] + MM/dt # assemble rhs
ff[num] = ff[num] + F # assemble load
错误似乎是因为num
是一个float
。 只需做:
dx = abs(np.diff(g_coord[np.int32(num)]))
但是,它会在几行后引发另一个错误num
因为它是一个 2 元素数组。你知道代码应该做什么,而我没有。如果您有更多问题,可以在下面发表评论或在解决第一个问题后编辑您的问题。
另外,我注意到您在 Matlab 中将所有;
都留在了行尾。在 Python 中你不需要这个。另外,我认为您在创建zeros
矩阵时无需指定float
,它们自然是float
的。