循环索引错误（索引100超出尺寸100的轴0的范围）

Question

我试图通过使用不同的Peclet数字来使用for循环来计算theta变量，但是当我将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

theta是一个ndarray或稀疏矩阵（请参见此处的文档：https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.spsolve.html#scipy.sparse.linalg.spsolve）。

当n = 100（for循环的第二次迭代）时，您要求theta矩阵中的第100个项目，该项目不存在。

我不太了解您要尝试做什么，但是，如果您只想遍历矩阵并获得第i个项目，则可以进行枚举for循环，如下所示：

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)