Polynomial Approximations: Why monomials are bad

Uncategorized

Define the concept of an inner product on some domain [a,b] as:

(f,g) = \int_a^b \! f(x)g(x) \, \mathrm{d}x = 0

If we consider vectors in R^n space the inner product is essentially the dot product of the two vectors. Vectors that are orthogonal to each other have a dot product of 0. The concept of orthogonality can be abstracted and extended to vectors in any vector space, so that vectors that have a zero inner product are orthogonal.

Some arbitrary function u(x,t) on the domain I_k =[x_{k-^1\!/_2},x_{k+^1\!/_2}] may not permit an analytical representation and so most be approximated in some way. An easy choice is a polynomial approximation of order M. A natural choice for the polynomial basis is the monomials, \psi_m(x) = x^{m-1}. These are satisfactory for low orders, but we encounter a problem at higher orders.

If we normalize the monomial basis with respect to the l^2 norm \|f\|=\sqrt{(f,f)} we arrive at \bar{\psi}_m(x) = \sqrt{2m+1}x^m . If the inner product of two basis functions is 0, they are orthogonal and therefore linearly independent. If many of our basis functions end up linearly dependent it is difficult to reconstruct the solution approximation.

Consider the inner product of two adjacent basis functions in our normalized monomial basis:

(\bar{\psi}_{m-1},\bar{\psi}_m) = \int_0^1 \! \sqrt{2m-1}x^{m-1}\sqrt{2m+1}x^{m} \, \mathrm{d}x = \sqrt{1-\frac{1}{4m^2}}

It is easy to see that for even moderate values of m the inner product is close to 1, indicating near linear dependence between basis functions.

It is therefore important to choose a polynomial basis of (ideally) orthogonal functions so that the solution approximation is well conditioned. One good choice is the Legendre polynomials P_m which are orthogonal on the interval [-1,1] (these are defined as the solutions to the Sturm–Liouville equation for a weighting function w(x)=1 ). They can be defined using a recurrence formula or the explicit Rodrigue’s formula.

P_m(x) = \frac{(2m-1)xP_{m-1}(x) - (m-1)P_{m-2}(x)}{m}

P_m(x) = \frac{1}{2^m m!}\frac{\mathrm{d}^2}{\mathrm{d}x^m}\left[(x^2-1)^m\right]

Orthogonality means that the inner product of two Legendre polynomials is zero if they are not of the same order, this can be expressed using the Kronecker product as

(P_n,P_m) = \int_{-1}^1 \! P_n P_m \, \mathrm{d}x = \delta_{nm} \left(\frac{2}{2m+1}\right)

Note that we could use a similar orthonormal version of the basis such that every inner product is simply the Kronecker product, however one nice property of the current form is that P_m(1)=1 and P_m(-1)=(-1)^m . In the orthonormal case we would need to explicitly evaluate P_m(x) at the endpoints for every m (it is important to be clear with what one means by normalization in this case: orthonormal with respect to the l^2 norm or normalized with respect to the maximal value of P_m(x) on [-1,1] ).

The Legendre polynomials are orthogonal on the domain [-1,1] however we would like to use them on the domain of interest. We can define a mapping that will transform the domain of interest to that of the Legendre polynomials, \tilde{x} = \frac{2(x-x_{k-^1\!/_2})}{x_{k+^1\!/_2}-x_{k-^1\!/_2}}-1 to transform I_k to [-1,1] . Using this we can define a variant of the shifted Legendre polynomials such that \tilde{P}_m(x) = P_m(\tilde{x}(x)) . The newly defined shifted Legendre polynomials are now orthogonal over I_k  and the inner product becomes

(\tilde{P}_n,\tilde{P}_m) = \int_{x_{k-^1\!/_2}}^{x_{k+^1\!/_2}} \! P_n(\tilde{x}) P_m(\tilde{x}) \, \mathrm{d}\tilde{x} = \delta_{nm} \left(\frac{\Delta x}{2m+1}\right)

Using our orthogonal Legendre polynomial basis we can now construct an approximating solution u for arbitrary order M using basis functions \overset{m}{\psi}(x) \in \tilde{P}_m and basis weights \overset{m}{a}(t)

u = \sum_{m=0}^M \overset{m}{a}(t)\overset{m}{\psi}(x)

This approximation is well conditioned for arbitrary orders and has the additional benefit that it is easy to evaluate integrals and derivatives of it.

Where the orthogonal basis really shines however is when one needs to find the integral of the product of the approximation (or it’s derivative) with some other function (or it’s derivative) using the same orthogonal basis for the approximation. Usually the integrand of the two approximations can be calculated analytically and presents a very succinct closed form.

For example to solve the 1D the scalar conservation law PDE

\frac{\partial u}{\partial t} + \frac{\partial f(u)}{\partial x} = 0

could be solved using a N^{th} order spatial modal Discontinuous Galerkin method, so the the weak form becomes

\sum_{m=0}^M \left[ \frac{\mathrm{d}\overset{m}{a}}{\mathrm{d} t} \int_{I_k}\! \overset{m}{\psi} \, \overset{n}{\phi} \,\mathrm{d}x \right] + [g(v^-\!,v^+) \; \overset{n}{\phi}] \Big\rvert_{x_{k-^1\!/_2}}^{x_{k+^1\!/_2}} \!-\! \int_{I_k}\! f(v) \,\frac{\mathrm{d} \overset{n}{\phi}}{\mathrm{d} x} \,\mathrm{d}x = 0\;\;for \; 0 \leq n \leq N

Using an orthogonal Legendre basis approximation for both the solution and test function (a Galerkin approach) simplifies the expression greatly to

\frac{\mathrm{d}\overset{n}{a}}{\mathrm{d} t}\left(\frac{\Delta x}{2n+1}\right)= \!\! \sum_{\underset{by \: 2}{m=n-1}}^0 \overset{m}{a}_k-\sum_{m=0}^M \left[ \overset{m}{a}_k - \overset{m}{a}_{k-1}(-1)^m \right]\;\;\; for \; 0 \leq n \leq N

Which is then readily solved using an explicit method of lines approach for time discretization. Without an orthogonal basis approximation the integrands would all need to be calculated using quadrature and the routine would be far more expensive to perform.

2 thoughts on “Polynomial Approximations: Why monomials are bad”

  2. Hello,

    All stuff related to DG-FEM implementation for 1D convection equation is really helpful for me.
    Could you please explain the part of ExactRHS in matlab code?
    I am not getting it.

    Reply

Leave a Reply

Your email address will not be published. Required fields are marked *