\section{A new scaling law for $m$ and $\ell$}
\label{sec:stab}

The selection of the degree $m$ of the PHS RBF requires more thought. On the one hand, $m$ is known to control the constant in the RBF-FD error estimate \eqref{eq:error_estimate}, with higher values of $m$ giving lower errors~\cite{FlyerPHS,BarnettPHS}. On the other hand, larger values of $m$ may cause greater instability, as the smallest eigenvalue of the block matrix from \eqref{eq:rbf_linsys} can be bounded from below by $q^m$, where $q$ is the separation distance of the stencil $P_1$; see~\cite{Wendland:2004} for some preliminary estimates, and~\cite{FlyerNS} for remarks about stability. A third more subtle issue is the question of whether $m$ and $\ell$ should be related. Historically, it was common to use the \emph{classical scaling law}~\cite{iske2002}:
\begin{align}
m = 2\ell+1,
\label{eq:scaling_law_classical}
\end{align}
where $\ell$ is selected to be the minimum to prove unisolvency for the RBF interpolant with $\phi(r) = r^m$. This approach was also recently used for PHS-based methods for advection on the sphere~\cite{ShankarWrightJCP2017}, and advection and reaction-diffusion equations on more general manifolds~\cite{SNKJCP2018}. However, modern PHS RBF-FD methods for Euclidean domains decouple $m$ and $\ell$~\cite{FlyerNS,FlyerPHS,BarnettPHS,ShankarJCP2017}. Our goal is to attempt to find a scaling law that produces the smallest spurious eigenvalue (if any) in the spectrum of the discrete Laplacian $L$. To that end, we explored two alternatives to the classical scaling law. The first alternative is:
\begin{align}
m = 
\begin{cases}
	\hfill \ell \hfill & \text{if $\ell$ is odd} \\
	\hfill \ell+1 \hfill & \text{if $\ell$ is even}\\
\end{cases}
\label{eq:scaling_law_1}
\end{align}
which serves to gently increase $m$ as $\ell$ is increased. The second alternative, which we ultimately settled on, is:
\begin{align}
m = 
\begin{cases}
	\hfill \ell \hfill & \text{if $\ell$ is odd} \\
	\hfill \ell-1 \hfill & \text{if $\ell$ is even}\\
\end{cases}
\label{eq:scaling_law_2}
\end{align}
The eigenvalues for the discrete Laplacian $L$ computed using the classical scaling law, these two alternatives, and the choice of $m=3$ are shown for $N=658$ nodes on the disk using $\ell=6$ in Figure \ref{fig:eigs_laws1}.
\begin{figure}[h!]
\centering
\subfloat[$m=2\ell+1$]
{
	\includegraphics[scale=0.6]{"Figures/Eig2lp1".png}
	\label{fig:c11}	
}
\subfloat[$m=\ell+1$]
{
	\includegraphics[scale=0.6]{"Figures/Eig7".png}
	\label{fig:c12}	
}

\subfloat[$m = \ell-1$]
{
	\includegraphics[scale=0.6]{"Figures/Eiglm1".png}
	\label{fig:c21}	
}
\subfloat[$m = 3$]
{
	\includegraphics[scale=0.6]{"Figures/Eig3".png}
	\label{fig:c22}	
}
\caption{Spectrum of the discrete Laplacian $L$ on a node set with $N=658$ nodes for different relationships between $m$ and $\ell$. The node set here was a set of Poisson disk samples generated on the disk using the method from~\cite{SFKSISC2017}.}
\label{fig:eigs_laws1}
\end{figure}
From Figure \ref{fig:c11}, it is clear that the classical scaling law $m=2\ell+1$ results in a large spurious eigenvalue. In general, all the other scaling laws appear to give very similar results, albeit with some subtle differences. Figure \ref{fig:c12} uses the scaling law $m=\ell+1$, while Figure \ref{fig:c21} uses the scaling law $m = \ell-1$; clearly, both produce almost identical results, with the exception that $m=\ell+1$ results in a slightly larger spurious eigenvalue. Continuing this trend, it may be tempting to always simply choose a very small reasonable value of $m$. However, Figure \ref{fig:c22} with $m=3$ shows that this may not be entirely advisable either, since the largest spurious eigenvalue here is actually larger than the cases of $m=\ell+1$ and $m=\ell-1$. 
\begin{figure}[h!]
\centering
\subfloat[$m=2\ell+1$]
{
	\includegraphics[scale=0.6]{"Figures/Eig2lp1_2".png}
	\label{fig:d11}	
}
\subfloat[$m=\ell+1$]
{
	\includegraphics[scale=0.6]{"Figures/Eiglp1_2".png}
	\label{fig:d12}	
}

\subfloat[$m = \ell-1$]
{
	\includegraphics[scale=0.6]{"Figures/Eiglm1_2".png}
	\label{fig:d21}	
}
\subfloat[$m = 3$]
{
	\includegraphics[scale=0.6]{"Figures/Eig3_2".png}
	\label{fig:d22}	
}
\caption{Spectrum of the discrete Laplacian $L$ on a node set with $N=919$ nodes for different relationships between $m$ and $\ell$. The node set here was a set of Poisson disk samples generated on the disk using the method from~\cite{SFKSISC2017}.}
\label{fig:eigs_laws2}
\end{figure}
The spectrum of $L$ for different scaling laws for $N=919$ nodes with $\ell=6$ is shown in Figure \ref{fig:eigs_laws2}. Once again, $m=2\ell+1$ can be rejected as extremely unstable, even after refinement (Figure \ref{fig:d11}). However, the gaps between the other scaling laws have widened slightly as well. For instance, $m=\ell+1$ (Figure \ref{fig:d12}) now results in a spurious eigenvalue that is twice as large as $m=\ell-1$ (Figure \ref{fig:d21}). Even more interestingly, $m=3$ (Figure \ref{fig:d22}) results in a spurious eigenvalue whose real part is much larger than the $m=\ell+1$ and $m=\ell-1$ cases. Once again, the most stable choice appears to be $m = \ell-1$. Picking a fixed value of $m=3$ (or any other number) regardless of the choice of $\ell$ may be inadvisable on Euclidean domains except when the node sets are very regular. These trends continue for different values of $\ell$ (not shown). It may be the case that the correct scaling relationship needs to be reevaluated for each type of node set used.

Thus far, we have only discussed the \emph{global} stability associated with the different scaling laws. It is reasonable to expect that these global patterns are reflected in some way at the stencil level. In other words, we expect that the global stability must be connected to \emph{local} stability in some way.  Before we can explore this connection, we must first define a few quantities. Returning to the stencil $P_1$, consider the \emph{local $\calL$-Lebesgue function} of the augmented RBF interpolant $s_1(\vx,\vy)$ used for approximating the operator $\calL$, defined as
\begin{align}
\Lambda^1_{\calL}(\vy) = \sum\limits_{j=1}^n |w^1_j(\vy)| = \|{\bf w}^1\|_1,
\end{align}
that is, the 1-norm of the RBF-FD weight vector at a point. Let $w^1_1(\vy)$ be the RBF-FD weight associated with the stencil point $\vy$. From~\cite{ShankarJCP2017}, we have the sufficient condition for all eigenvalues of the differentiation matrix $L$ (corresponding to $\calL$) having non-positive real parts as:
\begin{align}
2w^1_1(\vy) \geq - \Lambda^1_{\calL}(\vy).
\end{align}
This condition highlights the connection between the local Lebesgue function (itself connected to the stability of the local RBF interpolant~\cite{Wendland:2004}) and the \emph{global} stability of the RBF-FD method. In general, large values of $\Lambda^1_{\calL}(\vy)$ are associated with instability, and typically produce spurious eigenvalues with positive real parts~\cite{ShankarJCP2017}. To examine the impact of the scaling laws \eqref{eq:scaling_law_classical}--\eqref{eq:scaling_law_2} on the local Lebesgue functions (and therefore local stability), we will now study the local $\Delta$-Lebesgue functions for a set of $N = 4000$ Halton points in $[-1,1]^2$, where $\Delta$ is the Laplacian; evenly-spaced points are used on the boundary; these Halton points are typically even more scattered than Poisson disk samples, and represent a worst-case scenario. The results are shown in Figure \ref{fig:leb_c1}.
\begin{figure}[h!]
\centering
\subfloat[]
{
	\includegraphics[scale=0.4]{"Figures/LocalStab_9_4".pdf}
	\label{fig:leb_c13}	
}
\subfloat[]
{
	\includegraphics[scale=0.4]{"Figures/LocalStab_5_4".pdf}
	\label{fig:leb_c11}	
}
\subfloat[]
{
	\includegraphics[scale=0.4]{"Figures/LocalStab_3_4".pdf}
	\label{fig:leb_c12}	
}

\subfloat[]
{
	\includegraphics[scale=0.4]{"Figures/LocalStab_13_6".pdf}
	\label{fig:leb_c23}	
}
\subfloat[]
{
	\includegraphics[scale=0.4]{"Figures/LocalStab_7_6".pdf}
	\label{fig:leb_c21}	
}
\subfloat[]
{
	\includegraphics[scale=0.4]{"Figures/LocalStab_5_6".pdf}
	\label{fig:leb_c22}	
}

\subfloat[]
{
	\includegraphics[scale=0.4]{"Figures/LocalStab_17_8".pdf}
	\label{fig:leb_c33}	
}
\subfloat[]
{
	\includegraphics[scale=0.4]{"Figures/LocalStab_9_8".pdf}
	\label{fig:leb_c31}	
}
\subfloat[]
{
	\includegraphics[scale=0.4]{"Figures/LocalStab_7_8".pdf}
	\label{fig:leb_c32}	
}
\caption{Local $\Delta$-Lebesgue functions for different scaling laws on a logarthmic scale. The figure contrasts \eqref{eq:scaling_law_classical} (left), \eqref{eq:scaling_law_1} (middle), and \eqref{eq:scaling_law_2} (right). Lighter colors indicate larger values, which in turn indicate greater potential for instability.}
\label{fig:leb_c1}
\end{figure}
The first column in Figure \ref{fig:leb_c1} corresponds to the classical law \eqref{eq:scaling_law_classical}, the second to \eqref{eq:scaling_law_1}, and the third to \eqref{eq:scaling_law_2}. The colorbar scale is fixed across each row, but allowed to increase between rows to make clear that RBF-FD weights can increase in magnitude as the polynomial degree is increased. Immediately, it is clear that the classical scaling law results in the largest weights, and is consequently almost always unstable (Figures \ref{fig:leb_c13},\ref{fig:leb_c23}, \ref{fig:leb_c33}).  Careful observation also shows that the last column (corresponding to \eqref{eq:scaling_law_2}) results in more stable approximations than the middle column (corresponding to \eqref{eq:scaling_law_1}), which matches our conclusions from the eigenvalue plots (Figures \ref{fig:eigs_laws1} and \ref{fig:eigs_laws2}). In general, we find that the scaling law \eqref{eq:scaling_law_2} produces stable approximations across a wide range of parameters, and therefore we use it in this article.
{"config": "arxiv", "file": "1806.03798/Stability.tex"}