Ошибка индекса зацикливания (индекс 100 выходит за пределы оси 0 с размером 100)

Question

Я пытался использовать для l oop для вычисления тэта-переменных, используя разные числа Пеклета, но когда я установил число Pe на 100, мой l oop не работает. У вас есть идеи, что не так с моим кодом?

Вот мой код

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

первая часть кода работает, но когда я использую для Код l oop не работает. Кроме того, я могу вручную изменить число Pe для определенной функции comolokko без какой-либо ошибки индекса.

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 или разреженная (см. Do c здесь: https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.spsolve.html#scipy .sparse.linalg.spsolve ).

Когда n = 100 (вторая итерация для l oop), вы запрашиваете сотый элемент в тэта-матрице, которого не существует.

Я не совсем понимаю, что вы пытались сделать, но, если вы просто хотите пройти l oop через матрицу и получить i-тый элемент tge, вы можете выполнить перечисление для l oop, например: это:

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)