循环索引错误（索引100超出大小为100的轴0的界限）

Question

我试图使用for循环通过使用不同的Peclet数来计算θ变量，但当我将Pe数设置为100时，我的循环不起作用。你知道我的代码有什么问题吗

这是我的代码

import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as spl
import matplotlib.pyplot as plt
np.set_printoptions(precision=2)

rho = 1000
c = 4000
Pe = 10000
L = 2
N = 10
deta = 2 / N
k = 0.07

def comolokko(Pe):   
    x = np.linspace(-1,1,N+1)  
    y = np.linspace(-1.25+deta/2,1.25-deta/2,N)
    y2 = np.linspace(deta/2-deta/2,deta/2-deta/2,N)
    ux = -np.cos(np.pi*x/L)
    C = sp.diags([0.5*np.ones(N),0.5*np.ones(N)],[-1,0],(N+1,N)).tocsr() 
    G = sp.diags([-np.ones(N),np.ones(N)],[-1,0],(N+1,N)).tocsr() / deta
    G[0,0] = 0    # No flux at the boundary X direction
    G[-1,-1] = 0 
    C = sp.diags(ux) @ C
    I = sp.eye(N)
    I2 = np.zeros(N)
    Iy = sp.diags(np.sin(np.pi*y/L))
    Iy2 = sp.diags(np.sin(np.pi*y2/L))
    #print(Iy2.todense())
    F = sp.kron(Iy,Pe * C) -sp.kron(I,G) #flux operator for all vertical faces
    F2 = sp.kron(Iy2,Pe * C) -sp.kron(I,G)
    D = sp.diags([-np.ones(N),np.ones(N)],[0,1],(N,N+1)).tocsr()
    Ax = sp.kron(I,D) @ F
    Ax2 = sp.kron(I,D) @ F2
    ####y direc
    G = sp.diags([-np.ones(N),np.ones(N)],[-1,0],(N+1,N)).tocsr() / deta 
    G[0,0] *= 2    # Diriclet at boundary   cuz 2l ?
    G[-1,-1] *= 2  # Diriclet at boundary
    C = sp.diags([0.5*np.ones(N),0.5*np.ones(N)],[-1,0],(N+1,N)).tocsr()
    # Coordinates of the horizontal face centers
    y = np.linspace(-1,1,N+1)
    y2 = np.linspace(0,0,N+1)
    #print(y2)
    x = np.linspace(-1.25+deta/2,1.25-deta/2,N)
    x2 = np.linspace(deta/2-deta/2,deta/2-deta/2,N)
    uy = np.cos(np.pi*y/2) # Velocity part dependent on y only
    #print(uy)

    C = sp.diags(uy) @ C    # The operator u_y * phi on the faces
    C2 = sp.diags(y2) @ C
                            # to 2D
    I = sp.eye(N)
    Ix = sp.diags(np.sin(np.pi*x/2))  
    Ix3=  sp.diags(np.sin(np.pi*x2/2))
    F = sp.kron(Pe * C,Ix) - sp.kron(G,I)
    F3 = sp.kron(Pe * C2,Ix3) - sp.kron(G,I) # Flux operator for all the horizontal faces
    # Conservation operator 
    D = sp.diags([-np.ones(N),np.ones(N)],[0,1],(N,N+1)).tocsr()
    Ay = sp.kron(D,I) @ F
    Ay2 = sp.kron(D,I) @ F3
    A3 = Ay2
    A = Ax + Ay
    A2 = Ax2 + Ay2
    g = np.zeros(N)
    g[0] = G[0,0]
    b = np.kron(g,np.ones(N))
    theta = spl.spsolve(A,b)
    theta2= spl.spsolve(A2,b)
    theta3= spl.spsolve(A3,b)
    return theta, theta2

代码的第一部分工作正常，但当我使用for循环时，代码不工作。此外，我还可以手动更改comolokko定义函数的Pe编号，而不会出现任何索引错误

thetaq = list()
thetaq2 = list()
Pe = [10, 100, 1000, 10000]
for n in Pe:
    theta, theta2 = comolokko(n)
    thetaq.append(theta[int(n)])
    theta, theta2 = comolokko(n)
    thetaq2.append(theta2[int(n)])
thetaq = np.array(thetaq)
thetaq2 = np.array(thetaq2)
plt.semilogy(Pe, ((thetaq/thetaq2)*2))
plt.xlabel('Peclet Number')
plt.ylabel('q/q_0')
plt.grid()
plt.show()

有此错误

IndexError                                Traceback (most recent call last)
<ipython-input-38-3b63413cc0ab> in <module>
      4 for n in Pe:
      5     theta, theta2 = comolokko(n)
----> 6     thetaq.append(theta[int(n)])
      7     theta, theta2 = comolokko(n)
      8     thetaq2.append(theta2[int(n/2)])

IndexError: index 100 is out of bounds for axis 0 with size 100

Answer 1

θ是一个数组或稀疏矩阵（参见文档：https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.spsolve.html#scipy.sparse.linalg.spsolve）

当n=100（for循环的第二次迭代）时，您需要θ矩阵中不存在的第100项

我真的不明白你到底想做什么，但是，如果你只是想在矩阵中循环并得到第I项，你可以对循环进行枚举，如下所示：

for i, n in enumerate (Pe):
   #you code
   #where i is the i-th iteration (like a counter, i runs from 0 to len(Pe) - 1
   #and n is the Pe itself as in the list (10, 100, etc)