This is a variant of mg5ptdemo. Here we compute the solution via multigrid and print the error at each iteration. We then use multigrid as a preconditioner for conjugate gradients and print the error for that.
import numpy as np
def gridfn(n):
# create a zero grid function with shape (n+1, n+1)
return np.zeros((n+1, n+1))
def mg5pt(f, u):
n = f.shape[0] - 1 # deduce n from the shape of f
if n == 2:
# exact solve when there is only 1 unknown
v = gridfn(2)
v[1, 1] = f[1, 1]/16 # h = 1/2, so 4/h^2 = 16
else:
v = u.copy()
h = 1./n
n2 = n/2
# presmooth with Gauss-Seidel
for i in range(1, n):
for j in range(1, n):
v[i, j] = (v[i+1, j] + v[i-1, j] + v[i, j+1] + v[i, j-1] + h**2 * f[i, j])/4.
# compute residual
r = gridfn(n)
r[1:n, 1:n] = f[1:n, 1:n] - (4.*v[1:n, 1:n] - v[2:n+1, 1:n] - v[0:n-1, 1:n] - v[1:n, 2:n+1] - v[1:n, 0:n-1])/h**2
# transfer to next coarser mesh
r2 = gridfn(n2)
r2[1:n2, 1:n2] = ( 4.*r[2:n-1:2, 2:n-1:2] \
+ 2.*(r[1:n-2:2, 2:n-1:2] + r[3:n:2, 2:n-1:2] + r[2:n-1:2, 1:n-2:2] + r[2:n-1:2, 3:n:2]) \
+ r[1:n-2:2, 1:n-2:2] + r[3:n:2, 1:n-2:2] + r[1:n-2:2, 3:n:2] + r[3:n:2, 3:n:2] ) / 16.
# recursively apply multigrid iteration on coarse mesh with initial guess zero to estimate error
e2 = mg5pt(r2, gridfn(n2))
# transfer error to fine mesh
e = gridfn(n)
e[2:n-1:2, 2:n-1:2] = e2[1:n2, 1:n2]
e[1:n:2, 2:n-1:2] = (e2[0:n2, 1:n2] + e2[1:n2+1, 1:n2])/2.
e[2:n-1:2, 1:n:2] = (e2[1:n2, 0:n2] + e2[1:n2, 1:n2+1])/2.
e[1:n:2, 1:n:2] = (e2[0:n2, 0:n2] + e2[1:n2+1, 0:n2] + e2[0:n2, 1:n2+1] + e2[1:n2+1, 1:n2+1])/4.
# correct
v = v + e
# postsmooth with Gauss-Seidel in reverse order
for i in range(n-1, 0, -1):
for j in range(n-1, 0, -1):
v[i, j] = (v[i+1, j] + v[i-1, j] + v[i, j+1] + v[i, j-1] + h**2 * f[i, j])/4.
return v
Next we try it out. We define \(n\) (a power of 2), set up the grid, and define the rhs grid function f.
# set up the grid
n = 2**5
h = 1./n
x = np.linspace(0, 1, n+1)
y = np.linspace(0, 1, n+1)
[X, Y] = np.meshgrid(x, y, indexing='ij') # (X(i,j), Y(i,j)) = (i h, j h)
# specify the function f(x,y)
def ff(x, y):
return -190*x**2 * (1-y) + 26.4*np.exp((x-.25)**2 + (y-.25)**2)
# evaluate it to define the grid function f
f = ff(X, Y)
For comparison purposes only, we form the matrix and compute the solution with a direct solve.
from scipy.sparse import *
from scipy.sparse.linalg import spsolve
I = eye(n-1) # (n-1) x (n-1) identity matrix
e = np.ones(n-1) # vector (1, 1, ..., 1) of length n-1
e0 = np.ones(n-2) # vector (1, 1, ..., 1) of length n-2
A1 = diags([e0, -4*e, e0], [-1, 0, 1]) # (n-1) x (n-1) tridiagonal matrix
A2 = diags([e0, e0], [-1, 1]) # same with zeros on diagonal
A = -kronsum(A1, A2, format='csr')
b = h**2 * np.reshape(f[1:n, 1:n], (n-1)**2) # right hand side
U = spsolve(A, b)
uex = gridfn(n)
uex[1:n, 1:n] = np.reshape(U, (n-1, n-1))
Now compute the multigrid iterates and print the error.
# initial guess, here zero
u = gridfn(n) # zero grid function
print 'iteration 0. max norm error:{:12.8f}'.format(abs(u-uex).max())
niters = 10
for niter in range(niters):
u = mg5pt(f, u)
print 'iteration {:2d}. max norm error:{:12.8f}'.format(niter+1, abs(u-uex).max())
Next we iterate preconditioned conjugate gradient, using mg5pt as the preconditioner.
# initial guess, here zero
u = gridfn(n) # zero grid function
print 'iteration 0. max norm error:{:12.8f}'.format(abs(u-uex).max())
# compute residual and initial search direction
r = f
s = mg5pt(r, u)
r2 = sum(sum(r*s)) # corresponds to inner product of vectors
niters = 10
for niter in range(niters):
# compute t = A*s
t = gridfn(n)
t[1:n, 1:n] = (4*s[1:n, 1:n] - s[0:n-1, 1:n] - s[2:n+1, 1:n] - s[1:n, 0:n-1] - s[1:n, 2:n+1])/h**2
s2 = sum(sum(s*t))
lmbda = r2/s2
u = u + lmbda*s
print 'iteration {:2d}. max norm error:{:12.8f}'.format(niter+1, abs(u-uex).max())
r = r - lmbda*t
r2old = r2
z = mg5pt(r, gridfn(n))
r2 = sum(sum(r*z))
s = z + (r2/r2old)*s