This is the home page for the 18.336 course at MIT in Spring 2006, where the syllabus, lecture materials, problem sets, and other miscellanea are posted.

**Lectures**: Tuesday/Thursday 1–2:30pm (2-151).
**Office Hours:** Thursday 4:30–5:30 (2-388).

**Topics**: Advanced introduction to applications and theory of
numerical methods for solution of partial differential equations,
especially of physically-arising partial differential equations, with
emphasis on the fundamental ideas underlying various
methods. Discretization methods, including finite difference &
finite-volume schemes, spectral collocation, and Galerkin methods.
Uniform and non-uniform grids, finite-element and boundary-element
methods. Stability, consistency, and convergence analyses. Fast
solvers, including fast Fourier transforms, preconditioned conjugate
gradient, multigrid, and fast multipole methods. Adjoint methods and
sensitivity analysis. The course assumes familiarity with basic
(numerical) linear algebra and will involve some programming in
Matlab.

**Online notes**: at least for some portions, the course will closely
follow the syllabus from previous terms. See: 18.336
on OpenCourseWare.

**Grading**: 33% problem sets (about six). 33% mid-term exam
(Thurs. April 6), 34% final project (last month).

Introduced a few simple classes of PDEs: elliptic (Poisson, Laplace, Helmholtz), parabolic (diffusion), hyperbolic (waves). Broad overview of "discretization" of unknown function (finite difference, finite elements, spectral methods, boundary elements), "discretization" of PDE (finite differences, Galerkin, integral equations), imposition of boundary conditions, and solution of resulting linear (or nonlinear system) by iterative methods. Briefly discussed questions of convergence, stability, and efficiency.

**Handouts:** Problem Set 1 (due 23 Feb
2006), introduction to FFTs (from an IAP
seminar; you need not do the "homework" problem)

Introduced spectral methods, a.k.a. **spectral collocation**,
a.k.a. pseudospectral methods. See also: J. P. Boyd, *Chebyshev
and Fourier Spectral Methods* (Dover, 2000) online.
Some useful Matlab code is available at this page.

Considered Poisson's equation with periodic boundary conditions, starting with Fourier series solution. Discussed truncation of Fourier series by constrainng solution at collocation points. Trigonometric interpolation polynomials and the discrete Fourier transform (DFT). Aliasing. Solving PDE via DFT and IDFT (inverse DFT). Introduced the concept of fast Fourier transform algorithms (FFTs) and discussed their variety and significance.

**Note**: due date of pset 1 was changed (above) due to holiday.

**FFT algorithms**: Cooley-Tukey algorithm and one- to
multi-dimensional mapping of DFT. Implementation choices and
real-world performance (cache oblivious algorithms, factorization
choice, base cases, bit reversal); see also FFTW. Other FFT algorithms:
Prime-Factor Algorithm (a.k.a. Good-Thomas); Bluestein's and Rader's
algorithms for prime sizes (and large prime factors).

Convergence of Fourier series and DFT. Exponential vs. algebraic
convergence and Darboux's principle (singularities determine
convergence rate). Proved upper bound of O(1/|k|^{m})
convergence of series coefficients when function is differentiable m-1
times. Inferred bound of O(1/N^{m-1}) on Chebyshev norm
(maximum absolute error) for series with N terms.

Convergence of Fourier series and DFT continued. Considered
L_{2} (root-mean-square) error of truncated Fourier series and
showed that it is O(1/N^{m-1/2}) (m defined as above) for
algebraic convergence. Considered "discretization" error of DFT
vs. Fourier series and showed that it is equivalent to aliasing and
goes as O(1/N^{m}) in coefficients and contributes
O(1/N^{m-1/2}) to L_{2} error.

Considered convergence of DFT-spectral solution to Poisson's
equation, and showed that it is dominated by O(1/N^{m})
discretization error. Thus, for the homework where m=1 one might
expect O(1/N) convergence but instead you should get
O(1/N^{2})!! Explained how this arises in this special case
because we forced the k=0 term in the DFT to have zero error.

Addition to problem 2 in pset1: Since I explained 2(b) in class, and since you have two extra days, I would like you to repeat 2(b) for: (c) a ρ(x) that is continuous with discontinuous slope; (d) a ρ with discontinuous second derivative only; and (e) a ρ showing exponential convergence. That is, make up three new ρ functions (don't forget that the integral of ρ must be zero!) and plot the L_{2}error in φ vs. N on a log-log scale to find the convergence rate.

Discussed application of spectral methods when the operator isn't
diagonal as it is for Poisson, in particular for the example of the
Schrodinger equation. You do spatial operations on f(x) and
derivatives on c_{k}, converting back-and-forth via FFTs, and
thus can apply iterative methods in O(N log N) time.

Considered alternative basis of **Chebyshev polynomials** (of the first
kind) on finite intervals [-1,+1]. Showed relation to Fourier
(cosine) series via x=cosθ change of variables, and how they can
be computed quickly via an FFT (type-I
DCT). Discussed differentiation via FFT.

Brief overview of spectral elements and domain decomposition, and coordinate transformations to map to infinite domains (giving rise to rational Chebyshev functions).

**Handouts:** Pset 1 solutions, Problem Set 2 (due Thurs. Mar. 9).

**Announcement:** mid-term exam is scheduled for Thursday, April 6.

Brief overview of spectral methods in more than one dimension. O(N log N) complexity of multidimensional FFT, curvilinear spectral elements.

Introduction to **finite-difference methods**. See also 18.336
lecture notes on OpenCourseWare, chapter 1. As initial example,
focusing on one-way wave equation u_{x} + a u_{t} = 0.

Forward/backward/central differences (with linear/linear/quadratic accuracy). Defined/discussed consistency, stability, convergence, and well-posedness. Lax's (or Lax-Richtmeyer) equivalence theorem equating stability and convergence (sketched short proof, at least for smooth functions).

Introduced DTFT (discrete-time Fourier transform, although we apply it to the spatial dimensions). Introduced Von-Neumann stability analysis.

**Von-Neumann analysis** for several schemes:

- forward-time, backward-space: |g|
^{2}= 1 - 4aλ (1-aλ) sin^{2}(θ/2). Stable for aλ between 0 and 1, dissipative. - forward-time, center-space: unstable unless Δt ~ Δx
^{2} - Lax-Wendroff: stable for |aλ| less than 1, dissipative
- fully-implicit time-stepping: unconditionally stable, dissipative
- Crank-Nicolson: unconditionally stable, |g| = 1.

Discussed connection of Von-Neumann analysis to 18.03 stability analysis of coupled first-order equations: we are just finding the eigenvalues of the propagation matrix, which is easy for linear homogeneous problems (eigenvectors are Fourier modes, due to representation theory).

Proved necessary and sufficient condition for stability: |g| bounded above by 1+KΔt for some K.

Discussed solution of implicit schemes like Crank-Nicolson: requires solving sparse linear equations at every time step: either use iterative method, or exploit fact that matrix is tridiagonal in 1d (or product of tridiagonal, for higher-dimensional ADI = alternating-difference implicit schemes).

Crank-Nicolson was second-order in time but required implicit time-step. Introduced Leap-frog method (second-order in time also). Here, propagating 2 components so Von-Neumann analysis leads to eigenvalue analysis of 2×2 G matrix. Briefly discussed degenerate-eigenvalue (defective) case. Stability condition is then |aλ| < 1 (CFL condition).

Leap-frog for "real" 2-way wave equation u_{tt}=a^{2}u_{xx}. Proved stability under CFL condition, mentioned extension to higher dimensions.

Explained why dissipation can be a feature, not a bug: dissipative schemes attenuate high-frequency more than low-frequency, sifting the gold from the dross. Note that not all schemes with numerical dissipation are "dissipative" in this sense (example: Lax-Friedrichs).

Animation demo: animwave.m Matlab file for
one-way wave equation u_{t}+au_{x}=0 with periodic
boundary conditions. Showed numerical dissipation killing
high-frequency noise for square-wave pulse, but dissipating good stuff
too. Showed non-dissipative Crank-Nicolson giving crazy
high-frequency stuff for square-wave pulse, but reasonable behavior
for smooth pulse. Showed group vs. phase velocity and dispersion
(pulse spreading).

Defined **group and phase velocity**, and group-velocity dispersion (=
ω-dependent group-velocity). Discussed origin of dispersion in
physical systems: waveguide dispersion from inhomogeneities, material
dispersion from frequency-dependent (non-instantaneous) response, and
dispersion from discretization.

Derived group velocity by considering a wave packet (narrowly
peaked in Fourier domain) and Taylor-expanding ω(β) around
peak β_{0}. Discussed effects of dispersion and
dispersion parameter d^{2}β/dω^{2} on pulse
spreading, etc. Noted problems with group-velocity concept in
dissipative systems.

Analyzed phase and group velocity in Crank-Nicolson scheme, and showed why we observed negative group velocities for high-frequency components in the demo of the previous lecture.

Introduced the concept of absorbing boundary regions to simulate
infinite domains with a finite computational cell, and in particular
the concept of **Perfectly Matched Layers** (PML) developed by
Berenger (1994).

Considered the "real" wave equation u_{tt} = a^{2}
u_{xx}, and broke it into two coupled equations u_{t}
= b v_{x} and v_{t} = c u_{x}, with bc =
a^{2}. Showed that condition for zero reflection at an
interface was that √(c/b), the "impedance" of the media, is
matched. Concluded that we can create a reflectionless absorbing
medium by using b/s and c/s for any s, and in particular showed that
we get ω-independent absorption by using s=1+iσ/ω.

**Handouts:** Pset 2 Solutions, Problem Set 3 (due Thurs. 23 March).

Continued analysis of PML boundary regions. Showed that s=1+iσ/ω corresponds to simply adding -σu and -σv dissipation terms to the two equations.

Discretization of PML wave equation: introduced the staggered-grid leap-frog method. Noted that if we increase σ discontinuously in the PML regions, we will get numerical reflections — instead, one increases it e.g. quadratically or cubically with x.

For ease of generalization, introduced "**stretched-coordinate**"
viewpoint on PML (Chew and Weedon, 1994), where s is viewed as a
simple re-mapping into complex coordinates and hence is automatically
reflectionless. Showed PML for 2d scalar wave equation as example.
Noted applicability to other coordinate systems, other wave equations,
other numerical methods (e.g. spectral or finite elements).

Introduced **parabolic equations** (chapter 2 of OCW notes): the
heat/diffusion equation u_{t} = b u_{xx}. Discussed
how they arise physically, conservation of u, and the fact that the
L_{2} norm of u is monotonically decreasing. Discussed
consequences of this for stable consistent finite-difference schemes:
we must have |g|<1 for small θ ≠ 0.

Proved decreasing-slope property of heat-equation solution. Discussed forward-time, center-space discretization and analyzed vaia Von-Neuman. Discussed "finite-volume" methods based on discretizing integral conservation laws, and re-derived forward-time center-space scheme.

Analyzed steady-state solution, and proved that heat-equation approaches solution of Poisson's equation. Discusssed solution of Poisson's equation by center-difference scheme, and showed that it is a quadratic approximation to the spectral-collocation method we studied long ago. Discuss higher-order schemes, and why exponential-accurate (spectral) scheme corresponds to infinite-range finite-difference schemes.

Discussed time-stepping and Von-Neumann analysis for (spatial) spectral basis. Analyzed 4th-order Runge-Kutta scheme via Von-Neumann. Analyzed 2nd-order Crank-Nicolson scheme for diffusion equation.

Du-Fort Frankel scheme: unconditionally stable, but conditionally consistent. Discussed convection-diffusion equation.

Introduced a new way to get discretized PDEs: **weighted-residual
and Galerkin methods**. (See chapter 3 of the J. P. Boyd book,
above, which is available online.) Defined residual of PDE, defined
inner product, and defined weighed-residual approach in terms of N
weight functions. Showed how this gives N×N linear equations
for a linear PDE.

Continued discussion of weighted-residual methods. Showed how spectral collocation methods appear as the special case of delta-function weights. Showed least-squares method (which always leads to Hermitian positive-definite matrices for linear PDEs).

**Galerkin method**: weight functions equal basis functions. (Example:
using Fourier basis gives usual Fourier-series approach.) Showed that
if PDE is Hermitian and/or positive-definite then Galerkin matrices
are too. Showed how Galerkin for eigenproblems gives generalized
eigenproblems.

Showed how Galerkin methods require weaker smoothness conditions on basis than, e.g. collocation methods, via integration by parts.

Discussed physical and mathematical equivalence of Galerkin methods (for Hermitian problems) to minimization of "energy", both for linear equations and for eigenproblems (via the variational theorem). Gave examples of Euler-Bernoulli beam equation and Schrodinger eigen-equation.

**Handouts:** Pset 3 solutions
(Matlab files: pset3prob2.m, pset3prob3.m, pset3prob3b.m)

**Finite element methods** (FEM). Described general outlines, and gave
1d example of linear (first-order) elements ("tent functions").
Showed close connection of Galerkin FEM to finite-difference methods
for uniform grid (where gives 2nd-order method) and non-uniform grid
(where gives 1st-order method), in example of Poisson's equation.

Briefly discussed higher dimensions, e.g. 2d triangular meshes with tent-function bases. Discussed evaluation of element integrals via numerical quadrature.

**Numerical quadrature** (or cubature, in higher dimensions). Defined
the problem of integrating a function f(x) over an interval [a,b] with
a given weight function w(x) by evaluating f at N points. Showed how
to integrate via Lagrange polynomial interpolation through N points.

Introduced Gaussian quadrature: showed that all polynomials of
degree < 2N can be integrated exactly if we use the roots of the
degree-(N-1) orthogonal polnomial over the given interval. (See also
Arfken & Weber, *Mathematical Methods for Physicists*.)
Different weight functions give different orthogonal polynomials (most
common are Legendre polynomials for w(x)=1).

Clenshaw-Curtis quadrature: quadrature by using Chebyshev
polynomials at Chebyshev points (*not* at roots), via DFT/DCT to
get expansion coefficients.

Discussed adaptive quadrature: using embedded lower-order rule to get error estimate, and subdivide interval if error is too large. (via Clenshaw-Curtis or Gauss-Kronrod quadrature.)

**Handouts:** Problem set 4 (due Thursday April 20). You may find the following files useful for problem 2: schrodinger_fd.m, schrodinger_fd_err.m, makegrid.m.

**Iterative methods for solving linear equations.**

Stationary schemes: Jacobi, Gauss-Seidel, SOR. Showed how convergence relates to spectra of operators, and explained why diagonal dominance is required for Jacobi/Gauss-Seidel. Related Jacobi to "method of relaxation" for Laplace/Poisson problem. Compared convergence rates.

Non-stationary schemes: conjugate-gradient, steepest-descent, GMRES, biconjugate gradient (and variants like biCGSTAB(L)). Introduced steepest-descent and line minimization. Introduced conjugate-gradient method, defined conjugacy condition, and proved convergence in N steps via property of giving the global minimum over the span of all previous search directions.

Derived preconditioned conjugate-gradient method, using an approximate inverse for the operator to reduce the condition number. Discussed convergence rates and the possibility of superlinear convergence. Nonlinear conjugate gradient and preconditioning compared to Newton's method.

Methods for finding **preconditioners** (in general a "black art" that
requires problem-specific experimentation). Jacobi (diagonal)
preconditioners, especially for spectral methods. Solving a
simplified PDE. Incomplete Cholesky factorization (and incomplete LU
for nonsymmetric problems).

**Handouts:** Problem set 4 solutions (Matlab files for problem 2: schrodinger_galerkin.m, schrodinger_galerkin_err.m, schrodinger_galerkin_sparse.m.)

Online resources for iterative methods: The PETSc toolkit; the book, *Templates
for the Solution of Linear Systems: Building Blocks for Iterative
Methods* by Barrett *et al.*; the book, *Templates for
the Solution of Algebraic Eigenvlue Problems* by Bai *et
al.*.

**Operator splitting** and alternating-difference implicit (ADI)
methods. For ADI methods, see the course notes, §3.5, and also:
J. W. Thomas, *Numerical Partial Differential Equations:
Finite-Difference Methods* (Springer: New York, 1995), §4.4.
A useful review article is: R. I. McLachlan and G. R. W. Quispel, "Splitting
methods," *Acta Numerica* **11**, 341–434 (2002).

**Handout:** pset 4 contest results

Operator splitting continued: the split-step Fourier method for the
nonlinear Schrodinger equation. See also: G. P. Agrawal, *Nonlinear
Fiber Optics*, 3rd ed. (Academic Press: San Diego, 2001).

Discussed contest results from problem set 4 (see handout). Winner was Peng Yi, with an error about 1/2 of the uniform grid. (All but three contestants did worse than the uniform grid!)

Overview of grid generation and adaptation techniques. Discussed
structured vs. unstructured, and hybrid grids. Analyzed truncation
error of a structured grid for an arbitrary (sufficiently smooth)
coordinate mapping: for a center-difference approximation, nonuniform
and non-orthogonal grids still have quadratic convergence but with a
penalty in the constant factor. Discussed grid refinement,
heuristically balancing uniformity & orthogonality (a priori), and
adaptation to a particular solution (a posteriori, via an error
estimate or error indicator). Grid generation via quadtree/octree
refinement, Delaunay triangulation, solving elliptic PDEs (for
structured grids), etcetera. *In summary: a whole mess of
heuristics with a thin veneer of theory.* **Useful reference:**
*Handbook of Grid Generation*, J. F. Thompson, B. K. Soni, and
N. P. Weatherill, eds. (CRC Press: Boca Raton, 1999).

**Handout:** Notes on Adjoint Methods (by SGJ)

Basic introduction to **adjoint methods**, for evaluating the sensitivity
of some objective function g(**x**,**p**) to changes in design
parameters **p**, where **x** solves some equations depending on
**p**. Adjoint methods allow us to solve for the gradient
dg/d**p**, given **x**, about as quickly as solving for **x** once.

Went through the notes above: covered adjoint methods for linear
equations, nonlinear equations, eigenproblems, and simple
initial-value problems. Discussed applications to sensitivity
analysis and optimization ("inverse design"), and showed simple
example: inverse design of Schrodinger potential V(x) to obtain a
given ground-state wavefunction ψ_{0}(x).

Matlab code: schrodinger_fd_adj.m, schrodinger_fd_opt.m. You will also need the nonlinear conjugate-gradient code by H. B. Nielsen: conj_grad.m, linesearch.m.

A useful reference for more sophisticated adjoint problems is:
Y. Cao, S. Li, L. R. Petzold, and R. Serban, "Adjoint
sensitivity analysis for differential-algebraic equations: The adjoint
DAE system and its numerical solution," *SIAM
J. Sci. Comput.* **24** (3) 1076-1089 (2003).

Adjoint methods continued. Erratum: at the end of the last
lecture, I assumed the B and B_{p} commuted...gave the
more general solution today by showing how to properly differentiate
the matrix exponential (the Notes above contain the correct solution).

Considered more general time-dependent DAE problems and time-dependent/time-integrated objective functions, based on the treatment in Cao et al. (2003) above. (Only treated the special case of index 0 and 1 DAE systems, whereas the paper is more general.)

**Multigrid methods**. [Reference: U. Trottenberg, C. Ooseterlee, and a. Schüller, *Multigrid* (Academic Press: London, 2001).]

Basic idea: remove high-frequency errors with fine grid, low-frequency errors by interpolating to/from coarser grids. Reviewed why fine-grid methods (without preconditioning) converge low-frequency errors slowly (e.g. steepest descent, Jacobi, ...). Review stationary schemes, define iteration matrix M and relate its spectral radius (max. |eig.|) to convergence rate.

Focus on 1d Poisson's equation, as example. Von-Neumann analysis for Jacobi and ω-Jacobi. Introduced restriction R to coarse grid, and prolongation P back to fine grid, and considered several choices (injection, linear interpolation, full weighting, half weighting). Derived iteration matrix for 2-level multigrid.

Continued with multigrid. Von-Neumann analysis (a.k.a. Local Fourier Analysis) of 2-level multigrid. Showed how aliasiing means that we need 2×2 matrix (in one dimension). Derived eigenvalues of 2-level multigrid and showed that they are <1 for all Fourier components θ ... conclude that convergence occurs in fixed number of steps, independent of Δx.

Analyzed recursive application of multigrid. First, solve coarse
problems *completely* by recursive multigrid, and showed that
this leads to O(N^{2}) or worse performance. Instead,
consider *approximate* solution by one or two steps of multigrid,
recursively, leading to V-cycles and W-cycles.

Brief discussion of multigrid in higher dimensions, for finite elements, as preconditioner, etcetera.

**Boundary-element methods**. [Reference: Marc Bonnet,
*Boundary Integral Equation Methods for Solids and Fluids*
(Wiley: Chichester, England, 1999). Also L. Gaul, M. Kogl, and
M. Wagner, *Boundary Element Methods for Engineers and
Scientists* (Springer, 2003).] Focusing on example of Poisson's
equation for electrostatic potential given a charge distribution ρ
and a set of perfectly conducting regions and empty space Ω
(boundaries dΩ).

Started with "intuitive" 8.02 picture of problem: Green's function
is potential of point charge, discontinuity of potential gradient at
conductor interface gives a surface-charge density. The total
potential is then the convolution of the Green's function with ρ
and the surface charges, *plus* the potential energy of the
"induced" surface charge when you bring in a test charge.

Broke boundary-element methods (BEMs) into four steps:

- First, find an integral equation for the potential u(y) at a point
y in the interior of Ω, in terms of unknowns on
*surfaces*only. (Uses Green's functions + reciprocity relations.) - Second, take the limit as y approaches the surface dΩ (requires care due to singularities).
- Third, discretize the problem: write unknowns on
*surfaces*in terms of a finite basis to form N×N linear equations (for linear PDEs). (e.g. a finite-element discretization of the*surfaces*.) - Fourth, solve the equations. Since the matrices are dense, we need new tricks here to get a fast matrix-vector product: fast multipole methods, etc.

Outlined some advantages and disadvantages of BEMs.

(1) Derived reciprocity relation (just integration by parts) for
Laplacian, and used this to get integral equation for u(y) in terms of
the Green's function g_{y} and its gradient h_{y}.
Related to previous "intuitive" picture surface and induced charges.

(2) Began taking y to dΩ limit (following Bonnet book): put y on
dΩ, but then deform the surface slightly (by ε) to put y
*outside* the new volume Ω_{ε}.
Furthermore, showed how to weaken the h_{y} singularity by
adding zero in a clever fashion.

Continued BEM treatment. Finished step (2): showed that by assuming Hölder continuity locally, we can take the limit as ε goes to zero and get back an integral over the original surface. (The singularities are integrable.)

Showed that for exterior problems, where Ω is unbounded and part of our dΩ surface goes to infinity, we just need to add an additional u(y) term.

(3) Showed how to discretize the surface by weighted-residual
methods: as usual, expand unknowns in some finite basis (N terms) and
then integrate equations against N weight functions on the surface,
giving N×N linear equations. Discussed triangular meshes, shape
functions, and simple constant or linear basis functions. Compared
collocation methods (weight = delta function) and Galerkin methods
(weight = basis function). Discussed how to deal with singular
element integrals via Duffy transformations. [See M. G. Duffy, *SIAM
J. Num. Anal.* **19** (6), 1260 (1982). Also D. J. Taylor, *IEEE
Trans. Ant. Prop.* **51** (7), 1630 (2003).]

(4) Discussed how to get fast (O(N log N)) matrix-vector multiplies
for iterative solution of the equations. Briefly compared fast
multipole methods (FMM) [Greengard, *J. Comp. Phys.* **73**,
325 (1987)], pre-corrected FFT methods [e.g. J. R. Phillips
and J. K. White, *IEEE Trans. Comp.-Aided Design* **16**,
1059 (1997)], and wavelet matrix compression [e.g. Tausch & White,
*SIAM J. Sci. Comput.* **24**, 1610 (2003)].

Sketched basic ideas of FMM for Poisson problem with Dirichlet
boundaries and collocation approach. Potentials of nearby charges are
summed exactly, whereas potentials of far-away charges are summed
approximately using a multipole expansion (see e.g. Jackson,
*Classical Electrodynamics* for multipole expansions). In
particular, space is divided into a hierarchy of boxes, with multipole
expansions of larger and larger boxes taken for regions further and
further away, leading to an O(N log N) method. Multipole expansions
are computed recursively, coarser boxes from finer ones, using simple
relationships between the multipole coefficients around different
centers.