誰能解釋一下如何解決這個索引錯誤-有解無憂

目前我正在使用此代碼輸出二維隱式對流擴散矩陣。我無法弄清楚如何填充列的最后一行。我讓它從一個 txt 檔案中獲取所有這些資料，該檔案也包含在代碼下面。然后我嘗試使用范圍值讓它輸出一個矩陣，但我似乎無法讓矩陣正確輸出。任何有關為什么它不起作用的幫助或解釋將不勝感激。當我將 i 和 j 的范圍更改為 Nx 或 Ny 時，我收到一個索引錯誤，我也不太明白，但是為了解決這個問題，我將范圍設定為 Nx-1 和 Ny-1，我相信不正確。再次非常感謝任何建議。

import os
import numpy as np
#import data into variables for use

file = open("diffsolverin.txt","r")
line = file.readline() # rea the first line (remove comments from input)
tokens = line.split() # this splits the line into chunks of characters separated by spaces
Nx = int(tokens[0]) # assign the value of the first chunk of characters to Nx after interpreting them as an integer
Ny = int(tokens[1])
width = float(tokens[2])
height = float(tokens[3])
D = float(tokens[4])
t = float(tokens[5])
U_x = float(tokens[6])
U_y = float(tokens[7])
S_p = float(tokens[8])
line = file.readline() # read the second line
tokens = line.split()
left_bc_type = int(tokens[0]) #0 is Dirichlet, 1 is Neumann
right_bc_type = int(tokens[1])
top_bc_type = int(tokens[2])
bottom_bc_type = int(tokens[3])
line = file.readline() # read the third line
tokens = line.split()
left_bc_value = float(tokens[0]) 
right_bc_value = float(tokens[1])
top_bc_value = float(tokens[2])
bottom_bc_value = float(tokens[3])
        
        
N=Nx*Ny
Amat = np.zeros((N,N))
dx = Nx/width
dy = Ny/height
B = np.zeros(N)
phi = np.zeros(N)

       
#solve and establish boundary conditions

for i in range(Nx-1):
   
    for j in range(Ny-1):
        k = (j)*(Ny)   i
        print (k)

        if i == 0:    # this is a left boundary cell, deal with it given the BCs)
            if left_bc_type == 0:
                
                A_e = D*dx-U_x/2
                A_w = 0
                A_n = D*dy U_y/2
                A_s = D*dy-U_y/2
                A_p = A_e   A_w   A_s   A_n
                Amat[k,k] = A_p
                Amat[k,k-Nx]=A_n
                Amat[k,k Nx]=A_s
                Amat[k,k 1]=A_e
                B[k] = A_p*left_bc_value
            elif left_bc_type == 1:
                
                A_e = D*dx-U_x/2
                A_w = 0
                A_n = D*dy U_y/2
                A_s = D*dy-U_y/2
                A_p = A_e   A_w   A_s   A_n
                Amat[k,k] = A_p
                Amat[k,k-Nx]=A_n
                Amat[k,k Nx]=A_s
                Amat[k,k 1]=A_e
                B[k] = A_p
        elif i == Nx-1: #this is a right boundary cell, ...
            if right_bc_type == 0:
                
                A_e = 0
                A_w = D*dx U_x/2
                A_n = D*dy U_y/2
                A_s = D*dy-U_y/2
                A_p = A_e   A_w   A_s   A_n
                Amat[k,k] = A_p
                Amat[k,k-Nx]=A_n
                Amat[k,k Nx]=A_s
                Amat[k,k-1]=A_w
                B[k] = A_p*right_bc_value
                
            elif right_bc_type == 1:
                
                A_e = 0
                A_w = D*dx U_x/2
                A_n = D*dy U_y/2
                A_s = D*dy-U_y/2
                A_p = A_e   A_w   A_s   A_n
                Amat[k,k] = A_p
                Amat[k,k-Nx]=A_n
                Amat[k,k Nx]=A_s
                Amat[k,k-1]=A_w
                B[k] = A_p

        elif j == 0:  # this is a top boundary cell
               if top_bc_type == 0:
                   
                   A_e = D*dx-U_x/2
                   A_w = D*dy U_x/2
                   A_n = 0
                   A_s = D*dy-U_y/2
                   A_p = A_e   A_w   A_s   A_n
                   Amat[k,k] = A_p
                   Amat[k,k 1]=A_e
                   Amat[k,k Nx]=A_s
                   Amat[k,k-1]=A_w
                   B[k] = A_p*top_bc_value
                   
               elif top_bc_type == 1:
                   
                   A_e = D*dx-U_x/2
                   A_w = D*dx U_x/2
                   A_n = 0
                   A_s = D*dy-U_y/2
                   A_p = A_e   A_w   A_s   A_n
                   Amat[k,k] = A_p
                   Amat[k,k-1]=A_w
                   Amat[k,k Nx]=A_s
                   Amat[k,k 1]=A_e
                   B[k] = A_p

        elif j == Ny-1: # this is a bottom boundary cell
               if bottom_bc_type == 0:
                   
                   A_e = D*dx-U_x/2
                   A_w = D*dx U_x/2
                   A_n = D*dy U_y/2
                   A_s = 0
                   A_p = A_e   A_w   A_s   A_n
                   Amat[k,k] = A_p
                   Amat[k,k 1]=A_e
                   Amat[k,k-Nx]=A_n
                   Amat[k,k-1]=A_w
                   B[k] = A_p*bottom_bc_value
               elif bottom_bc_type == 1:
                   
                   A_e = D*dx-U_x/2
                   A_w = D*dx U_x/2
                   A_n = D*dy U_y/2
                   A_s = 0
                   A_p = A_e   A_w   A_s   A_n
                   Amat[k,k] = A_p
                   Amat[k,k 1]=A_e
                   Amat[k,k-Nx]=A_n
                   Amat[k,k-1]=A_w
                   B[k] = A_p
        else:
            A_e = D*dx-U_x/2
            A_w = D*dx U_x/2
            A_n = D*dy U_y/2
            A_s = D*dy-U_y/2
            A_p = A_e   A_w   A_s   A_n
            Amat[k,k] = A_p
            Amat[k,k 1]=A_e
            Amat[k,k-Nx]=A_n
            Amat[k,k-1]=A_w
            Amat[k,k Nx]=A_s
            B[k]=A_p
           
print (Amat)
#print(B)
#print (ans)
file = open("diffsolverout.txt", "w ")
filew = str(Amat)
filew1=str(B)
file.write(filew)
file.write(filew1)
file.close()

這是文本檔案的內容：

10 10 1.0 2.0 0.5 0.1 1 4 1
1 1 1 1
10 10 10 10

我再一次得到矩陣的最后一行空白，我認為這是不正確的，或者我因嘗試填充矩陣而得到索引錯誤。

uj5u.com熱心網友回復：

正如錯誤所暗示的那樣，您正在到達矩陣的邊界之外。不過，這不是來自價值i或j錯誤。你的邊界單元定義采用k Nxor k - Nx，但是當你檢查它時它是錯誤的，因為索引是如何為 python 作業的。

錯誤發生在時k = 90，您有Nx = 10，所以k Nx = 90 10 = 100。但是索引來自0-99not 1-100。您應該能夠將其中的每一個更改為k (Nx - 1)或k - (Nx - 1)并解決此特定問題。

你在說那個的時候已經能夠識別出這一點j == Ny - 1，這是相同的邏輯。

uj5u.com熱心網友回復：

好的，在進一步研究之后，您索引的Amat方式對于您嘗試做的事情是不正確的。您必須了解在 Python 中（更一般地在數學中）矩陣中的索引是如何作業的。我認為它值得一個新的解決方案，因為我想正確解釋它。

有i行j數和列數。您的方法是用整數從左到右、從上到下標記矩陣中的每個單元格k。這很好，并且對這種方法有很大幫助，但您不能直接使用它來改變矩陣中的值。（作為旁注，i是行j還是列，您的方法雖然設定為這種方式，但實際上在您定義時將其i作為列和j行k）

想象一下，我們有k = 5. 然后在矩陣中，我們將A[5,5]要處理的點作為中心點。向左或向右我們得到A[k, k-1]和A[k, k 1]，但向上或向下我們必須做A[k-1,k]和A[k 1,k]。

對于您的 else 塊（以顯示所有示例），這實際上應該是：

Amat[k, k] = A_p
Amat[k, k   1] = A_e
Amat[k - 1, k] = A_n
Amat[k, k - 1] = A_w
Amat[k   1, k] = A_s

您會注意到這只是k作為參考點，但我們仍然需要使用標準的索引方式來達到您想要的值，而不僅僅是添加或減去Nx單元格編號。如果這是有道理的。

最后，我只想添加一些建議，使您的代碼更具可讀性和可除錯性。A_e, A_w, A_n,的值A_s可以在開始時只定義一次，然后A_p可以使用特定情況所需的值進行定義。

轉載請註明出處，本文鏈接：https://www.uj5u.com/gongcheng/355126.html

標籤：Python 麻木的

上一篇：ValueError：在Pandas中移動一列時無法從重復軸重新索引

下一篇：如何遍歷和比較資料幀的值？