StackOverflow2013

Note that there are some explanatory texts on larger screens.

plurals

POSuccessive over-relaxation not converging (when not done in-place)
text
Body
copied!<p>I'm trying to find the potential given some boundary conditions using the successive over-relaxation method.</p> <p>I have 2 solutions:</p> <p>-One iterates over all elements and applies the formula <code>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</code> in place. This is slow because it doesn't access memory in a nice way.</p> <p>-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.</p> <p>With alpha = 1.9 the first algorithm converges while the second one doesn't. For alpha = 1.0 both converge but very slowly.</p> <p>Can anyone tell me what I'm doing wrong? And how can I fix the fast solution.</p> <p>Full code:</p> <pre><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 </code></pre>

Querying!

Guidance

An individual column

Larger individual text columns get their own page to allow for proper reading.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload