2da*_*aaa 2 python numpy scipy
我正在尝试使用scipy.ndimage.convolve计算二维场A的拉普拉斯算子。
stencil = numpy.array([[0, 1, 0],[1, -4, 1], [0, 1, 0]])
scipy.ndimage.convolve(A, stencil, mode='wrap')
Run Code Online (Sandbox Code Playgroud)
但这似乎并没有给我正确的答案。有什么想法我出错了,或者有更好的方法在 numpy 中计算拉普拉斯吗?
5年前我想知道还有人关心...
我认为这是因为卷积方法只是一种近似,模板基本上是通过有限离散微分来近似二阶导数。增量越小,数值近似值就越接近。
我根据别人的信息做了一些测试。下面是代码和结果图:
import numpy
import scipy.ndimage.filters as filters
import scipy.signal as signal
import matplotlib.pyplot as plt
stencil=numpy.array([[0,1,0],[1,-4,1],[0,1,0]])
fig=plt.figure(figsize=(10,10),dpi=100)
for ii,jj in enumerate([10,100,1000,10000]):
x=numpy.linspace(-5,5,jj)
xx,yy=numpy.meshgrid(x,x)
step=x[1]-x[0]
image=numpy.exp(-xx**2-yy**2)
lap1=filters.laplace(image)/step**2
#lap2=filters.convolve(image,stencil,mode='wrap')/step**2
#lap3=signal.convolve2d(image,stencil,mode='same')/step**2
lap4=4*image*(xx**2+yy**2)-4*image
ax=fig.add_subplot(2,2,ii+1)
img=ax.imshow(lap1-lap4)
ax.set_title('stencil - analytical (dx=%.4f)' %step)
plt.colorbar(img)
fig.tight_layout()
plt.show()
Run Code Online (Sandbox Code Playgroud)
、filters.laplace()和filters.convolve()都signal.convolve2d()给出了非常接近的结果(事实上,如果你查看filters.laplace()的源代码,它所做的事情与卷积模板内核本质上是相同的),所以我只包括一个filters.laplace()。请注意,以上所有内容都没有除以步长的平方。
该图显示增量步长越小,越接近解析解(即4z(x^2+y^2)-4z)。
| 归档时间: |
|
| 查看次数: |
7208 次 |
| 最近记录: |