我试图计算这个稀疏矩阵问题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)
字符串
2条答案
按热度按时间toiithl61#
在较小的
N
上运行函数:字符串
它警告我们,单独设置稀疏
csr
数组的元素是没有效率的。型
注意,矩阵不是特别稀疏,可能的121个值中有103个是非零的。
型
它可以被修剪,去掉许多0。
型
这可以在以下情况下转换为稀疏:
型
除非
N
太大而导致内存错误,否则我认为这种方法会更快。这并没有解决你的
nan
值的问题。我对你的问题了解不够,不知道A
是否有合理的值。编辑
我想知道这是否
型
如果我在
my_func
中忽略它,因为它已经是0,那么就会有一个不同的A[8,9]
值(-1)。型
从
our_func
使用A1
:型
更正稀疏警告会得到相同的结果:
型
我不知道
A
是否正确,至少它不是单数。看起来A
应该有3个对角线,'-1,?,-1'。ttisahbt2#
我不明白我错在哪里。
A
矩阵是奇异的,所以你试图求解的线性系统没有唯一的解。你可以看到这一点,因为对于N=10
和N=100
,倒数第二列的所有元素都是零(或基本为零)。例如:字符串
当
sparse.linalg.spsolve
遇到一个(数值上的)奇异矩阵时,它会发出警告并返回一个NaN数组,而不是引发错误。型
若要解决此问题,请检查用于生成矩阵的代码,并确保其按预期工作。如果您确定这是您要解决的线性系统,则需要确定系统没有唯一解意味着什么。