# [Numpy-discussion] need help with simple conjugate gradient Laplace solver

rob rob at pythonemproject.com
Sat Mar 30 09:00:03 CST 2002

```After I post, I always see the dumb error.  I am not including the 6x
term in my finite difference equation.  It now converges, but I get
wierd looking V map.  Rob.

Here is the fixed code

from math import *
from Numeric import *

#
#***  ENTER DATA
filename= "out"
#
bobfile=open(filename+".bob","w")

print "\n"*30

NX=30    # number of cells

NY=30

NZ=30

N=30   # size of box

#allocate arrays

##Ex=zeros((NX+2,NY+2,NZ+2),Float)
##Ey=zeros((NX+2,NY+2,NZ+2),Float)
##Ez=zeros((NX+2,NY+2,NZ+2),Float)

p=zeros((NX+1,NY+1,NZ+1),Float)
q=zeros((NX+1,NY+1,NZ+1),Float)
r=zeros((NX+1,NY+1,NZ+1),Float)
u=zeros((NX+1,NY+1,NZ+1),Float)

u[0:N,0:N,0]=0    # box at 1V with one side at 0V
u[0:N,0,0:N]=1
u[0,0:N,0:N]=1
u[0:N,0:N,N]=1
u[0:N,N,0:N]=1
u[N,0:N,0:N]=1

r[1:NX-1,1:NY-1,1:NZ-1]=(u[1:NX-1,0:NY-2,1:NZ-1]+            #initialize
r matrix
u[1:NX-1,2:NY,1:NZ-1]+
u[0:NX-2,1:NY-1,1:NZ-1]+
u[2:NX,1:NY-1,1:NZ-1]+
u[1:NX-1,1:NY-1,0:NZ-2]+
u[1:NX-1,1:NY-1,2:NZ]-
6*u[1:NX-1,1:NY-1,1:NZ-1])

p[...]=r[...]   #initialize p matrix

#
#**** START ITERATIONS
#
N=(NX-2)*(NY-2)*(NZ-2)   # left over from Jacobi solution, not used
KK=0     # iteration counter

res=0.0;    #set  residuals=0

while(1):

q[1:NX-1,1:NY-1,1:NZ-1]=(6*p[1:NX-1,1:NY-1,1:NZ-1]-

p[1:NX-1,0:NY-2,1:NZ-1]-       # finite difference eq

p[1:NX-1,2:NY,1:NZ-1]-

p[0:NX-2,1:NY-1,1:NZ-1]-

p[2:NX,1:NY-1,1:NZ-1]-

p[1:NX-1,1:NY-1,0:NZ-2]-

p[1:NX-1,1:NY-1,2:NZ])

#    Calculate r dot p and p dot q

rdotp = 0.0
pdotq = 0.0

p[1:NX-1,1:NY-1,1:NZ-1]))
q[1:NX-1,1:NY-1,1:NZ-1]))

#    Set alpha value

alpha = rdotp/pdotq

#   Update solution and residual

u[1:NX-1,1:NY-1,1:NZ-1] += alpha*p[1:NX-1,1:NY-1,1:NZ-1]
r[1:NX-1,1:NY-1,1:NZ-1] +=  - alpha*q[1:NX-1,1:NY-1,1:NZ-1]

#    calculate beta

rdotq = 0.0

rdotq =

beta = rdotq/pdotq

#   Set the new search direction

p[1:NX-1,1:NY-1,1:NZ-1] = r[1:NX-1,1:NY-1,1:NZ-1] -
beta*p[1:NX-1,1:NY-1,1:NZ-1]

res = sort(ravel(r[1:NX-1,1:NY-1,1:NZ-1]))[-1]    #find largest
residual
#    resmax = max(resmax,abs(res))

KK+=1
#
print "Iteration Number %d  Residual %1.2e" %(KK,abs(res))
if (abs(res)<=resmax): break     # if residual is small enough break
out

print "Number of Iterations ",KK

```