How to robustly and numerically expand a $k$-order polynomial in two variables defined on a polygon domain?

Question

Given a $k$-order polynomial in two variable $p(x, y)$ defined on a polygon domain $K$. And I want to numerically expand it to the following form

$$ p(x, y) = c_0 + c_1 x + c_2 y + c_3 x^2 + c_4 xy + c_5 y^2 + .... $$

Let $\{\phi_i\}_{i=0}^{n_k-1}$ are all the basis with form $x^my^n$, here $m, n$ are non-negative integers and

$$ \begin{aligned} m + n &\leq k,\\ n_k &= \frac{(k+1)(k+2)}{2} \end{aligned} $$

I have tried to do it by $L^2$ projection, namely, solve the following linear equation

$$ \begin{bmatrix} (\phi_0, \phi_0)_K & (\phi_0, \phi_1)_K & \cdots & (\phi_0, \phi_{n_k-1})_K \\ (\phi_1, \phi_0)_K & (\phi_1, \phi_1)_K & \cdots & (\phi_1, \phi_{n_k-1})_K \\ \vdots & \vdots & \ddots & \vdots \\ (\phi_{n_k-1}, \phi_0)_K & (\phi_{n_k-1}, \phi_1)_K & \cdots & (\phi_{n_k-1}, \phi_{n_k-1})_K \\ \end{bmatrix} \begin{bmatrix} c_0 \\ c_1 \\ \vdots \\ c_{n_k-1} \end{bmatrix} = \begin{bmatrix} (p(x, y), \phi_0)_K \\ (p(x, y), \phi_1)_K \\ \vdots \\ (p(x, y), \phi_{n_k-1})_K \end{bmatrix} $$ where $(\phi_i, \phi_j)_K = \int_K\phi_i\phi_j\mathrm d x\mathrm d y$

But when the area of domain $K$ is very small, or the degree $k$ is bigger, the left matrix is nearly singular, then this system can not be solved robustly.

Any suggestions?

Answer 1

The problem that you are seeing is a well known problem of these basis functions:

$$\phi_{m,n}(x,y) = x^{m}y^{n}$$

I quote from here:

The reason why the coefficient matrix is nearly singular and ill-conditioned is that our basis functions $\Psi_{i}(x) = x^{i}$ are nearly linearly dependent for large i. That is, $x^{i}$ and $x^{i+1}$ are very close for i not very small. This phenomenon is illustrated in Figure 4. There are 15 lines in this figure, but only half of them are visually distinguishable. Almost linearly dependent basis functions give rise to an ill-conditioned and almost singular matrix. This fact can be illustrated by computing the determinant, which is indeed very close to zero (recall that a zero determinant implies a singular and non-invertible matrix): $10^{−65}$ for $N = 10$ and $10^{−92}$ for $N = 12$. Already for $N = 28$ the numerical determinant computation returns a plain zero.

So, your best bet is to use a set of orthogonal polynomials for your $L^{2}$ projection, such as Legendre, Hermite, etc, which I believe strongly depends on the type of the problem and its boundary conditions that you are trying to solve.

Let's say you have an orthogonal basis functions as: $\Phi_{m,n}(x,y)$ on $K$: ($\int_{K} \Phi_{m,n} \Phi_{m^{'},n^{'}} dx dy = \delta_{mm^{'}} \delta_{nn^{'}}$). As a result, your set of $c_{m,n}$ would be calculated as:

$$p(x,y) = \sum_{m}\sum_{n} c_{m,n}\Phi_{m,n}(x,y)$$

$$c_{m,n} = \int_{K} \Phi_{m,n}(x,y) p(x,y) dx dy$$

Without knowing more about your problem that lead to solve a $L^{2}$ projection, it's not possible to suggest which orthogonal polynomial basis is suitable for your application.