Ill-conditioned Gram Matrix Assembly

Question

I'm trying to find the best approximation to the function $e^x$ in the finite dimensional polynomial space $P_4$ with respect to the standard basis vectors $B=\{1,x,x^2,x^3,x^4\}$ with inner product $$(f,g)=\int_0^1 \frac{f(x)g(x)}{\left(x(1-x)\right)^{\frac{1}{2}}}dx$$

Let $w=\sum_{j=0}^4a_j\phi_j $ where $\phi_j=x^j$. Then, the best approximation of $w$ to $e^x$ satisfies the orthogonality property:

$$(e^x-w,\phi_i)=0$$ for $i=0,...,4$. Hence,

$$\sum_{j=0}^4 a_j (\phi_j,\phi_i) = (e^x,\phi_i)$$ for $i=0,...,4.$. This generates linear system of equations which I'm trying to solve. Unfortunately, the integrals are a bit difficult to evaluate by hand. So, I tried using a three point gaussian quadrature rule to approximate the integrals and solve the system. Unfortunately, I keep getting a nearly singular system. When I observed this, I tried using a higher order quadrature rule, but the matrix remains ill-conditioned.

I looked at my code over and over again and I can't seem to find any errors in the implementation. I'm tempted to believe that these inner products simply can't be evaluated well with a gaussian quadrature rule. Is this common for Gram Matrices? Should I use a different quadrature rule? Or is it a bug in my code?

I have provided my matlab code assembing this matrix system below for reference. Any help would be greatly appreciated.

n=5;
%Quadrature points on [-1,1]
z1=-sqrt(3/5);
z2=0;
z3=sqrt(3/5);

%Quadrature points on [0,1]
x1=z1/2+1/2;
x2=z2/2+1/2;
x3=z3/2+1/2;

A=zeros(n,n);
b=zeros(n,1);

%assemble gramm matrix
for i=1:n
    for j=1:n
        F1=x1^(i-1+j-1)/sqrt(x1*(1-x1));
        F2=x2^(i-1+j-1)/sqrt(x2*(1-x2));
        F3=x3^(i-1+j-1)/sqrt(x3*(1-x3));
        A(i,j)=(1/2)*((5/9)*F1+(8/9)*F2+(5/9)*F3);
    end
    G1=x1^(i-1)*exp(x1)/sqrt(x1*(1-x1));
    G2=x2^(i-1)*exp(x2)/sqrt(x2*(1-x2));
    G3=x3^(i-1)*exp(x3)/sqrt(x3*(1-x3));
    b(i)=(1/2)*((5/9)*G1+(8/9)*G2+(5/9)*G3);
end

Answer 1

Preliminary observations

• If you have a basis $\{\phi_i(x)\}$ of size $n$, you must have at least $n$ quadrature points. You can see this by writing the inner product as $(u,v) = u^T B^T W B v$ where $W$ is a diagonal weight matrix and $B_{ij} = \phi_j(x_i)$ is the basis evaluation matrix. If there are fewer than $n$ quadrature points, $B$ has fewer than $n$ rows, thus $\DeclareMathOperator{\rank}{rank} \rank(B^T W B) < n$.

• Your weight function cannot be integrated exactly using Gauss quadrature, depending on the order and desired accuracy, may have unacceptable errors. You can find papers on special quadratures for all sorts of specific functional forms. Note that there is also quadrature error integrating the right hand side. Fortunately, these functions are smooth and the weight effectively asks for more points near the endpoints, so Gauss quadrature is not that bad. See Trefethen (2008) Is Gauss Quadrature Better Than Clenshaw-Curtis for more on generic quadratures.

• The representation $A = B^T W B$ is a more convenient way to write your code.

• It is considered good Matlab style to use matrix and vector notation when convenient. It substantially compresses and simplifies the code.

Refactored and working code

n = 5;
qpoints = 10;

% Quadrature points on [-1,1]; see Golub & Welsch (1969) and Trefethen (2008)
beta = .5 ./sqrt(1-(2*(1:qpoints-1)).^(-2)); % 3-term recurrence coefficients
T = diag(beta,1) + diag(beta,-1);            % Jacobi matrix
[V,D] = eig(T);                              % eigenvalue decomposition
[x,i] = sort(diag(D));                       % Gauss-Legendre points
w = 2*V'(i,1).^2;                            % weights
% Map quadrature to interval [0,1]
x = 0.5 + 0.5*x; w = 0.5*w;

% assemble Gram matrix
B = fliplr(vander(x,n));         % basis evaluation matrix
W = diag(w ./ sqrt(x.*(1-x)));   % diagonal quadrature weight matrix
A = B'*W*B;                      % Gram matrix
b = B'*W*exp(x);                 % exponential integrated against test functions

taylor = 1./factorial(0:n-1)';
ceoffs = [A \ b, taylor]

Further observations

• The monomial basis is notoriously ill-conditioned. If you would like to work with higher degree polynomials, you should choose a well-conditioned basis such as (a) Legendre/Chebyshev polynomials, (b) Lagrange interpolants on good point sets like Gauss-Lobatto, or (c) Newton-form (divided differences) polynomials associated with a good point set. Of these, (b) is the most popular in PDE computations.
• The computed coefficients do not match the leading terms in the exponential because they are weighted-$L^2$ optimal on the interval $[0,1]$ rather than only in the limit as $x \to 0$. In contrast, the Taylor polynomials only apply to the $x\to 0$ limit and provide no bounds on a finite interval in any norm.
• You can produce optimal polynomials with respect to non-smooth norms such as $L^\infty(\Omega)$. Approximation in $L^\infty$, sometimes called "Chebyshev approximation", gives you optimal pointwise error bounds across the entire interval. This blog post has some nice figures, an introduction to the algorithms, and further links to the literature.

Answer 2

The gram matrix of the powers is known to be very ill-conditioned. It is far better if you expand in terms of Chebyshev polynomials for an integral in [-1,1], and transform your integral to that range.