我目前正在从事一个涉及计算电场和 3D 空间中的梯度的项目。这需要对拉普拉斯方程进行数值求解，我已经编写了一个类来执行此操作（如下），它有效，但这里仅用于背景；
################################################################################
# class: ThreeDRectLaplaceSolver #
# #
# A class to solve the Laplace equation given the potential on the boundaries.#
################################################################################
class ThreeDCuboidLaplaceSolver:
#############################################################################
# Store the init variables as class varibales, calculate the grid spacing to#
# solve Laplace's equation on and create arrays to store the iterations and #
# results in. #
############################################################################
def __init__(self, xmin, xmax, ymin, ymax, zmin, zmax, gridstep):
self.xmin, self.xmax = xmin, xmax
self.ymin, self.ymax = ymin, ymax
self.zmin, self.zmax = zmin, zmax
self.xpoints = int((xmax-xmin)/gridstep) + 1
self.ypoints = int((ymax-ymin)/gridstep) + 1
self.zpoints = int((zmax-zmin)/gridstep) + 1
self.dx = (xmax-xmin)/self.xpoints
self.dy = (ymax-ymin)/self.ypoints
self.dz = (zmax-zmin)/self.zpoints
self.v = np.zeros((self.xpoints, self.ypoints, self.zpoints))
self.old_v = self.v.copy()
self.timeStep = 0
############################################################################
# Set constant values along the boundaries #
# #
# Top(bottom) is +ive(-ive) end of z-axis #
# Right(left) is +ive(-ive) end of y-axis #
# Front(back) is +ive(-ive) end of x-axis #
############################################################################
def setBC(self, frontBC, backBC, rightBC, leftBC, topBC, bottomBC):
self.v[-1, :, :] = frontBC
self.v[0 , :, :] = backBC
self.v[: ,-1, :] = rightBC
self.v[: , 0, :] = leftBC
self.v[: , :,-1] = topBC
self.v[: , :, 0] = bottomBC
self.old_v = self.v.copy()
def solve_slow(self, PercentageTolerance = 5):
PercentageError = PercentageTolerance + 1
while PercentageError > PercentageTolerance:
self.Iterate()
PercentageError = self.Get_LargestPercentageError()
print "Completed iteration number %s \n Percentage Error is %s\n" % (str(self.timeStep), str(PercentageError))
return self.v
def solve_quick(self, Tolerance = 2):
AbsError = Tolerance + 1
while AbsError > Tolerance:
self.Iterate()
AbsError = self.Get_LargestAbsError()
print "Completed iteration number %s \nAbsolute Error is %s\n" % (str(self.timeStep), str(AbsError))
return self.v
def Get_LargestAbsError(self):
return np.sqrt((self.v - self.old_v)**2).max()
def Get_LargestPercentageError(self):
AbsDiff = (np.sqrt((self.v - self.old_v)**2)).flatten()
v = self.v.flatten()
vLength = len(v)
Errors = []
i=0
while i < vLength:
if v[i]==0 and AbsDiff[i]==0:
Errors.append(0)
elif v[i]==0 and AbsDiff[i]!=0:
Errors.append(np.infty)
else:
Errors.append(AbsDiff[i]/v[i])
i+=1
return max(Errors)*100
# Perform one round of iteration (ie the value at each point is iterated by one timestep)
def Iterate(self):
self.old_v = self.v.copy()
print self.Get_vAt(0,5,0)
self.v[1:-1,1:-1,1:-1] = (1/26)*(self.v[0:-2, 2:, 2: ] + self.v[0:-2, 1:-1, 2: ] + self.v[0:-2, 0:-2, 2: ] +\
self.v[1:-1, 2:, 2: ] + self.v[1:-1, 1:-1, 2: ] + self.v[1:-1, 0:-2, 2: ] +\
self.v[2: , 2:, 2: ] + self.v[2: , 1:-1, 2: ] + self.v[2: , 0:-2, 2: ] +\
self.v[0:-2, 2:, 1:-1] + self.v[0:-2, 1:-1, 1:-1] + self.v[0:-2, 0:-2, 1:-1] +\
self.v[1:-1, 2:, 1:-1] + self.v[1:-1, 0:-2, 1:-1] +\
self.v[2: , 2:, 1:-1] + self.v[2: , 1:-1, 1:-1] + self.v[2: , 0:-2, 1:-1] +\
self.v[0:-2, 2:, 0:-2] + self.v[0:-2, 1:-1, 0:-2] + self.v[0:-2, 0:-2, 0:-2] +\
self.v[1:-1, 2:, 0:-2] + self.v[1:-1, 1:-1, 0:-2] + self.v[1:-1, 0:-2, 0:-2] +\
self.v[2: , 2:, 0:-2] + self.v[2: , 1:-1, 0:-2] + self.v[2: , 0:-2, 0:-2])
self.timeStep += 1
# Iterate through a certain number of time steps
def IterateSteps(self, timeSteps):
i = 0
while i < timeSteps:
self.Iterate()
i+=1
# Get the value of v at a point (entered as coordinates, NOT indices)
def Get_vAt(self, xPoint, yPoint, zPoint):
# Find the indices nearest to the coordinates entered
diff = [np.sqrt((x-xPoint)**2) for x in np.linspace(self.xmin,self.xmax,self.xpoints)]
xIndex = diff.index(min(diff))
diff = [np.sqrt((y-yPoint)**2) for y in np.linspace(self.ymin,self.ymax,self.ypoints)]
yIndex = diff.index(min(diff))
diff = [np.sqrt((z-zPoint)**2) for z in np.linspace(self.zmin,self.zmax,self.zpoints)]
zIndex = diff.index(min(diff))
# retun the value from of v at this point
return self.v[xIndex, yIndex, zIndex]
所以当我运行以下
solver = ThreeDCuboidLaplaceSolver(0, 20, 0, 10, 0, 20, 1)
TwoDsolver = TwoDRectLaplaceSolver(0,20,0,10,1)
TwoDsolver.setBC(1000,0,0,0)
v_0 = np.zeros((21,11))
v_0[4:6,4:6] = 1000
v_0[14:15, 9:10] = 2000
solver.setBC(0,0,0,0,0,v_0)
v = solver.solve_quick(0.00001)
我得到一个 3D numpy 数组，其中包含网格上每个点的电位 (V) 值。然而，为了用这个做更多有用的事情，我希望能够用一个连续函数来近似这个值数组，这样我就可以计算不在我的网格上的点的场及其梯度。
A）这可能吗？B) 你会怎么做？我见过基本的 scipy 拟合函数，但没有任何东西可以以这种方式处理 3D 数据（或者确实适合 2D 数组的类似函数）。
以前有人尝试过这样做可能是一个很长的机会，但任何帮助将不胜感激，
谢谢