Computation of shape functions

Shape functions are funcationals, approximating linear partial differential operators in given points. A shape function $\b{\chi}_{\mathcal{L}, \vec{p}}$ of a linear parial differential operator $\mathcal{L}$ in a point $p$ is such a functional that \[ (\mathcal{L} u)(\vec{p}) \approx \b{\chi}_{\mathcal{L}, \vec{p}} \b{u}, \] where $\b{u}$ is a vector of function values $u(x_i)$ evaluated in the support of $\vec{p}$. For more definitions and their usage see the section on Meshless Local Strong Form Method (MLSM).

Derivation

Let $\mathcal{L}$ be a linear partial differential operator. We wish to approximate $(\mathcal{L}u)(\vec{p})$ using only values of $u$ in support of $p$. To do so, we construct a MLS approximation (see Moving Least Squares (MLS)) and compute the coefficients $\b{\alpha}$, such that $\hat{u}(\vec{p}) = \b{b}(\vec{p})^\T \b{\alpha}$. The coefficients are computed by solving a system $WB \b{\alpha} = W\b{u}$ using the pseudoinverse \[ \b{\alpha} = (WB)^+W\b{u}, \] where $A^+$ denotes the Moore-Penrose pseudoinverse that can be calculated using SVD decomposition. Writing down the approximation function $\hat{u}$ we get \[ \hat u (\vec{p}) = \b{b}(\vec{p})^\T \b{\alpha} = \b{b}(\vec{p})^\T (WB)^+W\b{u} = \b{\chi}_{\vec{p}} \b{u}. \]

We have defined $\b{\chi}$ to be \[ \b{\chi}(\vec{p}) = \b{b}(\vec{p})^\T (WB)^+W. \] Vector $\b{\chi}$ is a row vector, also called a shape function. The name comes from being able to take all the information about shape of the domain and choice of approximation and store it in a single row vector, being able to approximate a function value from given support values $\b{u}$ with a single dot product. For any values $\b{u}$, value $\b{\chi}(p) \b{u}$ gives us the approximation $\hat{u}(\vec{p})$ of $u$ in point $\vec{p}$. Mathematically speaking, $\b{\chi}(\vec{p})$ is a functional, $\b{\chi}(\vec{p})\colon \R^n \to \R$, mapping $n$-tuples of known function values to their approximations in point $\vec{p}$.

The same approach works for any linear operator $\mathcal L$ applied to $u$. We approximate \[ (\mathcal{L}u)(\vec{p}) \approx (\mathcal{L} \hat{u})(p) = (\mathcal{L}\b{b})(\vec{p})^\T \b{\alpha} = (\mathcal{L}\b{b})(\vec{p})^\T (WB)^+W\b{u} = \b{\chi}_{\mathcal{L},\vec{p}} \b{u}. \]

For example, take a $1$-dimensional case for approximation of derivatives with weight equal to $1$ and $n=m=3$, with equally spaced support values at distances $h$. We wish to approximate $u''$ in the middle support point, just by making a weighted sum of the values, something like the finite difference \[ u'' \approx \frac{u_1 - 2u_2 + u_3}{h^2}. \] This is exactly the same formula as we would have come to by computing $\b{\chi}$, except that our approach is a lot more general. But one should think about $\b{\chi}$ as one would about the finite difference scheme, it is a rule, telling us how to compute the derivative. \[ u''(s_2) \approx \underbrace{\begin{bmatrix} \frac{1}{h^2} & \frac{-2}{h^2} & \frac{1}{h^2} \end{bmatrix}}_{\b{\chi}} \begin{bmatrix}u_1 \\ u_2 \\ u_3 \end{bmatrix} \]

The fact that $\b{\chi}$ is independent of the function values $\b{u}$ but depend only on domain geometry means that we can just compute the shape functions $\b{\chi}$ for points of interest and then approximate any linear operator of any function, given its values, very fast, using only a single dot product.

Numerical calculation of the shape functions

The expression \[ \b{\chi}_{\mathcal{L}, \vec{p}} = (\mathcal{L}\b{b})(\vec{p})^\T (WB)^+W \] can be evaluated directly, but this is not the most optimal approach. Numerically cheaper and more stable way is to translate the problem of inverting the matrix to solving a linear system of equations.

Invertible $B$ case: If $B$ is invertible, then $\b{\chi}_{\mathcal{L}, \vec{p}} = (\mathcal{L}\b{b})(\vec{p})^\T B^{-1}$, transposing the equation then multiplying it from the left by $B$, $\b{\chi}_{\mathcal{L}, \vec{p}}$ can be thought as a solution of a system $B^\T\b{\chi}_{\mathcal{L}, \vec{p}}^\T = (\mathcal{L}\b{b})$, which can be solved using LU or Cholesky decomposition for example.

General case: For a system written as $Ax = b$, where $A$ is a $n\times m$ matrix, $x$ is a vector of length $m$ and $b$ a vector of length $n$, a generalized solution is defined as $x$ that minimizes $\|A x - b\|_2^2$. If more $x$ attain the minimal value, $x$ with the minimal $\|x\|$ is chosen. Note that this generalizes the solution a general system ($A$ is invertible) and over-determined system ($n > m$ and $A$ has full rank). Such an $x$ can be computed using the pseudoinverse $x = A^{+} b$.

In our case, let us denote a part of the solution containing the pseudoinverse by $\tilde{\b{\chi}}_{\mathcal{L}, \vec{p}}$. \[ \b{\chi}_{\mathcal{L}, \vec{p}}(\vec{p}) = \underbrace{(\mathcal{L}\b{b})(\vec{p})^\T (WB)^+}_{\tilde{\b{\chi}}_{\mathcal{L}, \vec{p}}} W \] We have an expression $\tilde{\b{\chi}}_{\mathcal{L}, \vec{p}} = (\mathcal{L}\b{b})(\vec{p})^\T (WB)^+$ which after transposition takes the form $\tilde{\b{\chi}}_{\mathcal{L}, \vec{p}}^\T = ((WB)^\T)^+(\mathcal{L}\b{b})(\vec{p})$, the same as $x = A^+b$ above. Therefore, $\tilde{\b{\chi}}_{\mathcal{L}, \vec{p}}^\T$ is the solution of an (indeterminate) system $(WB)^\T \tilde{\b{\chi}}_{\mathcal{L}, \vec{p}}^\T = (\mathcal{L}\b{b})(\vec{p})$. After solving that, we can get the shape function $\b{\chi}_{\mathcal{L}, \vec{p}} = \tilde{\b{\chi}}_{\mathcal{L}, \vec{p}} W$ by multiplying by matrix $W$.

The system before can be solved using any decomposition of matrix $(WB)^\T = B^\T W$, not neceserally the SVD decompostion, but depending on our knowledge of the problem, we can use Cholesky ($B^\T W$ is positive definite), $LDL^\T$ if it is symmetric, $LU$ for a general square matrix, $QR$ for full rank overdetermined system and SVD for a general system. If more shapes need to be calculated using the same matrix $B^\T W$ and only different right hand sides, it can be done efficiently by storing the decomposition of $B^\T W$.

Another view on derivation of shape functions

Solving the system $B^\T\b{\chi}_{\mathcal{L}, \vec{p}}^\T = (\mathcal{L}\b{b})$ (in invertible $B$ case) has more than just algebraic meaning. It can be interpreted as using the method of undetermined coefficients. We wish to find such a combination of function values evaluated in support points that approximates the value of an operator $(\mathcal{L}u)(p)$, ie., we are looking for coefficients $\chi_i$, such that \[ \sum_{i=1}^n \chi_i u_i \approx (\mathcal{L}u)(p). \]

We require that $\chi_i$ be such, that above approximations holds on a certain set of functions, $\{b_j\}_{j=1}^m$ and we write the equations: \[ \begin{bmatrix} b_1(x_1) & \cdots & b_1(x_n) \\ \vdots & \ddots & \vdots \\ b_m(x_1) & \cdots & b_m(x_n) \\ \end{bmatrix} \begin{bmatrix} \chi_1 \\ \vdots \\ \chi_n \end{bmatrix} \begin{bmatrix} (\mathcal{L} b_1)(\vec{p}) \\ \vdots \\ (\mathcal{L} b_m)(\vec{p}) \end{bmatrix} \] The solution os the sought shape function $\b{\chi}$. By choosing the basis functions $b_j$ we decide for which functions the aproximation $\b{\chi}_{\mathcal{L}, \vec{p}} \b{u} \approx (\mathcal{L}u)(p)$ should be exact. If we do not specify enough functions $b_j$, we minimize $\|\b{\chi}_{\mathcal{L}, \vec{p}}\|$ instead.

Linearity and reduction to only computing the shapes of the operator space basis

Observe that the expression for the shape functions \[ \b{\chi}_{\mathcal{L}, \vec{p}} = (\mathcal{L}\b{b})(\vec{p})^\T (WB)^+W \] is linear in $\mathcal{L}$. Say that your equation constitutes of operators up to a certain $k$-th order. Then it is sufficient to compute only shape functions for the basis of that operator space, namely \[ \{\b{\chi}_{\frac{\partial}{\partial x^\alpha}, p}, 1 \leq |\alpha| \leq k \}, \text{ where } \frac{\partial}{\partial x^\alpha} = \frac{\partial^|\alpha|}{\partial x_1^{\alpha_1}\cdots x_d^{\alpha_d}} \text{ and } |\alpha| = \sum_{i=1}^d \alpha_i. \] Having the shape functions $\b{\chi}_{\frac{\partial}{\partial x^\alpha}, p}$ available, we write $\mathcal{L}$ in the above basis in point $p$, $\mathcal{L} = \sum_{1 \leq |\alpha| \leq k} a_\alpha(p) \frac{\partial}{\partial x^\alpha}$ and evaluate it's shape function as \[ \b{\chi}_{\mathcal{L}, \vec{p}} = \sum_{1 \leq |\alpha| \leq k} a_\alpha(p) \b{\chi}_{\frac{\partial}{\partial x^\alpha}, p}. \] This saves us space and gives certain elegance to the implementation.