python 矩阵运算
本帖最后由 lidaxin 于 2020-6-13 10:33 编辑import numpy as np# 导入numpy模块
import matplotlib.pyplot as plt# 导入matplotlib库
timeout = 0.0045# physical time to compute,s
length = 10.0# length of the doamin in x direction,m
nnod = 101# number of nodes
dx = length/(nnod-1)# space interval 步长
X = np.mat(np.arange(0,length+dx,dx)) # x-coordinates of nods
T0 = 300.0# initial temperature,K
P0L = 10.0# left initial pressure,bar
P0R = 1.0 # right initial pressure,bar
v0 = 0.0# initial velocity,m/s
gamm = 1.4# Cp/Cv
R = 287.0# gas constant,J/kg.K
CFL= 0.98
# 1bar=100000Pa
D0L = 100000.0*P0L/(R*T0)# left density,kg/m^3
D0R = 100000.0*P0R/(R*T0)# right density,kg/m^3
E0 = R*T0/(gamm-1.0)# initial total energy,J/kg
H0 = E0+R*T0# initial enthalpy,J/kg
c0 = (gamm*R*T0)**0.5# initial sound speed,m/s
Den = np.mat(np.zeros((1,nnod))) # density at nodes
Tem = np.mat(np.zeros((1,nnod)))# temperature at nodes
Pre = np.mat(np.zeros((1,nnod)))# pressure at nodes
Vel = np.mat(np.zeros((1,nnod)))# velocity at nodes
M = np.mat(np.zeros((1,nnod)))# mach number at nodes
Ene = np.mat(np.zeros((1,nnod)))# total energy at nodes
Ent = np.mat(np.zeros((1,nnod)))# enthalpy at nodes
c = np.mat(np.zeros((1,nnod)))# sound speed at nodes
dt = CFL*dx/c0# initial CFL condition
# construction and initialization of U matrix
matL = np.mat(np.ones((1,int((nnod-1)/2))))
matR = np.mat(np.ones((1,int((nnod-1)/2+1))))
U1 = np.hstack((matL*D0L,matR*D0R))# rho
U2 = np.mat(np.zeros((1,nnod)))# rho*u
U3 = np.hstack((matL*D0L*E0,matR*D0R*E0))# rho*E
Un = np.array()# Un
Un = np.mat(Un)
Un1 = Un
Unm = Un1
print(Unm)
# construction and initialization of F matrix
F1 = np.mat(np.zeros((1,nnod)))
F2 = np.hstack((matL*P0L,matR*P0R))# rho*u^2+p
F3 = np.mat(np.zeros((1,nnod)))# rho*u*h
Fn = np.array()
Fn = np.mat(Fn)
a = Fn-Fn
b = Un-a
print(b)
print(Fn)
time = 0.0
nlter = 1
while (time<timeout):
time = time+dt
for j in range(0,nnod):
if (j==nnod-1):
Unm = Un
print(Unm)
elif (j >= 0)and(j <= nnod-2):
Unm = Un-(dt/dx)*(Fn-Fn)
Den = Unm
Vel = np.divide(Unm,Den)
Ene = np.divide(Unm,Den)
Tem = np.subtract(Ene,0.5*np.multiply(Vel,Vel))*((gamm-1.0)/R)
Pre = np.multiply(R*Den,Tem)
Ent = np.add(Ene,R*Tem)
c = np.sqrt(gamm*R*Tem)
M = np.divide(Vel,c)
Fn_ = np.array((np.multiply(Den,Vel),np.add(np.multiply(np.multiply(Den,Vel),Vel),Pre),np.multiply(np.multiply(Den,Vel),Ent)))
Fn_ = np.mat(Fn)
Un1 = 0.5*np.add(Un,Unm)-0.5*(dt/dx)*(np.subtract(Fn_,Fn))
Den = Un1
Vel = np.divide(Un1,Den)
Ene = np.divide(Un1,Den)
Tem = np.subtract(Ene,0.5*np.multiply(Vel,Vel))*((gamm-1.0)/R)
Pre = np.multiply(R*Den,Tem)
Ent = np.add(Ene,R*Tem)
c = np.sqrt(gamm*R*Tem)
M = np.divide(Vel,c)
dt = CFL*dx/np.max(np.add(c,Vel)) # calculate minimum time step
nlter += 1
Un = Un1 # update Un and Fn
Fn = np.array((np.multiply(Den,Vel),np.add(np.multiply(np.multiply(Den,Vel),Vel),Pre),np.multiply(np.multiply(Den,Vel),Ent)))
Fn = np.mat(Fn)
else:
print(Den)
print(Vel)
print(Pre)
print(Tem)
print(M)
错误：
Traceback (most recent call last):
File "D:/py/python/lianxi/lianxi/Mac.py", line 84, in <module>
Unm = Un-(dt/dx)*(Fn-Fn)
File "D:\py\python\lib\site-packages\numpy\matrixlib\defmatrix.py", line 195, in __getitem__
out = N.ndarray.__getitem__(self, index)
IndexError: index 1 is out of bounds for axis 1 with size 1
求问大佬，这是哪里出了问题？
这行报错，索引超出 1轴的范围：
Unm = Un - (dt / dx) * (Fn - Fn)
Twilight6 发表于 2020-6-13 11:01
这行报错，索引超出 1轴的范围：
我怎么感觉没超出？有什么办法解决吗？
lidaxin 发表于 2020-6-13 14:11
我怎么感觉没超出？有什么办法解决吗？
没办法只能你自己看看呗...不是自己的代码看上去感觉有点晕，而且我还不怎么娴熟 numpy 模块 Twilight6 发表于 2020-6-13 14:14
没办法只能你自己看看呗...不是自己的代码看上去感觉有点晕，而且我还不怎么娴熟 numpy 模块
好的好的，谢谢
页:
[1]