View source for 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}.
\]
This saves us space and gives certain elegance to the implementation.