I'm trying to find the potential given some boundary conditions using the successive over-relaxation method.
I have 2 solutions:
-One iterates over all elements and applies the formula 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 in place. This is slow because it doesn't access memory in a nice way.
-The other one, I create 4 matrices shifted in the 4 directions by 1. I Apply the same formula by then adding the matrices up. This however doesn't take into account modifications done during the current iteration. This is significantly faster then the previous one.
With alpha = 1.9 the first algorithm converges while the second one doesn't. For alpha = 1.0 both converge but very slowly.
Can anyone tell me what I'm doing wrong? And how can I fix the fast solution.
Full code:
#! 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
