Role of weight function in Galerkin methods

Question

I have difficulties in understanding the role of the weight function $w(x)$ that occurs in the solution of PDEs via the Galerkin approach. Consider a linear differential equation of the form $$ \partial_t \,u(x,t) \ = \ L \, u(x,t) $$ where $L= L[x,\partial_x,\partial_{xx}]$ is a differential operator. In the Galerkin approach one uses the ansatz $u(x,t) = \sum_k a_k(t) \,p_k(x)$ to obtain $$ \sum_k {\dot a}_k(t) \,p_k(x) \ = \ \sum_k a_k(t) L\, p_k(x) $$ The next step is an application of the functional $\int p_i(x) w(x) dx$: $$ \sum_k {\dot a}_k(t) \,(p_i,p_k)_w \ = \ \sum_k a_k(t) \;(p_i, L \, p_k)_w $$ The resulting equation is solved by the usual methods for time-propagation.

Now to the point: the functions $p_k(x)$ are often chosen as orthogonal polynomials (with appropriately built-in boundary conditions via basis recombination). For this, the book of Hesthaven and Gottlieb states (e.g. page 118)

Since we are using (orthogonal) polynomials $p_k(x)$ as a basis it is natural to choose the weight function $w(x)$ such that $(p_k,p_j)_w =\gamma_k \delta_{kj}$

Questions:

1. Isn't it so that for each different $w(x)$, one solves a different problem?

2. Why should one a priory attribute particular importance to different regions of $[-1,1]$ and relatively neglect others (as it is done for all weights $w(x)\neq1$). Put differently, why should one ever use the Chebyshev polynomials and weight function $w_\text{Cheb}(x) = 1/\sqrt{1-x^2}$ if I can use the Legendre polynomials and $w_{\text{Leg}}(x)=1$ (--one reason is the FFT, but please neglect that here, as it's meant more as a general question)?

3. And if one does so, i.e. applies a weighting in $[-1,1]$, which weight function is "optimal" (in the sense, that the original PDE is approximated "best")?

4. And, just to confirm my understanding: if I use Chebyshev polynomials with a weight $w(x)=1$ (which is effectively Clenshaw-Curtis integration), I lose all the advantages of the Gaussian grid, yes?

I'd appreciate also partial answers, thanks in advance.

Answer 1

Chebyshev are orthogonal wrt to a weight function which is singular at the end-points. When you approximate a function f(x) with Chebyshev, the convergence of the approximations is not affected by the values of f or its derivatives at the end points. The convergence rate depends only on the smoothness of the function f and not on its boundary values. This is a general property of any orthogonal basis that comes as eigenfunctions of a sturm-liouville problem with singular weight functions. See ([1], section 3). This would not be true with uniform weight.

[1] Gottlieb and Orszag: Numerical analysis of spectral methods

Answer 2

First, Galerkin in his article from 1915 does not discuss any weighting other the $\omega(x) = 1$. The Galerkin method is a direct generalization of the Rayleigh-Ritz method, and a variational procedure cannot be constructed with any other weight. That approach was the only one used until the FFT came about. If you want to use FFT, then you must use Chebyshev polynomials and hence the weighting $\omega(x) = 1/\sqrt{1-x^2}$. I would call it a Chebyshev-Galerkin method to distinguish it form a true Galerkin method. But to answer your first question, yes, you are using a different approximation. Everyone seems to get hung up on the basis. What you need to look at is not the basis, but the residual of the differential equation. We know from uniqueness that you can use any basis you want - orthogonal polynomials, Lagrange polynomials or monomials. You'll get the same answer, just in a different form. I'll simplify your equation further to a BVP: $Lu = f$. The residual is $R(x,\boldsymbol{a}) =Lu - f$. Let's also say we have homogenous Dirichlet boundary conditions on the interval [-1,1]. Since the Galerkin method makes the residual orthogonal to the basis functions, it does in effect use a weighting of $1-x^2$, because the basis functions obey the boundary conditions which are zero at both ends. This sort of makes sense, because you know the solution at the boundaries, so it should receive less weight. If you use a weighting of unity it would be the method of moments which is equivalent to the Legendre-tau method. In all cases the residual is made orthogonal to the first $n$ polynomials with some weighting. So, the residual can be represented by a series of Jacobi polynomials: $$R(x,\boldsymbol{a}) =\sum_{i=n}^{m} \tau_i P_{i}^{(\alpha,\alpha)}(x)$$ where $n$ is the number of basis functions, $m$ is some upper limit which depends on the differential equation and $\alpha = 0,\frac12,1$ for moments, Chebyshev-Galerkin, and Galerkin methods respectively. You could also consider a Chebyshev-tau or Chebyshev-moments method with $\alpha = -\frac12$. I recently published a paper Orthogonal Collocation Revisited (https://doi.org/10.1016/j.cma.2018.10.019) where I compute the values of $\tau$ for a simple problem. There are a few studies of the tau method which have looked at which choice of $\alpha$ is better, but mostly you have numerical experiments. You will find a more detailed discussion and other examples in the unfinished monograph at my website (link below). Which I’ve updated recently.

Answer 3

The basis functions $\{p_n(x): n \in \mathbb{N}\}$ can be orthogonal in $[a, b]$ with respect to a weight function $w(x)$. For example, Hermite polynomials are orthogonal in $(-\infty, \infty)$ with respect to $e^{-x^2}$,

$$\int_{-\infty}^\infty H_m(x) H_n(x)\, e^{-x^2}\, \mathrm{d}x = \sqrt{ \pi}\, 2^n n! \delta_{nm}\, .$$

The weight function is needed for the orthogonality. In this case the integral would not even be bounded without the weight function.

Update

Let's just consider [-1,1], I'll repeat the question: how could one ever come up with a different weight function than w(x)=1? What is the motivation a priori?

That's an interesting question, but this is somewhat the reverse thought that one normally has. I would say that you don't want to have a different (from $w(x)=1$) weight function, it just happens to be the case that for some set of functions in a certain vector space to be orthogonal you have to have it. There might be several reasons for picking one set of functions over other as basis to approximate the solution of a PDE.

The first time I ended up using orthogonal functions was when I was using separation of variables to solve PDEs, where was common to write the spatial part as a Sturm-Liouville. When you have a differential equation written in this form you know that you have a set of orthogonal functions with a certain weight.

Answer 4

I will answer you in order:

1) The problem to solve is independent of the weighing function you choose, e.g. you could solve for a PDE using standard finite element basis (Lagrange) without using special functions that has special properties (orthogonal functions).

2) The weight function $w(x)$ for certain polynomials to be orthogonal arise naturally, e.g. the functions $f(x) = 1$ and $g(x) = x$ are orthogonal in the interval $[-1,1]$ w.r.t. the weights $w(x) = 1$ and $w(x) = 1/\sqrt{1-x^2}$, but if sequence continues this is no longer true. For example, basis functions like $f_n(x) = \cos(nx)$ are orthogonal for different $n$ in $[0,2\pi]$: $$\int_0^{2\pi}{f_nf_mdx}=\int_0^{2\pi}{\cos{(nx)}\cos{(mx)}dx}=\pi\delta_{nm}$$ This functions, are orthogonal w.r.t. the weight function $w(x) = 1$. But, if we change variables to $y=\cos{x}$ you will discover actually that: $$\int_0^{2\pi}{\cos{(nx)}\cos{(mx)}dx}=\int_{-1}^{1}\frac{T_n(y)T_m(y)}{\sqrt{1-y^2}}dy = \pi\delta_{nm}$$ Being $T_n(x) = \cos{(n\cos^{-1}{(y)})}$

3) The optimal basis is a basis that spans the whole space (your interval) and its "basis vectors" are as independent as possible among them. The perfect basis is of course an orthogonal one (think about cartesian axis, where everything is easy e.g. scalar products). If you choose this kind of basis, your solution will be best representated. There micht be reasons to choose a basis before another, e.g. Legendre polynomials are preferred in situations when extremal points of the interval are of importance and thus they are nodes.

4) The only advantage of the use of an orthogonal polynomial is that the scalar product leads to a diagonal matrix mass matrix, and inversion is then trivial.