MIT 18.336J/6.7340J, Fall 2022
Fast Methods for Partial Differential and Integral Equations
This course broadly covers modern numerical methods for solving large-scale partial differential and integral equations. The focus varies year to year and is often updated to include state-of-the-art techniques. This semester we will focus in particular on Fourier and modern polynomial spectral methods, the fast multipole method, boundary integral equations, and applications to fluid dynamics and electromagnetism.
Catalog description: A unified introduction to the theory and practice of modern, near linear-time, numerical methods for large-scale partial differential and integral equations. Topics include: preconditioned iterative methods; generalized Fast Fourier Transform and other butterfly-based methods; multiresolution approaches including multigrid algorithms, hierarchical low-rank matrix decompositions, and low and high frequency Fast Multipole Methods. Example applications include: aircraft design, cardiovascular system modeling, electronic structure computation, and tomographic imaging.
Syllabus
Lectures: Tuesday/Thursday 9:30-11:00 AM in room 2-147.
Office Hours: Check Canvas calendar.
Prerequisites: This course covers advanced techniques for discretizing and solving PDEs. Some familiarity with ordinary differential equations, partial differential equaitons, Fourier transforms, linear algebra, and basic numerical methods for PDEs is assumed. It is strongly recommended that you have taken a previous course on basic numerical methods, such as 2.096/6.336/16.910, 2.097/6.339/16.920, 18.085, or 6.337/18.335. Problem sets will involve extensive coding and are required to be completed in Python or Julia notebooks.
Textbooks & Other Reading: Recommended reading will be posted as the class progresses. There is no textbook for the course, but the following books may be useful:
- Strauss "Partial Differential Equations: An Introduction". An advanced undergrad introduction to PDEs.
- Boyd "Chebyshev and Fourier Spectral Methods". Very readable and available online.
- Martinsson "Fast Direct Solvers for Elliptic PDEs". Modern and concise.
- LeVeque "Finite difference methods for ordinary and partial differential equations".
Grading: 50% problem sets (approximately biweekly), 50% final project report and presentation. Unless previous arrangements are made, late problem sets will be accepted for 1 week after the initial due date with a 50% penalty.
Collaboration Policy: Make a strong effort to solve problems on your own before discussing with any classmates. You must write up your own code and solutions, and indicate your collaborators on your assignments.
Problem Sets
Problem Set 1: Fast Fourier Transforms
Problem Set 2: Fourier and Finite Difference Helmholtz Solvers
Problem Set 3: Chebyshev Collocation
Problem Set 4: Ultraspherical Method
Problem Set 5: Low-Rank Methods
Final Projects
The final project is a 10-12 page paper and a 15 minute presentation during the last week of classes. The final project is broad in scope, but must include the implementation of a fast algorithm in Python or Julia along with performance and error analyses. The project can take the form a "literature review project" discussing a published algorithm or a "research project" attempting to implement a new solver for a problem of your choice.
Lit. review projects
One possibility is to review and implement an algorithm that was mentioned briefly or not covered in the course. Such a project should follow one or several published research papers describing the algorithm, along with a new implementation. Possible topics and suggested papers include:
- Advanced optimizations for production-level FFT implementations.
- Hierarchical Poincare-Steklov schemes using Chebyshev polynomails.
- Rectangular collocation with aliasing analysis.
- Alternating direction implicit scheme for multidimensional Chebyshev Poisson solvers.
- Bivariate orthogonal polynomials on triangles.
- Non-uniform Fast Fourier Transforms.
- Sparse spectral elements.
Research projects
Another option is to use the methods covered in class to implement a fast solver for a research problem in your field. This should include:
- A discussion of the scientific problem and a brief derivation of the model PDE.
- Mention of the current commonly used methods in the field for the problem.
- A fast implementation of a new solver for the problem, or a related first-step.
- Discussion of the prospects and limitations of fast/high-order techniques for this problem in the future.
The goal should be producing a functional solver matching or improving on existing techniques in certain cases. However, it is understood that this may not turn out to be feasible (or even possible), and that's what makes it research! In that case, an implementation for a related toy model is expected, along with an analysis of the barriers to making a solver for the original problem.
Report & presentation format
Your report should be written in the style a SIAM article, using the SIAM article LaTeX templates. Your report should be between 10 and 12 pages. The SIAM layout is very spacious so this is not a lot of space. Make sure your presentation and notation is concise, but also has enough content and is not just padded out with images!
Along with your report, you must submit the code implementing your algorithms, preferably in the form of a Jupyter notebook. Make good use of headings, text, comments, and descriptive function/variables names. Good code is code that is easily understandable by others!
Your report and presentation should both include:
- Background information for your algorithm / physical application.
- A concise mathematical description of the algorithm you're using.
- Performance and error analysis of your implementation.
Lecture Material and Summaries
Lecture 1: Introduction to fast methods, PDEs, IEs
Summary
- Course policies
- Why fast algorithms? History of fast algorithms for the Fourier transform.
- Why PDEs? Models for physical systems. Classes of PDEs. Elliptic regularity theorem.
- Why integral equations? Better conditioning from using exact solution formulae.
- Continuous and discrete Fourier transforms.
Related Reading
Lecture 2: Fast Fourier transforms
Summary
- History of FFTs. Facts that make FFTs possible.
- Radix-2 Cooley-Tukey algorithm.
- Radix-3 and algorithms for prime N.
- Sine, cosine, and other Fourier-related transforms.
Related Reading
Lecture 3: PDE discretization and preconditioning
Summary
- Discretizing linear boundary value problems.
- Discretizing nonlinear / initial value problems.
- Direct solvers and their complexity.
- Iterative solvers and their convergence.
- Fourier preconditioning using FFTs.
- Circulant matrices: diagonalization using DFT.
- Toeplitz matrices: embedding into circulant form.
- Hankel matrices: converison to Toeplitz form.
Lecture 4: Finite differences and fast Poisson solvers in 1D
Summary
- Review of finite difference methods.
- Spectral preconditioning for fast Poisson solvers in 1D:
- Periodic BCs using FFTs.
- Dirichlet BCs using DSTs.
- Neumann BCs using DCTs.
- Gauge fixing.
Lecture 5: Fast finite difference solvers in multiple dimensions
Summary
- Sylvester equations, Bartels-Stewart algorithm.
- Kronecker products and block matrices.
- Fast linear algebra with structured block matrices.
- Fast Poisson solvers in multiple dimensions.
- Extensions to other PDEs, e.g. Helmholtz.
- Limitations and alternatives for non-constant coefficients.
Lecture 6: Domain decomposition methods
Summary
- Extensions to other domains.
- Schur complement / Poincare-Steklov method for domain decomposition.
- Connecting two 1D segments.
- Connecting two 2D boxes.
- Heirarchical Poincare-Steklov method for multiple connections.
- Distorted domains.
Related Reading
Lecture 7: Introduction to spectral methods
Summary
- Spectral representations of functions.
- Rates of convergence.
- Truncation and discretization error.
- Convergence of Fourier series.
Related Reading
- Boyd chapters 1 & 2.
Lecture 8: Fourier spectral methods I
Summary
- Weighted residual method for discretizing PDEs.
- Fourier solver for Poisson equation with periodic BCs.
- Comparison to finite difference solver. Eigenvalue corrections.
- Fourier spectral collocation.
- Sine and cosine solvers for Poisson equation with Dirichlet/Neumann BCs.
- Conditions on forcing for exponential convergence.
- Parity mixing.
Related Reading
Lecture 9: Fourier spectral methods II
Summary
- Extending to multiple dimensions via direct products.
- Extending to systems of equations. Maintaining bandedness.
- Pseudospectral method for evaluating nonlinearities.
- Aliasing errors, dealiasing rules for arbitrary-order nonlinearities.
- Examples: Dedalus code, Yeung group turbulence simulations.
Related Reading
Lecture 10: Polynomial interpolation
Summary
- Polynomial interpolation.
- Lagrange representation.
- Lebesgue constant of equispaced points. Runge phenomenon.
- Lebesgue constant of Chebyshev nodes.
- Overview of orthogonal polynomials.
- Chebyshev polynomials.
- Convergence of Chebyshev series. Singularities in elliptical coordinates.
- Gaussian quadrature.
- Deriving nodes and weights.
- Discretization error from polynomial interpolation on Gaussian quadrature nodes.
Related Reading
- Boyd chapter 4.
Lecture 11: Chebyshev collocation methods
Summary
- Convergence of Gegenbauer polynomials. Legendre polynomials and finite elements.
- Collocation matrices for differentiation.
- Collocation matrices for multiplication. Aliasing errors.
- Roots vs extrema grid.
- Boundary bordering. Generalized tau equivalence.
- Rectangular collocation. Generalized tau equivalence.
Related Reading
Lecture 12: Dense Chebyshev spectral methods
Summary
- Newton's method in function space for nonlinear ODEs.
- Chebyshev recursion formulae.
- T-to-T matrices for differentiation.
- T-to-T matrices for multiplication. Band-limited expansions.
- Classical tau method for boundary conditions.
Lecture 13: Sparse Chebyshev spectral methods
Summary
- Chebyshev-T derivatives are Chebyshev-U polynomials.
- T-to-U differentiation and conversion matrices.
- Differentiation and conversion between ultraspherical polynomials.
- Ultraspherical method for arbitrary-order linear ODEs. Generalized tau equivalence.
- Woodbury matrix identity for banded + multi-bordered matrices.
Related Reading
Lecture 14: Multidimensional sparse methods
Summary
- Block-banded formulation for systems of equations.
- Pseudospectral method for nonlinear IVPs.
- Polynomial spectral methods in multiple dimensions:
- Chebyshev-Chebyshev, alternating direction implicit scheme
- Fourier-Chebyshev for periodic layers, cylinders, etc.
- Example: incompressible Navier-Stokes without splitting.
- Examples using the Dedalus code.
Related Reading
Lecture 15: Sparse methods for curvilinear domains
Summary
- Spectral methods with coordinate singularities:
- Non-smoothness of vector/tensor components.
- Radial regularity of Fourier components in polar coordinates.
- Radial regularity of vector/tensor components in polar coordinates.
- Modified Jacobi bases for tensors in the disk.
- Spin-weighted spherical harmonics for tensors on the sphere.
- Examples using the Dedalus code.
Related Reading
- Vasil et al. "Tensor calculus in polar coordinates using Jacobi polynomials"
- Vasil et al. "Tensor calculus in spherical coordinates using Jacobi polynomials. Part-I: Mathematical analysis and derivations"
- Lecoanet et al. "Tensor calculus in spherical coordinates using Jacobi polynomials. Part-II: Implementation and examples"
Lecture 16: Introduction to low-rank methods
Summary
- Scale separation in gravity.
- Direct N-body using Newton's law for point masses.
- Poisson equation for gravitational potential using Gauss's law for continuous distributions.
- Reformulating continuous problem as integral equation using Green's function.
- Both are limited in between: many point masses or well-separated distributions.
- Singular value decomposition.
- Basic definition and derivation.
- Optimal low-rank approximations using truncated SVDs.
- Fast matrix-vector products using low-rank approximations.
- Low-rank approximations to functions. Fast inner products for integral equations.
Related Reading
Lecture 17: Approximating low-rank interactions
Summary
- Numerically approximating the SVD using Gaussian elimination for matrices and functions.
- Low-rank nature of gravitational interaction.
- Projection: merging sources to compress the operator domain.
- Interpolation: merging measurements to compress the operator range.
- Duality of projection and interpolation.
- Analytical low-rank expansions using Taylor series.
- 2D and 3D multipole expansions. Error estimates.
Lecture 18: Fast multipole methods and boundary integral equations
Summary
- Hierarchical decompositions with error control.
- Barnes-Hut and FMM in 2D. Comparison to direct methods.
- Integral formulations of PDEs for complicated geometries, exterior problems.
- Basic ideas from potential theory for Laplace equation:
- Fictitious charges.
- Single-layer potentials.
Related Reading
Lecture 19: Discretizing boundary integral equations
Summary
- Double-layer potentials.
- Global quadrature discretizations.
- Quadrature for weakly singular kernels.
- Panel-based discretizations.
Related Reading
Lecture 20: Advanced boundary integral equations
Summary
- Domains with corners. Diadic refinement.
- BIEs for other equations:
- Stokes flow. Stokeslets and Stresslets.
- Linear elasticity.
- Helmholtz. Spurious resonances and combined field forms.
- Time-harmonic Maxwell.
- Summary and comparison to PDE solvers.
Related Reading