如何解决python上的RBF-FD拉普拉斯求解器
我正在尝试使用 RBF-FD 方法显式地找到函数的拉普拉斯算子 功能是
def func(x,z):
return np.sin(np.pi*x)*np.cos(np.pi*z/2)
f = func(points[:,0],points[:,1]).reshape(nx,nz)
其拉普拉斯算子的解析解为
def laplace_abs(x,z):
return -5/4*(np.pi**2)*np.sin(np.pi*x)*np.cos(np.pi*z/2)
laplace_abs = laplace_abs(points[:,nz)
参数为:
nx = 101
nz = nx
dx = 0.01
dz = dx
x = np.linspace(0,1,nx)
z = np.linspace(0,nz)
d2px = np.zeros((nz,nx))
d2pz = np.zeros((nz,nx))
X,Z = np.meshgrid(x,z)
points = np.column_stack([X.ravel(),Z.ravel()])
stencil = 5
size = len(points)
然后,我试图为 FDM 中的 5 点模板找到 RBF-FD 系数。找到系数后,拉普拉斯值应表示为加权和。相关公式如下:
代码的主要部分是:
def get_fd_indices():
nbrs = NearestNeighbors(n_neighbors=stencil,algorithm='ball_tree',metric='euclidean').fit(points)
return nbrs.kneighbors(points)
distances,indices = get_fd_indices()
r = distances
RHS = np.where(distances != 0,r**2*(16*np.log(r)+8),0)
def LHS(x,xi):
R = np.zeros((size,stencil,stencil))
for i in range(size):
R[i] = distance.cdist(x[i],x[i],'euclidean')
LHS = np.where(R != 0,R**4*np.log(R),0)
return LHS
LHS = LHS(points[indices],points[indices])
def get_coef(LHS,RHS):
c = np.zeros((size,1))
for i in range(size):
Left = LHS[i]
Right = RHS[i].reshape((stencil,1))
c[i] = LA.pinv(Left).dot(Right)
return c
c = get_coef(LHS,RHS)
values = f.ravel()
laplaceRBF = np.zeros(size)
for i in range(size):
index = indices[i]
laplaceRBF[i] = np.sum(values[index]*c[i])
laplaceRBF = laplaceRBF.reshape(nx,nz)
我希望获得类似于 FDM 解决方案的结果
def FDM(f,dx):
for i in range(1,nx - 1):
d2px[i,:] = (f[i - 1,:] - 2 * f[i,:] + f[i + 1,:]) / dx ** 2
for j in range(1,nz - 1):
d2pz[:,j] = (f[:,j - 1] - 2 * f[:,j] + f[:,j + 1]) / dz ** 2
return d2px+d2pz
FDM = FDM(f,dx)
现在 RBF-FD 的结果与 FDM 完全不同。 RBF-FD 实现中的一些错误。请告诉我,出了什么问题。
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。