Python - Google Colab - IndexError：索引 2 超出尺寸为 1 的轴 0 的范围

Question

我正在尝试求解稀疏矩阵方程

u_derivative_1 = A * u

（A为稀疏矩阵）

但是我收到以下错误：-

  IndexError                                Traceback (most recent call last)
  <ipython-input-24-f4af80e4ae52> in <cell line: 1>()
  ----> 1 trial1 = discretise_delta_u_v4(1000, 'implicit')
  <ipython-input-23-731d13e4ddf7> in discretise_delta_u_v4(N, method)
       61         for i in range (1 , N-1):
       62           for j in range (1 , N-1):
  ---> 63             A[i,j] = (u[i-1,j] + u[i+1,j] + u[i,j-1] + u[i,j+1] - (4*u[i,j]))/(h**2)
       64 
  IndexError: index 2 is out of bounds for axis 0 with size 1

我很困惑为什么会收到此错误以及如何解决此错误。这是我的代码 -

import numpy as np
import scipy
import scipy.sparse
from scipy.sparse import csr_matrix
from scipy.sparse import coo_matrix

def discretise_delta_u_v4(N, method):

    i = np.arange(0,N)
    j = np.arange(0,N)
    h = 2/N

    A = csr_matrix((N, N), dtype = np.float64).toarray()
    u = np.array([[(i*h), (j*h)]])

    #u[ih,jh] = 
    u[:,0] = 5    #Boundary
    u[:,-1] = 0   #Boundary
    u[0,:] = 0    #Boundary
    u[-1,:] = 0   #Boundary

    #Implicit
    if (method=='implicit') :
        A[0,:] = 0
        A[-1,:] = 0
        A[:,0] = 0
        A[:,-1] = 0
        for i in range (1 , N-1):
          for j in range (1 , N-1):
            A[i,j] = (u[i-1,j] + u[i+1,j] + u[i,j-1] + u[i,j+1] - (4*u[i,j]))/(h**2)

    # u_der_1 = A * u
    for i in range (0 , N):
      for j in range (0 , N):
        u_der_1 = scipy.sparse.linalg.spsolve(A,u)


trial1 = discretise_delta_u_v4(1000, 'implicit')

Answer 1

您具体看到的错误是由构建数组的方式引起的

u

。

做

u = np.array([[(i*h), (j*h)]])

会导致

u

：

[[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.], [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]]

这是一个形状为 (1, 2, 10) 的数组。因此，当您执行

u[i+1,j]

时，

i+1

的唯一有效值为 0，其他所有内容都超出范围。

鉴于您试图解决的问题，该方法还存在一些问题。

1. 您希望 A 是一个稀疏矩阵，但您通过将其转换为稠密数组来将其创建为稠密矩阵（即

   A = csr_matrix((N, N), dtype = np.float64).toarray()

   ）。您应该删除对

   .toarray()

   的呼叫。

2. 您还错误地初始化了

   u

   ，您想要一个由

   N**2

   元素组成的一维数组，而不是初始化一个二维数组。

3. 循环中有

   u_der_1 = scipy.sparse.linalg.spsolve(A,u)

   ，这也是不正确的。您需要先构建矩阵 A，然后再求解系统。

def discretise_delta_u_v4(N, method):
    
    # Initialise the matrix A as a list of lists first
    h = 2 / N
    A = list_of_lists_matrix((N ** 2, N ** 2), dtype=np.float64)

    if method == 'implicit':
        for i in range(N):
            for j in range(N):
                index = i * N + j
                if i == 0 or i == N - 1 or j == 0 or j == N - 1:
                    A[index, index] = 1
                else:
                    A[index, index] = -4 / (h ** 2)
                    A[index, index - 1] = 1 / (h ** 2)
                    A[index, index + 1] = 1 / (h ** 2)
                    A[index, index - N] = 1 / (h ** 2)
                    A[index, index + N] = 1 / (h ** 2)

    A = A.tocsr(). # You conver to a CSR matrix here

    u = np.zeros(N ** 2)
    u[0:N] = 5  # Per your boundary condition

    u_der_1 = scipy.sparse.linalg.spsolve(A, u_flat)
    return u_der_1.reshape(N, N)

trial1 = discretise_delta_u_v4(10, 'implicit')