This example by Aron Ahmadia and David Ketcheson is licensed under a Creative Commons Attribution 4.0 International License. All code examples are also licensed under the MIT license.
from IPython.core.display import HTML
css_file = './example.css'
HTML(open(css_file, "r").read())
Multigrid is one of the great algorithms for numerically solving PDEs. It is also one of the few essentially optimal algorithms for solving a linear system of equations, since it computes the solution of an $N\times N$ system in ${\mathcal O}(N\log(N))$ -- or even just ${\mathcal O}(N)$ -- operations.
This notebook is meant to accompany a reading of Section 4.6 of Randall LeVeque's text on finite difference methods. Other good resources are cited there.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
from clawpack.visclaw.JSAnimation import IPython_display
Let's use Jacobi's method to solve the one-dimensional boundary value problem
\begin{align} u''(x) & = f(x) & 0<x<1 \end{align}with boundary conditions $u(0) = 1, u(1) = 3$ where
\begin{align} f(x) & = -20 + \frac{1}{2} \phi''(x) \cos(\phi(x)) - \frac{1}{2} (\phi'(x))^2\sin(\phi(x)) \\ \\ \phi(x) & = 20\pi x^3. \end{align}This problem appears on page 103 of the text. In the box below, we set up the grid, the right hand side, and the boundary values.
m=2**8-1 # Number of grid points
h=1./(m+1) # Mesh width
x=np.linspace(0,1,m+2); x=x[1:-1] # grid
alpha=1.; beta=3. # boundary values
phi = lambda x: 20.* np.pi * x**3
f = lambda x: -20 + 0.5*120*np.pi*x * np.cos(phi(x)) - \
0.5*(60*np.pi*x**2)**2 * np.sin(phi(x)) # RHS
plt.plot(x,f(x)); plt.xlabel('x'); plt.ylabel('f(x)')
A centered-difference method for this problem gives us the system of equations
\begin{align} \frac{U_{i-1} - 2 U_i + U_{i+1}}{h^2} & = f(x_i). \end{align}Solving for $U_i$ gives \begin{align} U_i & = \frac{1}{2} ( U_{i+1} + U_{i-1} ) - \frac{h^2}{2} f(x_i). \end{align}
Jacobi's method consists of using the equation above as a fixed-point iteration:
\begin{align} U_i^{[k+1]} & = \frac{1}{2} ( U_{i+1}^{[k]} + U_{i-1}^{[k]} ) - \frac{h^2}{2} f(x_i). \end{align}In the cell below, implement Jacobi's method for the BVP. Don't forget about the boundary conditions. Take 100 iterations with it and plot the result.
U=np.linspace(alpha,beta,m) # Just use a straight line as initial guess
UU=[U]
F=0.5*h**2*f(x) # Construct the RHS
F[0]-=alpha/2.; F[-1]-=beta/2. # boundary conditions
e=np.ones(m-1)
G=0.5*(np.diag(e,-1)+np.diag(e,1)) # Jacobi matrix (average neighbors)
for i in range(100):
U=np.dot(G,U) - F
UU.append(U)
If you can, try storing the solution after each Jacobi iteration in a list, and make an animated plot of the solution after each iteration.
xf = np.linspace(0,1,1000);
uf = 1.+12.*xf-10.*xf**2 + 0.5*np.sin(phi(xf)) # Exact solution on fine grid
u = 1.+12.*x-10.*x**2 + 0.5*np.sin(phi(x)) # Exact solution on computational grid
fig = plt.figure()
ax = plt.axes(xlim=(0, 1), ylim=(-3, 6))
ax.plot(xf,uf)
line1, = ax.plot([],[],'-')
line2, = ax.plot([],[],'-o')
plt.legend(['Exact','Jacobi','Error'],loc=2)
def fplot(i):
line1.set_data(x,UU[i])
line2.set_data(x,UU[i]-u)
ax.set_title('# of iterations: '+str(i))
return line1,
animation.FuncAnimation(fig, fplot, frames=range(len(UU)))
Notice how the high-frequency components of the error are damped very quickly, but the low-frequency error barely changes. This is because the matrix $G$ appearing in Jacobi's method is tridiagonal, so in a single iteration each solution point only 'feels' the values of its immediate neighbors. It takes very many iterations for low-frequency information to propagate. As we'll see, part of the power of multigrid lies in enabling information to propagate rapidly between distant points in the solution.
At each Jacobi iteration, the error in the $p$th Fourier mode
$$u^p = \sin(p\pi x)$$
is reduced by a factor of
Which modes of the error converge most slowly?
The under-relaxed Jacobi method is obtained by taking, at each iteration, a weighted average of the previous iterate and the new iterate:
\begin{align} U_i^{[k+1]} & = (1-\omega)U^{[k]}_i + \omega\left(\frac{1}{2} ( U_{i+1}^{[k]} + U_{i-1}^{[k]} ) - \frac{h^2}{2} f(x_i)\right). \end{align}Now at each iteration the $p$th Fourier mode is reduced by a factor of
$$|\hat{\gamma}_p| := |1-\omega + \omega\cos(p\pi h)| = 1 + \omega(\cos(p \pi h) - 1)|.$$The code below plots this range of values. Try changing $\omega$ to $\frac{2}{3}$ and notice how the second half of the eigenvalues (for $p\ge m/2$) are all shifted to have smaller absolute value.
M=30; H=1./(M+1)
omegas=np.linspace(0.,1.)
p=np.arange(1,M+1)
fig = plt.figure()
ax = plt.axes(xlim=(0, M), ylim=(0,1))
ax.set_xlabel('p', fontsize=20)
ax.set_ylabel('$\widehat{\gamma}_p$', fontsize=20)
line1, = ax.plot([],[],'o')
def fplot(i):
omega = omegas[i]
gamma = abs(1.+omega*(np.cos(p*np.pi*H)-1))
line1.set_data(p,gamma)
ax.set_title('$\omega = %02.2f $' % omega, fontsize=20)
return line1,
animation.FuncAnimation(fig, fplot, frames=range(50))
In the cell below, implement under-relaxed Jacobi and run 100 iterations for the BVP with $\omega=2/3$. Does it make the results any better? Why or why not?
omega = 2./3.
U=np.linspace(alpha,beta,m) # Just use a straight line as initial guess
UU=[U]
F=0.5*h**2*f(x)
F[0]-=alpha/2.; F[-1]-=beta/2. # Construct the RHS, including boundary conditions
e=np.ones(m-1)
G=0.5*(np.diag(e,-1)+np.diag(e,1))
niter=100
for i in range(niter):
U=(1.-omega)*U+omega*(np.dot(G,U)-F)
UU.append(U)
xf = np.linspace(0,1,1000);
uf = 1.+12.*xf-10.*xf**2 + 0.5*np.sin(phi(xf)) # Exact solution on fine grid
u = 1.+12.*x-10.*x**2 + 0.5*np.sin(phi(x)) # Exact solution on computational grid
fig = plt.figure()
ax = plt.axes(xlim=(0, 1), ylim=(-3, 6))
ax.plot(xf,uf)
line1, = ax.plot([],[],'-')
line2, = ax.plot([],[],'-o')
plt.legend(['Exact','Jacobi','Error'],loc=2)
def fplot(i):
line1.set_data(x,UU[i])
line2.set_data(x,UU[i]-u)
ax.set_title('# of iterations: '+str(i))
return line1,
animation.FuncAnimation(fig, fplot, frames=range(niter))
The code below plots a mode of a certain frequency, and then shows how it is represented on grids with decreasing numbers of points. Notice how a low frequency mode (i.e., one with many points per wavelength) on a fine grid becomes a high frequency mode (one with very few points per wavelength).
Eventually, on grids with very few points, the mode is aliased to a completely different mode.
from matplotlib import animation
import matplotlib.pyplot as plt
from clawpack.visclaw.JSAnimation import IPython_display
fig = plt.figure()
ax = plt.axes(xlim=(0, 1), ylim=(-1.2, 1.2))
p=5
xf=np.linspace(0,1,1000) # fine grid
line1, = ax.plot(xf, np.sin(p*np.pi*xf), '-r', lw=2)
line2, = ax.plot([],[],'o-',lw=2)
def fplot(m):
x=np.linspace(0,1,m+2); # grid
line2.set_data(x,np.sin(p*np.pi*x))
ax.set_title('m='+str(m))
return line2,
animation.FuncAnimation(fig, fplot, frames=(63,31,15,7,3,1))
The multigrid idea is to solve successively on different grids, so that every mode has a small damping factor $\hat{\gamma}$ on some grid.
In order to implement multigrid, we need a way to take a function with values on a fine grid and restrict it to its values on a coarser grid (restriction). We also need a way to take a function with values on a coarse grid and interpolate onto a fine grid (prolongation). The functions below do just this, using linear interpolation in the case of the prolongation operator. Make sure you understand what they are doing. You may find it helpful to look up slicing notation in the NumPy help online.
def coarsen(f):
return f[1::2] # This slicing just takes the odd-numbered points
def interpolate(f,alpha,beta):
m_coarse=len(f)
m_fine =2*m_coarse+1
f_interp = np.zeros(m_fine)
f_interp[1::2]=f #Set the values of the odd numbered points
f_interp[2:-1:2]=0.5*(f[:-1]+f[1:]) #Set the values of the (interior) even numbered points
f_interp[0]=0.5*(f_interp[1]+alpha) #Set the values of the endpoints
f_interp[-1]=0.5*(f_interp[-2]+beta)
return f_interp
Now we compute the residual $$r = f-AU$$ on the fine grid, coarsen it, and solve $$Ae=-r$$ for the approximate error on the coarse grid. Then interpolate and subtract the interpolated error from the fine grid solution.
#compute the residual
A=2./h**2*(G-np.eye(m))
F=f(x); F[0]-=alpha/h**2; F[-1]-=beta/h**2
r=F-np.dot(A,U)
print max(abs(r))
omega=2./3.
m2=(m-1)/2. # Number of points in coarse grid
h2=1./(m2+1) # coarse mesh width
x2=coarsen(x) # coarse grid points
r2=coarsen(r) # residual restricted to coarse grid
err2=np.zeros(m2) # initial guess for the error
e2=np.ones(m2-1); G2=0.5*(np.diag(e2,-1)+np.diag(e2,1))
for i in range(3): #Solve Ae=-r by Jacobi iteration
err2=(1.-omega)*err2+omega*(np.dot(G2,err2)+0.5*h2**2*r2)
err=interpolate(err2,0,0) # interpolate the error
print err[0]
U=U-err # and use it to correct our solution
plt.clf()
plt.plot(x,u,x,U,x,U-u)
plt.legend(['Exact','Jacobi','Error'],loc='best')
Finally, in order to keep our multigrid code clean, we'll write a function to set up and take a fixed number of under-relaxed Jacobi iterations. The function returns the approximate solution and the residual.
def Jacobi(U,f,alpha,beta,m,nu):
"""Perform nu Jacobi iterations on a grid with m points, with initial guess U, right hand side function f and
Dirichlet boundary conditions with values alpha and beta. The function returns both the approximate
solution and the residual."""
omega=2./3.
h=1./(m+1)
F=0.5*h**2*f.copy()
F[0]-=alpha/2.; F[-1]-=beta/2.
e=np.ones(m-1)
G=0.5*(np.diag(e,-1)+np.diag(e,1))
for i in range(nu):
U=(1.-omega)*U + omega*(np.dot(G,U)-F)
A=2./h**2*(G-np.eye(m))
FF=f.copy(); FF[0]-=alpha/h**2; FF[-1]-=beta/h**2
rr=FF-np.dot(A,U)
return U,rr
Now let's do a full V-cycle. Look carefully through the code below until you understand what it does.
k=8;
m=2**k-1
rdep=k-1 # Recursion depth; this is how many grids down we want to go
# rdep=k-1 gives a full V-cycle
nu=3 # Number of Jacobi iterations to take at each step
U=np.linspace(alpha,beta,m) # Initial guess
x=np.linspace(0,1,m+2); x=x[1:-1] # grid
phi = lambda x: 20.* np.pi * x**3
u = 1.+12.*x-10.*x**2 + 0.5*np.sin(phi(x))
f = lambda x: -20 + 0.5*120*np.pi*x * np.cos(phi(x)) - 0.5*(60*np.pi*x**2)**2 * np.sin(phi(x))
alpha=1.; beta=3.
F=f(x)
r=[None]*(rdep+1); error=[None]*(rdep+1) # This just initializes these lists to have the right length
U,rr=Jacobi(U,F,alpha,beta,m,nu) # Initial iteration on fine grid
for i in range(1,rdep+1): # Going down
m=(m-1)/2 # = 2**(k-i) - 1
r[i]=coarsen(rr) # residual restricted to next coarser grid
error[i],rr=Jacobi(np.zeros(m),-r[i],0.,0.,m,nu)
for i in range(1,rdep): # Coming up
m=2*m+1
err=error[rdep-i]-interpolate(error[rdep+1-i],0,0) # Interpolate and subtract the correction
error[-i-1],rr=Jacobi(err,-r[rdep-i],0.,0.,m,nu)
m=2*m+1
U=U-interpolate(error[1],0,0) # final solution correction
U,rr=Jacobi(U,F,alpha,beta,m,nu) #Final iterations on original grid
plt.clf()
plt.plot(x,u,x,U,'o-',x,U-u)
plt.legend(['Exact','Jacobi','Error'],loc='best')
Since we have saved the errors on each grid, we can plot them to see how they represent the accumulation of errors at different scales.
fig = plt.figure()
ax = plt.axes(xlim=(0, 1), ylim=(-5,5))
line1, = ax.plot([],[],'-ok',lw=2)
def fplot(i):
y = error[-i]
m = len(y)
x=np.linspace(0,1,m+2); x=x[1:-1]
line1.set_data(x,y)
return line1,
animation.FuncAnimation(fig, fplot, frames=range(1,len(error)))
Now you have the basic ideas of multigrid. Of course, for this 1-dimensional problem (where the system to be solved is tridiagonal) it does not provide a big advantage over more naive methods. The power of multigrid is that it can be used in multiple dimensions, where the algebraic system is much more challenging for other methods. Furthermore, multigrid can be applied to nonlinear and time-dependent problems.
Modify the V-cycle code above to answer the following questions. Try to explain your results.
(a) How does the accuracy change as we change the number of Jacobi iterations performed at each step?
(b) Is it better to use a finer grid, or more Jacobi iterations if we want to improve the solution accuracy?
(c) What happens if we don't perform any Jacobi iterations in the "up" part of the V-cycle?
(d) What happens if we don't recurse all the way down to the 1-point grid?
(e) What happens if we use the original Jacobi method, or some other value of $\omega$?
Use the V-cycle code above to implement the full multigrid algorithm discussed in the text. You will probably find it helpful to write a Vcycle() function first.
from IPython.display import display, HTML
from IPython.html import widgets
from IPython.html.widgets import interact, interactive
def Jacobi(U,f,alpha,beta,m,nu,omega=2./3):
"""Perform nu Jacobi iterations on a grid with m points, with initial guess U, right hand side function f and
Dirichlet boundary conditions with values alpha and beta. The function returns both the approximate
solution and the residual."""
h=1./(m+1)
F=0.5*h**2*f.copy()
F[0]-=alpha/2.; F[-1]-=beta/2.
e=np.ones(m-1)
G=0.5*(np.diag(e,-1)+np.diag(e,1))
for i in range(nu):
U=(1.-omega)*U + omega*(np.dot(G,U)-F)
A=2./h**2*(G-np.eye(m))
FF=f.copy(); FF[0]-=alpha/h**2; FF[-1]-=beta/h**2
rr=FF-np.dot(A,U)
return U,rr
def Vcycle(k=8,nu=3,omega=2./3,smooth_after_interp=True):
m=2**k-1
rdep=k-1 # Recursion depth; this is how many grids down we want to go
# rdep=k-1 gives a full V-cycle
alpha=1.; beta=3.
U=np.linspace(alpha,beta,m) # Initial guess
x=np.linspace(0,1,m+2); x=x[1:-1] # grid
u = 1.+12.*x-10.*x**2 + 0.5*np.sin(phi(x))
F=f(x)
r=[None]*(rdep+1); error=[None]*(rdep+1) # This just initializes these lists to have the right length
U,rr=Jacobi(U,F,alpha,beta,m,nu,omega) # Initial iteration on fine grid
for i in range(1,rdep+1): # Going down
m=(m-1)/2 # = 2**(k-i) - 1
r[i]=coarsen(rr) # residual restricted to next coarser grid
error[i],rr=Jacobi(np.zeros(m),-r[i],0.,0.,m,nu,omega)
for i in range(1,rdep): # Coming up
m=2*m+1
err=error[rdep-i]-interpolate(error[rdep+1-i],0,0) # Interpolate and subtract the correction
if smooth_after_interp:
error[-i-1],rr=Jacobi(err,-r[rdep-i],0.,0.,m,nu,omega)
else:
error[-i-1] = err
m=2*m+1
U=U-interpolate(error[1],0,0) # final solution correction
U,rr=Jacobi(U,F,alpha,beta,m,nu,omega) #Final iterations on original grid
fig, ax = plt.subplots(figsize=(8,6),subplot_kw={'axisbg':'#EEEEEE',
'axisbelow':True})
ax = plt.axes(xlim=(0, 1), ylim=(-5,6))
ax.plot(x,u,x,U,'o-',x,U-u)
plt.legend(['Exact','Jacobi','Error'],loc='best')
interact(Vcycle, k=(4,10,1), nu=(1,4,1), omega=(0,2,1./6))