numpy python中带边界条件的矩阵向量乘法

3phpmpom  于 2023-11-18  发布在  Python
关注(0)|答案(1)|浏览(102)

我的任务是建立一个带边界条件的矩阵向量乘法,用迭代法求解波动方程,但我在建立这个问题上很吃力。
Sparse matrix problem
我该如何实现边界条件呢?这是我到目前为止所知道的:

def matrix_problem(N):
    """Generate the matrix and rhs associated with the wave equation."""

    k = (29 * np.pi) / 2
    h = 1 / N
    
    nelements = (2-((h**2)*(k**2))) -2

    row_ind = np.empty(nelements, dtype=np.float64)
    col_ind = np.empty(nelements, dtype=np.float64)
    data = np.empty(nelements, dtype=np.float64)

    f = np.empty(N * N, dtype=np.float64)

    count = 0
    for j in range(N):
        for i in range(N):
                if i==j:
                    row_ind[count] = 1
                    col_ind[count] = 1
                else:
                    if j==i:
                        row_ind[count] = 2-(h**2)*(k**2)
                        col_ind[count] = 2-(h**2)*(k**2)
                    if j==i+1 or j==i-1:
                        row_ind[count] = -1
                        col_ind[count] = -1
                    else:   
                        row_ind[count] = 0
                        col_ind[count] = 0
                    
    return coo_matrix((data, (row_ind, col_ind)), shape(N**2, N**2)).tocsr(), f

字符串
我试过设置矩阵向量问题,但我错误地实现了边界条件。我是新的编码,所以我真的不明白我在做什么,为什么,或者怎么会错。我试过以不同的方式实现if部分,但我遇到了同样的问题,那就是没有创建矩阵。

N = 100

A, f = discretise_poisson(N)
sol = spsolve(A, f)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/var/folders/d4/xpz4321n1k927dx91d1ntp6w0000gn/T/ipykernel_71373/3717722815.py in <module>
      1 N = 100
      2 
----> 3 A, f = matrix_problem(N)
      4 sol = spsolve(A, f)
      5 

/var/folders/d4/xpz4321n1k927dx91d1ntp6w0000gn/T/ipykernel_71373/3681280231.py in matrix_problem(N)
      7     nelements = (2-((h**2)*(k**2))) -2
      8 
----> 9     row_ind = np.empty(nelements)
     10     col_ind = np.empty(nelements)
     11     data = np.empty(nelements)

TypeError: 'float' object cannot be interpreted as an integer

的数据
不幸的是,我不知道我在做什么,所有的代码都是纯粹的猜测。我已经向我的讲师寻求帮助,但当然没有得到他们的回复,所以我不得不在这里问。一个例子会很好,这样我就可以理解它是如何设置的,所以我可以实际上破解作业(这是使用GPU加速来加速迭代求解器)。

tvmytwxo

tvmytwxo1#

关于编程,你有很多东西要学,pythonnumpyscipy.sparse,但首先要学的是调试。
我从来不只是在编辑器中编写代码,我总是有一个交互式的python会话,在那里我可以测试代码片段。
让我们测试一下nelement计算行:

In [229]: N=100; k = (29 * np.pi) / 2
     ...: h = 1 / N
     ...: nelements = (2-((h**2)*(k**2))) -2

In [230]: N,h,k,nelements
Out[230]: (100, 0.01, 45.553093477052, -0.20750843253290374)

字符串
因此,即使N=100是整数,h=1/100也是浮点数。k更大,但nelements仍然是浮点数,而且不仅仅是浮点数,0.2也很小。你不能用浮点数定义数组,更不用说小于1的数组了。更糟糕的是,它是负数!
我不知道你是从哪里得到这些计算的,但很明显他们错了。
我怀疑你从一些源代码中复制粘贴代码,而不知道发生了什么。这样做会给你留下很多调试问题。

相关问题