使用 numpy 实现微分矩阵时出错答案

【问题标题】：Error implementing differentiating matrix using numpy使用 numpy 实现微分矩阵时出错
【发布时间】：2020-01-15 20:23:45
【问题描述】：

我正在尝试使用微分矩阵求导数。现在设法获得微分矩阵，我无法计算导数。注意我之前问过： Constructing a Multidimensional Differentiation Matrix

我目前的代码如下

import numpy as np 
import matplotlib.pyplot as plt
import math


N = 20
x = np.linspace(-1,1,N)


f1 = np.sin(3*x) # exact function 
df1 = 3*np.cos(3*x) # exact derivative for comparison sake 

def Dmatrix(N,f): 
    m_ij = np.zeros((N,N,N))
    fprime = np.zeros(N)
    for i in range(0,N-1):
        x = np.cos([(2*i + 1)/2*N])
        for j in range(0,N-1):
            for k in range(0,N-1):
                m_ij[i,j,k] = -(2/N)*((k*np.sin(k*np.pi*(2*i + 1)/2*N))*(np.cos(k*np.pi*(2*j +1))/2*N)/(np.sin(np.pi*(2*i + 1)/2*N)))
                fprime[j] += f[x[j]]*m_ij[i,j,k]
    return m_ij,fprime

dij,fprime = Dmatrix(N,f1) 

plt.plot(x,f1,'b')
plt.show()
plt.plot(x,fprime,'k')
plt.show()
plt.plot(x,df1,'r')
plt.show()

不幸的是，我得到了错误：

  File "/home/~", line 20, in Dmatrix
    fprime[j] += f[x[j]]*m_ij[i,j,k]

IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices

注意修改后的Dmatrix函数没有搭配点x，代码如下：


def Dmatrix(N,f): 
    m_ij = np.zeros((N,N,N))
    fprime = np.zeros(N)
    for i in range(N):
        for j in range(N):
            for k in range(N):
                m_ij[i,j,k] = -(2/N)*((k*np.sin(k*np.pi*(2*i + 1)/2*N))*(np.cos(k*np.pi*(2*j +1))/2*N)/(np.sin(np.pi*(2*i + 1)/2*N)))
                fprime[j] += f(x[j])*m_ij[i,j,k]
    return m_ij,fprime

这为 fprime 提供了精确解 f1 的类似图，这是错误的，因为我想获得精确解的导数，即 df1。

我不知道为什么会这样。欢迎任何想法。

【问题讨论】：

那么您检查问题行中的索引了吗？ j、x[j] 和 i,j,k 是什么？您有运行到该点的代码。我也许可以复制粘贴并自己找出来，但如果你自己检查就好了！
好的，我能够运行您的代码并进行诊断打印。由此可见，x[j] 不是有效索引！
还要注意x 的形状为(1,)，因此对于任何大于2 的N，即使进行更正也会失败

标签： python numpy for-loop matplotlib numeric

【解决方案1】：

变化

fprime[j] += f[x[j]]*m_ij[i,j,k]

到

fprime[j] += f[x[0]]*m_ij[i,j,k]

输出以下三张图片

enter image description here

不确定这是否是您要找的东西？

【讨论】：

感谢@Radonic 发现这个错误！你还做过别的吗？我已根据您的建议更改了此设置，但仍然出现相同的错误。
我得到的唯一错误是索引超出范围，因为 x 是一个元素的数组，并且在执行上述更改后，代码有效。我附上的图片是你想要的输出吗？

【解决方案2】：

当我在您的代码中添加 print 时：

            print(i,j, k, x[j])
            fprime[j] += f[x[j]]*m_ij[i,j,k]

1531:~/mypy$ python3 stack59759031.py 
0 0 0 -0.8390715290764524
Traceback (most recent call last):
  File "stack59759031.py", line 25, in <module>
    dij,fprime = Dmatrix(N,f1) 
  File "stack59759031.py", line 22, in Dmatrix
    fprime[j] += f[x[j]]*m_ij[i,j,k]
IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices

所以x[j] 是一个浮点值，它不是一个有效的索引！

【讨论】：

谢谢，我会进一步考虑。我什至不确定如何将原始函数索引到我的导数矩阵。公式表明它必须是 ``` f(x_j) *Dmatrix{i,j}``` 其中 f 是原始函数。