Vedant Puri
PhD student at Carnegie Mellon University
Mechanical Engineering, Carnegie Mellon University
Agenda
- Motivation: PDE Solvers
- Polynomial Interpolation
- Orthogonal Polynomials
- Properties
- Generating polynomials
PDEs solvers are needed for predictive modeling
\partial_t \vec{v} + (\vec{v}\cdot\nabla)\vec{v} = -\nabla p + \frac{1}{\text{Re}}\Delta \vec{v} + f\\
\nabla\cdot\vec{v} = 0
Navier-Stokes Equations
(Flow past bluff body \( Re = 3900 \))
Need high quality function representation over (complex) geometry.
Main operations: \(\nabla, \, \int_\Omega\)
High-order interpolation is the underlying technology
With a basis over \(\Omega\),
\{
\phi_0(x), \, \phi_1(x), \cdots \phi_{N-1}(x)
\}
\frac{\mathrm{d}}{\mathrm{d} x} f(x) = \Sigma_{n=0}^{N-1} f_i \frac{\mathrm{d}}{\mathrm{d} x} p_n(x)
\int_\Omega f(x) \mathrm{d} x = \Sigma_{n=0}^{N-1} f_n \int_\Omega p_n(x) \mathrm{d} x
Functions can be written as expansions
f(x) = \Sigma_{n=0}^{N-1} f_i p_n(x)
\(\nabla, \, \int\) can be done on the basis
(Tensor-product basis)
(Linear Finite Element basis functions)
Exact
Fast
Properties of the basis are endowed on to the function
Well-conditioned
Lagrange interpolants and Gauss Quadrature are all we need
For a set of nodes
\{
x_0, \, x_1, \, \cdots, x_{N-1}
\}
Lagrange Polynomials can be defined
\ell_i(x) = \frac{\prod_{j\neq i}x - x_j}{\prod_{j\neq i} x_i - x_j}
\ell_i(x_j) = \delta_{ij}
(Lagrange Polynomials over Gauss-Lobatto-Legendre nodes)
Differentiation
Interpolation
What about integration?
For a special set of nodes
\int_\Omega w(x)f(x) \mathrm{d} x = \sum_{n=0}^{N-1} w_n f(x_n)
Integration
Gauss Quadrature
\(\{x_n\}_n\) are roots of the \(n^\text{th}\) orthogonal polynomial generated by \(w(x)\) in \(\Omega\).
Conditioning
Orthogonal Polynomials
\{
1, x, x^2, x^3, \cdots, x^{N-1}
\}
\langle p_i, p_j \rangle
=
\int_a^b p_i(x) p_j(x) \underbrace{w(x) \mathrm{d}x}_{\mathrm{d}\alpha}
Orthogonalization
(Gram-Schmidt)
Start with the monomials
Orthogonal WRT inner product on \( C(\mathbb{R}) \)
\{
p_0(x), p_1(x), \cdots, p_{N-1}(x)
\}
\(w(x)\) generates a unique family of orthogonal polynomials (up to normalization)
(Classical Orthogonal Polynomials)
Sturm-Louisville Theory
\frac{\mathrm{d}}{\mathrm{d} x}\left( p(x) \frac{\mathrm{d} }{\mathrm{d} x}\phi(x) \right)
+ q(x) \phi(x) = -\lambda w(x)\phi(x),\\
x\in\Omega, + \,\text{BC}
Classical Orthogonal Polynomials emerge as solutions to the Sturm-Louisville problem
Example: Fourier (trigonometric) polynomials are the eigenfunctions of the Laplace equation
-\frac{\mathrm{d}^2}{\mathrm{d}x^2}\phi(x) = \lambda \phi(x), \,
x\in(-1,1)\\
u|_{-1} = u|_{1}
\phi_{n}(x) = \sin(n\pi x) + \cos(n\pi x)\\
\lambda_n = n^2
But boundary condition is embedded in the interpolants. So they are not general purpose.
Trigonometric polynomials have fast transforms: \( \mathcal{O}(N\log(N)) \)
GOAL: Finding orthogonal polynomials is arbitrary geometries
Want to do calculus in arbitrary geometries
Need orthogonal basis over \(\Omega\):
- Gauss quadrature nodes and weights
- Interpolation, differentiation operators
- Matrices should be sparse, well-conditioned
- Method to impose boundary condition
Aim
- Develop method to generate orthogonal polynomials in 1 dimension
- Extend to multiple dimensions
Problems
- High-order polynomials (except in special cases) are nasty for numerical computation
- No known recurrence relation in 2D/3D
The three-term recurrence relation
(Favad, 1935) All families of orthogonal polynomials follow a 3-term recurrence relation
T_0(x) = 1\\
T_1(x) = x\\
T_{n+1}(x) = 2x T_n(x) - T_{n-1}(x)
Example: Chevyshev polynomials (first kind)
x p_n(x) = c_n p_{n-1}(x) + a_n p_n(x) + b_n p_{n+1}(x)
In general, for \( \{p_n\}_{n=0}^{N-1} \)
In matrix form, we have a tridiagonal system
(Golub & Welsch, 1968) eigenvalues of \(J\) correspond to zeros of \(p_n(x)\).
(Mysovskih, Gautschi) \(J\) is uniquely determined by \(w(x)\)
(Mysovskih, Gautschi) \(J\) is uniquely determined by \(w(x)\)
The moments of the \( w(x) \)
\mu_k = \int_\Omega w(x) x^k \mathrm{d}x
(Golub & Welsch) eigenvalues of \(J\) correspond to zeros of \(p_n(x)\)
Graphs Corresponding to Orthogonal Polynomials
\(p_n(x)\) is the Characteristic poly
Graphs Corresponding to Orthogonal Polynomials
Method
w(x)
Recurrence relation
Gauss Quadrature
Calculus
Conclusions
Implementation Problems
- Implementation
- High-order polynomials are notoriously ill-conditioned
- Gram-Schmidt failing for N > 20
- There are methods to get past numerical instability
- No guarantee on conditioning of interpolation, differentiation matrix
-
Mysovskih, Gautschi - algorithm requires moments computation
Chicken and egg problem - to get Gauss quadrature for integration, you have to do integration
\mu_k = \int_\Omega w(x) x^k \mathrm{d}x
Thanks
orthogonal_polynomials
By Vedant Puri
orthogonal_polynomials
Orthogonal Polynomials
