現在、電場の計算と 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)
グリッド上の各ポイントでポテンシャル (V) の値を含む 3D numpy 配列を取得します。ただし、これでより便利なことを行うために、この値の配列を連続関数で近似できるようにして、グリッド上にない点でフィールドとその勾配を計算できるようにしたいと考えています。
A) これは可能ですか？B) これを行うにはどうすればよいですか? 基本的な scipy フィッティング関数を見てきましたが、この方法で 3D データを処理するものはありません (または実際に 2D 配列の類似関数に適合するものはありません)。
誰かがこれをやろうとしたことがあるかもしれませんが、助けていただければ幸いです。
ありがとう