How does one approximate a second derivative with ATPS interpolation

Question

When using the Dual Reciprocity Boundary Element Method ( or any radial basis function method ) to solve a nonlinear differential equation it is necessary to approximate some derivatives of a potential field using radial basis functions.

If you have a potential field $\mathbf{u}_i$ with coordinates $\mathbf{x}_i=(x_i,y_i)$. The field can be approximated with radial basis functions such that $\mathbf{u}=\mathbf{F}\mathbf{\alpha}$. Where $\mathbf{F}$ is a matrix of radial basis functions based on the following: $$ f(\mathbf{x_i})=\Sigma_{j=1}^{N}\alpha_j]\phi(||\mathbf{x_i}-\mathbf{x_j}||)_2$$ $$ \mathbf{F}=\left[ \begin{matrix} \phi(||\mathbf{x_1}-\mathbf{x_1}||)_2 & \cdots & \phi(||\mathbf{x_1}-\mathbf{x_N}||)_2 \\ \vdots & \ddots & \vdots \\ \phi(||\mathbf{x_N}-\mathbf{x_1}||)_2 & \cdots & \phi(||\mathbf{x_N}-\mathbf{x_N}||)_2 \\ \end{matrix}\right] $$ $$ \mathbf{\alpha}=\left[\begin{matrix} \mathbf{\alpha_1} \\ \vdots \\ \mathbf{\alpha_N} \\ \end{matrix}\right]$$

For ATPS (Augmented Thin Plate Splines) the matrix of radial basis functions becomes: $$ \mathbf{x_i}=(x_i,y_i) $$ $$ r=\sqrt{ (x_i-x_j)^2 + (y_i+y_j)^2 }=||\mathbf{x_i}-\mathbf{x_j}|| $$ $$f(\mathbf{x_i})=\Sigma_{j=1}^{N}\alpha_j r^2 log(r) + \beta_1+\beta_2x_j+\beta_3y_j $$

$$\Sigma_{j=1}^N \alpha_j=\Sigma_{j=1}^N \alpha_j x_j=\Sigma_{j=1}^N \alpha_j y_j=0 $$

$$ \mathbf{P}= \left[ \begin{matrix} 1 & 1 & \cdots & 1 \\ x_1 & x_2 & \cdots & x_N \\ y_1 & y_2 & \cdots & y_N \\ \end{matrix} \right] $$ $$ \mathbf{F^*}= \left[ \begin{matrix} \mathbf{F} && \mathbf{P^T} \\ \mathbf{P} && \mathbf{0} \\ \end{matrix} \right] $$

$$ \left[ \begin{matrix} \mathbf{u} \\ \mathbf{0} \\ \end{matrix} \right]= \left[ \begin{matrix} \mathbf{F} && \mathbf{P^T} \\ \mathbf{P} && \mathbf{0} \\ \end{matrix} \right]\left[ \begin{matrix} \mathbf{\alpha} \\ \mathbf{\beta} \\ \end{matrix} \right]=\mathbf{F^*\alpha^*} $$

As $\alpha$ and $\beta$ are constants, a spatial derivative of $\mathbf{u}$ becomes a spatial derivative of the radial basis function series. Thus:

$$ \frac{\partial\mathbf{u}}{\partial x}=\frac{\partial\mathbf{F^*}}{\partial x}\mathbf{\alpha^*}=\frac{\partial\mathbf{F^*}}{\partial x}\mathbf{{F^*}^{-1}u} $$

similarly the second derivative would be $$ \frac{\partial^2\mathbf{u}}{\partial x^2}=\frac{\partial^2\mathbf{F^*}}{\partial x^2}\mathbf{{F^*}^{-1}u} $$

However when taking the second derivative of the ATPS function I obtain: $$\frac{\partial^2f(\mathbf{x_i})}{\partial x^2}=\Sigma_{j=1}^{N}\alpha_j \left[2log(r)+\frac{(y_i-y_j)^2}{r^2}+1 \right] $$

If you then take the limit of this as $\mathbf{x_i}$ approaches $\mathbf{x_j}$ you obtain $-\infty$ for the diagonals of the matrix $\frac{\partial^2\mathbf{F^*}}{\partial x^2}$.

Thus i do not understand how to obtain a second derivative approximation of a potential field u using Augmented Thin Plate Splines.

Answer 1

It is worth remembering that approximation theory is a branch of numerical analysis. In numerical analysis when you're dealing with negligible terms, you don't look closely into how smooth those terms are: once they are sufficiently close to zero they are zero.

As you program TPS evaluation and want to avoid a runtime error you shouldn't attempt evaluating $r^2 \log r$ at $r=0$: although the limit value is zero, part of that expression is still $\log r$ which is simply undefined at $0$.

Instead, TPS is programmed as a piecewise function $\cases{r^2\log r, r>\varepsilon \\0}$.

So what you're trying to calculate is the second derivative of a constant $0$. It is zero.

Numerics aside, in pure analysis, we are not even allowed to differentiate a functions where its undefined; and if you make limit a part of definition, then you won't be allowed to just swap differentiation and limit at the singularity.

Answer 2

A "classic TPS spline $(\sigma)$ is a continuously differentiable function $\sigma \in C^1(\mathbb{R}^2)$". If you need to approximate second derivatives of the function then you have to use a twice continuously differentiable version of the Duchon spline.