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/Nm-1) on Chebyshev norm (maximum absolute error) for series with N terms.
Convergence of Fourier series and DFT continued. Considered L2 (root-mean-square) error of truncated Fourier series and showed that it is O(1/Nm-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/Nm) in coefficients and contributes O(1/Nm-1/2) to L2 error.
Considered convergence of DFT-spectral solution to Poisson's equation, and showed that it is dominated by O(1/Nm) discretization error. Thus, for the homework where m=1 one might expect O(1/N) convergence but instead you should get O(1/N2)!! 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 L2 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 ck, 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 ux + a ut = 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:
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 utt=a2uxx. 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 ut+aux=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 d2β/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 utt = a2 uxx, and broke it into two coupled equations ut = b vx and vt = c ux, with bc = a2. 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 ut = b uxx. Discussed how they arise physically, conservation of u, and the fact that the L2 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/dp, 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 Bp 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(N2) 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:
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 gy and its gradient hy. 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 hy 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.)