我正在尝试使用连续的过度松弛方法在给定一些边界条件的情况下找到潜力。
我有 2 个解决方案：
-One 遍历所有元素并将公式应用field[y,x] = (1-alpha)*field[y,x] + (field[max(y-1,0),x] + field[min(y+1,field.shape[0]-1),x] + field[y,max(x-1,0)] + field[y,min(x+1,field.shape[1]-1)]) * alpha/4到位。这很慢，因为它不能以一种好的方式访问内存。
-另一个，我创建了 4 个矩阵，在 4 个方向上移动了 1。我应用相同的公式，然后将矩阵相加。然而，这并没有考虑在当前迭代期间所做的修改。这比前一个要快得多。
使用 alpha = 1.9，第一个算法收敛，而第二个算法不收敛。对于 alpha = 1.0，两者都收敛但非常缓慢。
谁能告诉我我做错了什么？以及如何解决快速解决方案。
完整代码：
#! python3
import numpy
import math
import time
def solve_laplace(boundary, mask, file = None, alpha = 1.0, threshold = 0.0001):
"""
We are using the successive over-relaxation method. We iterate until our solution changes less than some threshold value.
Vm+1(x,y,...) = alpha*( ((Vm(x-1,y,...) + Vm(x+1,y,...) + Vm(x,y-1,...) + Vm(x,y+1,...) + ...)/(2*nr dimensions) ) + (1-alpha)*Vm(x,y,...)
"""
dim = boundary.ndim
threshold = 0.0001
field = numpy.zeros_like(boundary)
numpy.copyto(field, boundary, casting = "safe", where = mask)
last_diff = float("infinity")
for iter_nr in range(10000):#max number of iterations
prev = field.copy() #make a copy of the field at the start of the iteration (python always stores pointers unless you explicitly copy something)
for d in range(dim): #can be scaled to arbitrary dimensions, using 2D for testing
#these 2 blocks are hard to follow but they work, read the comments
front = prev[tuple(0 if i==d else slice(None) for i in range(dim))] #select front face of cube/whatever
front = front[tuple(numpy.newaxis if i==d else slice(None) for i in range(dim))] #prepare it for next step
front = numpy.concatenate((front,prev),d) #add it the previous iteration's result
front = front[tuple(slice(-1) if i==d else slice(None) for i in range(dim))] #remove the back side of the previous iteration's result
#we now have the volume shifted right by 1 pixel, x now corresponds to the x-1 term
back = prev[tuple(-1 if i==d else slice(None) for i in range(dim))] #select back face of cube/whatever
back = back[tuple(numpy.newaxis if i==d else slice(None) for i in range(dim))] #prepare it for next step
back = numpy.concatenate((prev,back),d) #add it the previous iteration's result
back = back[tuple(slice(1,None) if i==d else slice(None) for i in range(dim))] #remove the front side of the previous iteration's result
#we now have the volume shifted left by 1 pixel, x now corresponds to the x+1 term
field += (front + back) * alpha/(2*dim) #this part of the formula: alpha*( ((Vm(x-1,y,...) + Vm(x+1,y,...) + Vm(x,y-1,...) + Vm(x,y+1,...))/(2*nr dimensions)
#numpy.copyto(field, boundary, casting = "safe", where = mask)
field -= alpha*prev #this part of the formula: (1-alpha)*Vm(x,y,...)
#reset values at boundaries
numpy.copyto(field, boundary, casting = "safe", where = mask)
#check if the difference is less than threshold
average = math.sqrt(numpy.average(field**2)) #sqrt of average of squares, just so i get a positive number
diff = math.sqrt(numpy.average((field-prev)**2)) #standard deviation
if last_diff < diff/average:
print("Solution is diverging.")
break
if diff/average < threshold:
print("Found solution after", iter_nr,"iteratiorn.")
break
last_diff = diff/average
if file is not None:
numpy.save(file,field)
return field
def solve_laplace_slow_2D(boundary, mask, file = None, alpha = 1.9,threshold = 0.0001):
"""
We are using the successive over-relaxation method. We iterate until our solution changes less than some threshold value.
Vm+1(x,y,...) = alpha*( ((Vm(x-1,y,...) + Vm(x+1,y,...) + Vm(x,y-1,...) + Vm(x,y+1,...) + ...)/(2*nr dimensions) ) + (1-alpha)*Vm(x,y,...)
"""
assert boundary.ndim == 2
field = numpy.zeros_like(boundary)
numpy.copyto(field, boundary, casting = "safe", where = mask)
last_diff = float("infinity")
start_time = time.time()
for iter_nr in range(10000):#max number of iterations
prev = field.copy()
for y in range(field.shape[0]):
for x in range(field.shape[1]):
if not mask[y,x]:
field[y,x] = (1-alpha)*field[y,x] + (field[max(y-1,0),x] + field[min(y+1,field.shape[0]-1),x] + field[y,max(x-1,0)] + field[y,min(x+1,field.shape[1]-1)]) * alpha/4
#check if the difference is less than threshold
average = math.sqrt(numpy.average(field**2)) #sqrt of average of squares, just so i get a positive number
diff = math.sqrt(numpy.average((field-prev)**2)) #standard deviation
if last_diff < diff/average:
print("Solution is diverging.")
break
if diff/average < threshold:
print("Found solution after the", iter_nr,"iteratiorn.")
break
if time.time() - start_time > 3600:
print("Completed in an hour time at iteration:", iter_nr)
break
last_diff = diff/average
#print(time.time() - start_time, iter_nr, last_diff)
if file is not None:
numpy.save(file,field)
return field
def test():
boundary = numpy.zeros((51,51))
boundary[25,25] = 1
for i in range(51):
boundary[0,i] = -1
boundary[50,i] = -1
boundary[i,0] = -1
boundary[i,50] = -1
mask = (boundary != 0)
print("Trying fast method:")
solve_laplace(boundary,mask,alpha = 1.5) #diverges
print("Trying slow method:")
solve_laplace_slow_2D(boundary,mask,alpha = 1.5) #converges but is very slow