python 矩阵运算

lidaxin

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([U1,U2,U3])  # 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([F1,F2,F3])
Fn = np.mat(Fn)
a = Fn[0:3,1]-Fn[0:3,0]
b = Un[0:3,0]-a
print(b)
print(Fn[0:3,0])


time = 0.0
nlter = 1
while (time<timeout):
   
    time = time+dt

    for j in range(0,nnod):
        if (j==nnod-1):
            Unm[0:3,j] = Un[0:3,j]
            print(Unm)
            
        elif (j >= 0)and(j <= nnod-2):
      
            Unm[0:3,j] = Un[0:3,j]-(dt/dx)*(Fn[0:3,j+1]-Fn[0:3,j])

        Den = Unm[0,0:nnod]
        Vel = np.divide(Unm[1,0:nnod],Den)  
        Ene = np.divide(Unm[2,0:nnod],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[0:3,j],Unm[0:3,j])-0.5*(dt/dx)*(np.subtract(Fn_[0:3,j],Fn[0:3,j-1]))

        Den = Un1[0,0:nnod]
        Vel = np.divide(Un1[1,0:nnod],Den)  
        Ene = np.divide(Un1[2,0:nnod],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[0:3,j] = Un[0:3,j]-(dt/dx)*(Fn[0:3,j+1]-Fn[0:3,j])
  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

求问大佬，这是哪里出了问题？

Twilight6

这行报错，索引超出 1轴的范围：

Unm[0:3, j] = Un[0:3, j] - (dt / dx) * (Fn[0:3, j + 1] - Fn[0:3, j])

lidaxin

Twilight6 发表于 2020-6-13 11:01
这行报错，索引超出 1轴的范围：

我怎么感觉没超出？有什么办法解决吗？

Twilight6

lidaxin 发表于 2020-6-13 14:11
我怎么感觉没超出？有什么办法解决吗？

没办法  只能你自己看看呗...不是自己的代码看上去感觉有点晕，而且我还不怎么娴熟 numpy 模块

lidaxin

Twilight6 发表于 2020-6-13 14:14
没办法  只能你自己看看呗...不是自己的代码看上去感觉有点晕，而且我还不怎么娴熟 numpy 模块

好的好的，谢谢