numpy 在Python中输出“nan”而不是矩阵元素?

inn6fuwd  于 11个月前  发布在  Python
关注(0)|答案(2)|浏览(131)

我试图计算这个稀疏矩阵问题Au=f的近似解,但是我在矩阵的输出中得到了“nan”而不是矩阵元素,尽管输入矩阵都包含np.float64。我不明白我错在哪里。
这是我的代码:

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

#Defining a function that takes N as an input

def our_func (N) :

  #Taking the given values of parameters h and k
  h =  1/N
  k = (29*(np.pi))/2
   
  #Defining initial input matrix with domensions N+1
  I = np.empty([N+1,N+1],  dtype=np.float64) 
  #[REF :- https://www.geeksforgeeks.org/how-to-create-an-empty-matrix-with-numpy-in-python/]
  print("datatype of I is :-")
  print(I.dtype)

  #Generating the matrix A in coo sparse format
  A_coo = coo_matrix(I).reshape((N+1,N+1))
  print("datatype of A_coo is :-")
  print(A_coo.dtype)

  #Converting the matrix A from coo sparse format to CSR sparse format 
  A = A_coo.tocsr().reshape((N+1,N+1))
  print("datatype of A is :-")
  print(A.dtype)
 
  # Creating an empty storage array for the right side of our equation Au=f
  f = np.empty((N+1), dtype=np.float64)

  #The rows of matrix A and f are defined as :-
  A[0,0] = 1    #for i==j==0
  A[N,N] = 1    #for i==j==N
  f[N] = 1      #for i==N
  f[0] = 0      #for i==0
  for i in range (int(N)) :
    for j in range (int(N)) :
        A[i,1] = 2 - ((h**2)*(k**2)) 
        A[i,i+1] = -1
        A[i,i-1] = -1
        A[i,j] = 0      #for i!==j
  
  #Our function Returns matrix A and f, with matrix A in sparse storage format
  return A, f

# Defining a new function using the function scipy.sparse.linalg.spsolve 
# to solve the sparse matrix-vector problem: Au=f

def our_solution(N):
  A,f = our_func(N)
  our_sparse_matrix = scipy.sparse.linalg.spsolve(A,f)
  return our_sparse_matrix

# Using this defined function to compute the approximate solution for our 
# problem, for N=10, N=100, and N=1000.
approx_sol_1 = our_solution(10)
print("approx_sol_1 is :- ")
print(approx_sol_1)

字符串

toiithl6

toiithl61#

在较小的N上运行函数:

In [4]: A,f = our_func(10)
datatype of I is :-
float64
datatype of A_coo is :-
float64
datatype of A is :-
float64
C:\Users\14256\miniconda3\lib\site-packages\scipy\sparse\_index.py:103: SparseEfficiencyWarning: Changing the sparsity structure of a csr_matrix is expensive. lil_matrix is more efficient.
  self._set_intXint(row, col, x.flat[0])

字符串
它警告我们,单独设置稀疏csr数组的元素是没有效率的。

In [5]: A
Out[5]: 
<11x11 sparse matrix of type '<class 'numpy.float64'>'
    with 103 stored elements in Compressed Sparse Row format>


注意,矩阵不是特别稀疏,可能的121个值中有103个是非零的。

In [6]: A.A
Out[6]: 
array([[  0.        ,  -1.        ,   0.        ,   0.        ,
          0.        ,   0.        ,   0.        ,   0.        ,
          0.        ,   0.        ,  -1.        ],
       [ -1.        , -18.75084325,  -1.        ,   0.        ,
          0.        ,   0.        ,   0.        ,   0.        ,
          0.        ,   0.        ,   0.        ],...


它可以被修剪,去掉许多0。

In [10]: A.eliminate_zeros()
In [11]: A
Out[11]: 
<11x11 sparse matrix of type '<class 'numpy.float64'>'
    with 28 stored elements in Compressed Sparse Row format>

In [17]: def my_func (N) :
    ...: 
    ...:   #Taking the given values of parameters h and k
    ...:   h =  1/N
    ...:   k = (29*(np.pi))/2
    ...: 
    ...:   #Defining initial input matrix with domensions N+1
    ...:   A = np.zeros([N+1,N+1],  dtype=np.float64)
    ...: 
    ...:   # Creating an empty storage array for the right side of our equation Au=f
    ...:   f = np.zeros((N+1), dtype=np.float64)
    ...: 
    ...:   #The rows of matrix A and f are defined as :-
    ...:   A[0,0] = 1    #for i==j==0
    ...:   A[N,N] = 1    #for i==j==N
    ...:   f[N] = 1      #for i==N
    ...:   #f[0] = 0      #for i==0
    ...:   for i in range (int(N)) :
    ...:     for j in range (int(N)) :
    ...:         A[i,1] = 2 - ((h**2)*(k**2)) 
    ...:         A[i,i+1] = -1
    ...:         A[i,i-1] = -1
    ...:         A[i,j] = 0      #for i!==j
    ...: 
    ...:   #Our function Returns matrix A and f, with matrix A in sparse storage format
    ...:   return A, f
    ...:   

In [18]: A,f = my_func(10)

In [19]: A.shape
Out[19]: (11, 11)


这可以在以下情况下转换为稀疏:

In [21]: M = coo_matrix(A).tocsr()
In [22]: M
Out[22]: 
<11x11 sparse matrix of type '<class 'numpy.float64'>'
    with 30 stored elements in Compressed Sparse Row format>


除非N太大而导致内存错误,否则我认为这种方法会更快。
这并没有解决你的nan值的问题。我对你的问题了解不够,不知道A是否有合理的值。

编辑

我想知道这是否

A[i,j] = 0      #for i!==j


如果我在my_func中忽略它,因为它已经是0,那么就会有一个不同的A[8,9]值(-1)。

In [51]: A,f = my_func(10)

In [52]: scipy.sparse.linalg.spsolve(A,f)
C:\Users\14256\miniconda3\lib\site-packages\scipy\sparse\linalg\_dsolve\linsolve.py:214: SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix format
  warn('spsolve requires A be CSC or CSR matrix format',
Out[52]: 
array([ 0.8987385, -0.1012615,  1.       ,  0.1012615,  0.8987385,
        1.797477 ,  1.       ,  0.1012615,  0.8987385,  1.797477 ,
        1.       ])


our_func使用A1

In [53]: scipy.sparse.linalg.spsolve(A1,f)
C:\Users\14256\miniconda3\lib\site-packages\scipy\sparse\linalg\_dsolve\linsolve.py:276: MatrixRankWarning: Matrix is exactly singular
  warn("Matrix is exactly singular", MatrixRankWarning)
Out[53]: array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan])


更正稀疏警告会得到相同的结果:

In [54]: scipy.sparse.linalg.spsolve(coo_matrix(A).tocsr(),f)
Out[54]: 
array([ 0.8987385, -0.1012615,  1.       ,  0.1012615,  0.8987385,
        1.797477 ,  1.       ,  0.1012615,  0.8987385,  1.797477 ,
        1.       ])


我不知道A是否正确,至少它不是单数。看起来A应该有3个对角线,'-1,?,-1'。

ttisahbt

ttisahbt2#

我不明白我错在哪里。
A矩阵是奇异的,所以你试图求解的线性系统没有唯一的解。你可以看到这一点,因为对于N=10N=100,倒数第二列的所有元素都是零(或基本为零)。例如:

A, f = our_func(10)
A.todense()[:, -2].T
# matrix([[0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
#          0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
#          0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
#          0.00000000e+000, 4.69807001e-294]])
# 4.69807001e-294 is essentially 0

字符串
sparse.linalg.spsolve遇到一个(数值上的)奇异矩阵时,它会发出警告并返回一个NaN数组,而不是引发错误。

from scipy import sparse
A = sparse.csc_matrix((3, 3))
b = np.zeros(3)
sparse.linalg.spsolve(A, b)
# array([nan, nan, nan])


若要解决此问题,请检查用于生成矩阵的代码,并确保其按预期工作。如果您确定这是您要解决的线性系统,则需要确定系统没有唯一解意味着什么。

相关问题