# Least SQuares Discretizations (LSQD): a robust and versatile meshless
paradigm for solving elliptic PDEs

Anna Kucherova  Gbocho M. Terasaki  Selma Strango  Maxime Theillard

###### Abstract

Searching for numerical methods that combine facility and efficiency, while
remaining accurate and versatile, is critical. Often, irregular geometries
challenge traditional methods that rely on structured or body-fitted meshes.
Meshless methods mitigate these issues but oftentimes require the weak
formulation which involves defining quadrature rules over potentially
intricate geometries. To overcome these challenges, we propose the Least
Squares Discretization (LSQD) method. This novel approach simplifies the
application of meshless methods by eliminating the need for a weak formulation
and necessitates minimal numerical analysis. It offers significant advantages
in terms of ease of implementation and adaptability to complex geometries. In
this paper, we demonstrate the efficacy of the LSQD method in solving elliptic
partial differential equations for a variety of boundary conditions,
geometries, and data layouts. We monitor h-P convergence across these
parameters/settings and construct an a posteriori built-in error estimator to
establish our method as a robust and accessible numerical alternative.

###### keywords:

meshless , least squares , h-P convergence , arbitrary geometries , complex
boundary conditions

††journal: Journal of Computational Physics

##  1 Introduction

In the ever-evolving landscape of computational science and engineering, the
quest for accurate, efficient, and versatile numerical methods remains
paramount. Traditional numerical approaches, such as finite element methods,
finite difference methods, and finite volume methods, have long dominated the
sphere of simulating complex physical phenomena. However, their reliance on
elements, structured grids, or finite volumes often presents challenges when
handling irregular geometries, dynamic domains, or intricate boundary
conditions. In response to these limitations, a trailblazing paradigm emerged:
meshless methods [18, 35].

Gaining traction almost two decades post-emergence [1, 2, 4, 5, 27, 34],
meshless methods have witnessed a surge in popularity across a multitude of
disciplines, including computational fluid dynamics, structural mechanics,
electromagnetics, and biomedical engineering [39, 31, 22, 26]. Their inherent
attributes (e.g. mesh generation simplicity, capacity to handle large
deformations/distortions and complex geometries [32], etc.), have proven a
versatile alternative to the traditional discretization techniques.

Some of the pioneering meshless methods fall into the following categories:
kernel-based [18, 34, 57], Galerkin [4, 6, 12, 1, 54, 2, 19, 56, 20, 10, 3] ,
and Finite Difference methods [33, 42]. Meshless Galerkin methods such as the
Element-free Galerkin [4, 6], h-p clouds and h-p methods [12], Meshless Local
Petrov-Galerkin [1, 2], the Natural Element Method [54, 19, 56], the Voronoi
finite volume method [20, 10], and partition of unity method [3], all rely on
the weak formulation of the partial differential equation (PDE). Thus,
although they are meshless, these methods are not truly “geometry-free” as the
weak formulation involves spatial integrals. This requires approximations of
definite integrals, or quadrature rules, over - potentially - intricate
geometries. These geometries, themselves, may reap extravagant computational
costs (e.g. Voronoi partition for the Meshless Natural Element method [54, 19,
56] or any Voronoi Finite Volume Method [20, 10]). Undeniably, the theoretical
foundation for the aforementioned methods allows for a comprehensive
mathematical study but may serve as a deterrent for interested non-experts.
Clearly, it would be a quixotic undertaking to create an encyclopedic
inventory for the trajectory of all meshless methods, so for an exhaustive
analysis, we direct interested readers to the following reviews [41, 44, 13].

Although commonly used in scattered data interpolation [15], Kansa became one
of the first to utilize radial basis functions (RBF) to solve fluid dynamics
problems [27]. Theoretical condition number estimates indicated that some
meshless methods using RBFs may be unstable [11]. Fortunately, Kansa was able
to successfully circumvent the ill-conditioned system by using improved multi-
quadrics and adaptive techniques [28]. Kansa’s method represents the solution
as a sum of global radial functions, centered at a set of points. The
coefficients in this decomposition are obtained by interpolating the PDE at
these points. While in general, the interpolation system can be ill-posed
(Mairhuber-Curtis theorem [38, 9]), when using appropriate radial functions,
this system can be proven to be non-singular [40, 46]. Not solely is it
straightforward to implement, this method is particularly well-suited for
problems involving complex curved surfaces (e.g. such as reaction-diffusion
patterns over biological entities such as frog species [45], red blood cells
[16], or 3-D spherical shells [55] used in geophysical calculations). It has
since been applied in the computational fluid dynamics realm [48] such as the
solution of Navier–Stokes equations [36], porous media flow [49], and the
solution of solid-liquid phase change problems [29]. Unfortunately, the
convergence of Kansa’s method is still unproven. Because the representation is
global, the interpolation system is full and consequently its inversion
rapidly becomes prohibitive. To address this limitation, domain decomposition
techniques can be employed [37], but these abandon the meshless nature of the
domain. Using local RBFs [21], the local Kansa method [49, 50] offers an
alternative approach to reduce the complexity. However, this approach still
relies on a domain decomposition to distinguish between inner and border
nodes, further requiring functional iterations on either node sets.

One of the most pervasive methods in computational fluid dynamics is the
finite point method - a meshfree method for solving PDEs on scattered
distributions of points. This method is considered a generalization of the
finite difference method, created for arbitrary irregular grids [33, 8].
Despite the fact that this method requires a much lower mathematical
background, it entails a more tedious implementation as the discretization
coefficients are obtained from high-order Taylor expansions over non-uniform
neighborhoods [42]. With the Radial Basis Function Finite Difference (RBF-FD
[14]), these coefficients are computed at each point by interpolating the
solution over a small neighborhood using a radial basis enriched with
polynomial functions. The RBF ensures good conditioning while the polynomial
enrichment grants convergence. The radial basis function-based differential
quadrature approach (RBF-DQ [52]) proposes a similar strategy and has been
employed to simulate compressible flows [52] and heat conduction [53]
problems.

Despite offering significant advantages (e.g. flexibility in handling complex
geometries, avoiding computational overhead associated with mesh generation),
meshless methods possess considerable limitations. One notable limitation is
that most meshless methods must use background cells to integrate a weak form
over the entire problem domain. The requirement of background cells for
integration makes the method not “truly” meshless. Another challenge lies in
enforcing boundary conditions accurately, especially so, as the degree of
freedom may not be adapted to the problem geometry. As pointed out by Atluri
and Zhu, in their work on meshless local Petrov-Galerkin methods, ensuring the
satisfaction of boundary conditions in meshless methods can be non-trivial,
potentially leading to inaccuracies in the solution; they attempt to mitigate
this by implementing a penalty function [1]. Additionally, the lack of a
structured mesh results in reduced convergence rates and increased
computational costs; this occurs particularly in problems with high gradients
or regions of varying solution behavior [4]. These limitations underscore the
importance of carefully assessing the suitability of meshless methods for
specific numerical analysis tasks and considering alternative approaches where
appropriate.

In amalgamating the desired advantages of meshless methods: accessibility,
simplicity, h-P convergence, and low complexity, while simultaneously
relinquishing the need for rigorous theoretical validation, we present the
Least Squares Discretization method (LSQD). This approach only assumes that
around each point, the solution can be decomposed on a local basis. No weak
formulation of the PDE, nor numerical analysis theory is needed to discretize
the problem. Instead, the PDE and boundary conditions are point-wise
evaluated, and the local expansions are connected through pointwise continuity
conditions. All of these conditions are expressed as linear equations
involving evaluations of the known local basis functions and the unknown
decomposition coefficients. They are regrouped in a rectangular linear system,
from which the vector of unknown coefficients is defined as the least squares
solution. The method is therefore entirely meshfree, requires minimal
mathematical and programming background, and scales with the number of grid
points.

Our computational exploration demonstrates how this method can be easily
implemented to approximate solutions of elliptic PDEs completed with arbitrary
boundary conditions enforced on non-trivial geometries, and using advanced
data structures. We observe the solution converging both as the spatial
resolution and the maximum polynomial order increase (h-P convergence), on
both uniform and arbitrary highly non-graded tree structures. Because our
method is meshless in nature, the underlying set of points can chosen as
virtually anything. In particular, the points do not have to conform to the
problem geometry or be aligned in any specific pattern. This facilitates the
construction and manipulation of the data structure and offers unbounded
solution adaptivity and error control strategies.

This manuscript is organized as follows. In section 2, we present the LSQD
method in the general context and discuss its practical implementation. The
method is illustrated with a simple one-dimensional example. In section 3 we
discuss the extension of the method to higher dimensions and propose an
implementation on non-graded quadtree grids. Convergence results are provided
in section 4 for a variety of grid and boundary configurations. We conclude in
section 5.

##  2 The LSQD Method

Figure 1: We represent the solution of a differential equation as a set of
local expansions, denoted by $u_{i}(x)$, which are centered at arbitrary
positions $x_{i=1..N}$. We build an overdetermined linear system by evaluating
the differential equation around each neighborhood and enforcing continuity
conditions between neighboring local representations. The coefficients
$\alpha_{i}^{q}$ are the least squares solution to the overdetermined system.

###  2.1 General Presentation

For the LSQD method, we begin by considering the following elliptic problem:

$\displaystyle au-\mu\triangle
u=f(\mathbf{x})\qquad\forall\mathbf{x}\in\Omega\subseteq\mathbb{R}^{d},\quad
f(\mathbf{x})\in\mathcal{C}^{0}(\Omega)$ (1)

with generalized (e.g. Dirichlet, Neumann, Robin, etc.)111For Dirichlet
boundary conditions: $\gamma=0$. For Neumann boundary conditions: $\beta=0$.
For Robin boundary conditions: $\beta\neq 0$ and $\gamma\neq 0$. boundary
conditions on the contour $\partial\Omega$, i.e.

$\displaystyle\beta u(\mathbf{x})+\gamma\nabla
u(\mathbf{x})\cdot\mathbf{n}=g(\mathbf{x}),\qquad\forall\mathbf{x}\in\partial\Omega,\quad
g(\mathbf{x})\in\mathcal{C}^{0}(\Omega),$ (2)

where $\mathbf{n}$ is the outward point normal to $\partial\Omega$. The
parameters of the elliptic problem, $a$ and $\mu$, are positive and the
solution is assumed to be continuously differentiable. We are interested in
constructing a discrete approximation of the solution $u(\mathbf{x})$, using a
set of points $(\mathbf{x}_{1},...,\mathbf{x}_{N})\in\Omega$. To each point
$\mathbf{x}_{i}$, we associate a set of neighbors $V_{i}$ and define the
vicinity $\Omega_{i}$ as the convex hull of this set. For each vicinity
$\Omega_{i}$, we approximate the solution using an expansion over a local
smooth basis $(\phi^{1}_{i}(\mathbf{x}),...,\phi^{Q}_{i}(\mathbf{x}))$

$\displaystyle
u(\mathbf{x})\approx\sum_{q=1}^{Q}\alpha_{i}^{q}\phi^{q}_{i}(\mathbf{x}),\qquad\forall\mathbf{x}\in\Omega_{i}.$
(3)

Thus, by estimating the coefficients $\alpha_{i}^{p}$, at each $\Omega_{i}$,
we hope to construct a family of local approximations to the solution
$u(\mathbf{x})$.

To find the coefficients $\alpha_{i}^{p}$ we evaluate the PDE (1), the
continuity conditions (for $\mathbf{u}$ and its first order derivatives), and
the boundary condition (2) for the above local representation. For each node
$\mathbf{x}_{i}$, we consider each neighbor pair
$(\mathbf{x}_{i},\mathbf{x}_{j})$ with $\mathbf{x}_{j}\in V_{i}$ and construct
the linear equations for $\alpha_{i}^{p}$ using the following procedure with
three main steps:

  1. 1.

Partial Differential Equation: we choose a point
$\mathbf{x}^{PDE}_{ij}\in\Omega_{i}$ at which we evaluate the PDE for the
expansion (3).

$\displaystyle
a\sum_{q=1}^{Q}\alpha_{i}^{q}\phi^{q}_{i}(\mathbf{x}^{PDE}_{ij})-\mu\sum_{q=1}^{Q}\alpha_{i}^{q}\triangle\phi^{q}_{i}(\mathbf{x}^{PDE}_{ij})=f(\mathbf{x}^{PDE}_{ij}).$
(4)

  2. 2.

Continuity conditions: we choose a point $\mathbf{x}^{C}_{ij}\in\Omega_{i}$ at
which we enforce that the expansions centered at $\mathbf{x}_{i}$ and
$\mathbf{x}_{j}$ are identical

$\displaystyle\sum_{q=1}^{Q}\alpha_{i}^{q}\phi^{q}_{i}(\mathbf{x}^{C}_{ij})=\sum_{q=1}^{Q}\alpha_{j}^{q}\phi^{q}_{j}(\mathbf{x}^{C}_{ij}),$
(5)

and that the expansions for all first-order derivatives $\partial_{d}$ are
also identical at that point

$\displaystyle\sum_{q=1}^{Q}\alpha_{i}^{q}\partial_{d}\phi^{q}_{i}(\mathbf{x}^{C}_{ij})=\sum_{q=1}^{Q}\alpha_{j}^{q}\partial_{d}\phi^{q}_{j}(\mathbf{x}^{C}_{ij}).$
(6)

  3. 3.

Boundary conditions: if $\mathbf{x_{i}}$ is close to the contour
$\partial\Omega$, we find a point $\mathbf{x}^{BC}_{i}\in\partial\Omega$, at
which we enforce the boundary condition for the local expansion (3).

$\displaystyle\beta\sum_{q=1}^{Q}\alpha_{i}^{q}\phi^{q}_{i}(\mathbf{x}^{BC}_{i})+\gamma\sum_{q=1}^{Q}\alpha_{i}^{q}\nabla\phi^{q}_{i}(\mathbf{x}^{BC}_{i})\cdot\mathbf{n}=g(\mathbf{x}^{BC}_{i}).$
(7)

Grouping all these equations, we form a rectangular system of width $Q\times
N$ and height $M$ (the total number of equations).

$\displaystyle\begin{bmatrix}\ldots\quad\ldots\quad\ldots\\\
\ldots\quad\ldots\quad\ldots\\\ \vdots\\\ \vdots\\\ \vdots\\\
\ldots\quad\ldots\quad\ldots\end{bmatrix}\begin{bmatrix}\alpha_{1}^{1}\\\
\alpha_{1}^{2}\\\ \vdots\\\ \alpha_{N}^{Q}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}b_{1}\\\ b_{2}\\\ \vdots\\\ \vdots\\\ \vdots\\\
b_{M}\end{bmatrix}\qquad\Leftrightarrow\qquad A\alpha=b.$ (8)

If the matrix $A$ is full rank, we construct the coefficient vector $\alpha$
as the least squares solution:

$\displaystyle\alpha=\underset{x}{\mathrm{argmin}}||Ax-b||_{2}\qquad\Leftrightarrow\qquad
A^{T}A\alpha=A^{T}b.$ (9)

If the matrix is not full-rank, we can consider the stabilized problem

$\displaystyle\left(A^{T}A+\epsilon I\right)\alpha=A^{T}b,$ (10)

which is well-known to be invertible for any $\epsilon>0$. Nevertheless, we
must be cognizant of two extremities: choosing a large $\epsilon$ will affect
the accuracy of the method, but choosing its value to be too small may leave
the system singular to machine precision. We address this topic in section
3.5.

###  2.2 Practical Implementation

####  2.2.1 Basis Selection

Formally, the method does not impose any restriction on the basis selection,
yet, it is well known that a naive choice could lead to singular least squares
systems, making it impossible to construct the solution (without
stabilization). For example, allow us to consider the standard interpolation
problem (i.e. choosing $\alpha=1$ and $\mu=0$ in (1), with appropriate
boundary conditions); from the Mairhuber–Curtis theorem [38], we know that
choosing global non-centered basis functions may result in a singular
interpolation matrix. For the rest of this study, at each point
$\mathbf{x}_{i}$ we will use polynomial functions, centered at the point
$\mathbf{x}_{i}$ and rescaled by the size of the neighborhood $V_{i}$. In
doing so, we aim to limit singularities in the final least squares system, as
well as provide a framework for achieving high-order convergence. We will
denote by $P$ the maximum polynomial order of the basis, and continue to
denote by $Q$, the total number of elements in the basis.

####  2.2.2 Neighborhood Construction

From the method description, we see that it is important to construct a full-
rank matrix A, while trying to keep it as compact and sparse as possible for
computational efficiency. The matrix structure and size are primarily dictated
by the neighborhoods configurations. It is therefore imperative to understand
how many neighbors are needed in each neighborhood, how to construct them, and
what properties they should satisfy. In particular, we identify three key
properties the neighborhood must satisfy

  1. i

System Size -  The total number of equations must be greater or equal to the
total number of basis functions (i.e. the number of unknown coefficients).

  2. ii

Connectivity -  The neighbors’ adjacency matrix must be invertible, or in
other words the neighboring graph must be irreducible. In practice, it means
that traveling along neighboring connection, you should be able to travel from
any point $\mathbf{x}_{i}$ to any other point $\mathbf{x}_{j}$.

  3. iii

Independence -  The neighborhood should contain points with at least $P+1$
unique coordinates in each directions.

The first condition is necessary for the matrix $A$ to be full rank. The
second condition ensures that all the neighborhoods are connected, and
therefore, that the computation of each basis’s coefficients depends on all
other coefficients. If this condition is not satisfied, there will be clusters
of nodes and, consequently, basis coefficients that are separated from each
other, violating the global nature of the elliptic equation (1). The third
condition reinforces linear independence between the coefficient equations. It
is constructed from the observation that if the problem is invariant in all
directions except one, the matrix A must still be full rank. In that case, the
number of linearly independent basis functions is simply $P+1$ and so we need
at least as many unique coordinates in the non-invariant direction.

We point out that the neighborhoods do not have to be symmetric, meaning that
if $\mathbf{x}_{j}$ is a neighbor of $\mathbf{x}_{i}$ (i.e. $\mathbf{x}_{j}\in
V_{i}$), $\mathbf{x}_{i}$ does not have to be a neighbor of $\mathbf{x}_{j}$.
They do not all have to be the same size, nor respect specific geometric
configurations. In practice though, it is natural to keep them local and
compact.

####  2.2.3 Evaluation Points

As mentioned in 2.1, the evaluation points $\mathbf{x}^{PDE}_{ij}$ and
$\mathbf{x}^{C}_{ij}$ for the PDE and continuity conditions must be chosen in
the domain $\Omega_{i}$ where the local $i^{th}$ expansion is formally
defined. Typically, the evaluation points are selected to be on the line
segment $[\mathbf{x}_{i},\mathbf{x}_{j}]$ as this ensures locality. We note
that choosing different evaluation points for the PDE and continuity
conditions ($\mathbf{x}^{PDE}_{ij}\neq\mathbf{x}^{C}_{ij}$) may reinforce
linear independence.

The boundary conditions are only evaluated when the point $\mathbf{x}_{i}$ is
close to the interface $\partial\Omega$. As well as other method parameters,
proximity of the interface must be defined in each context. For example, in
1D, we can place boundary points on the interface and evaluate the boundary
condition at these points only. In higher dimension, provided that we have a
smooth implicit representation (e.g. a level-set representation), we can
always detect whether $\Omega_{i}$ crosses $\partial\Omega$, and if so, find
at least one point $\mathbf{x}^{BC}_{i}$ on $\Omega_{i}\cap\partial\Omega$.

####  2.2.4 System Construction & Inversion

The main advantage of our method is that the matrix $A$ is straightforward to
construct. Once the neighborhoods have been constructed, we need to loop over
all neighboring pairs and evaluate the PDE, continuity, and boundary
conditions. The matrix assembly only requires evaluation of the basis function
at the evaluation points. No quadrature rule, finite difference formula, or
any other numerical analysis algorithms are required.

Assuming that the Least Squares problem is non-singular, the solution vector
can be computed using a preconditioned conjugate gradient method (since
$A^{T}A$ is symmetric), or using the LSQR method [43, 51] to avoid forming
$A^{T}A$. Forming $A^{T}A$ may allow us to construct a wider range of
preconditioners, and so it is, a priori, unclear which of these options will
be computationally more efficient. It should be noted that the matrices $A$,
$A^{T}$, and $A^{T}A$ are sparse and can be constructed in
$\mathcal{O}(NQ^{2})$ and $\mathcal{O}(NQ^{3})$ respectively. Therefore,
leaving aside the neighborhood construction, the total computational cost will
scale as the complexity of the linear solver, typically
$\mathcal{O}(NQ^{3}\sqrt{\kappa})$, with $\kappa$ being the condition number
of $A^{T}A$.

###  2.3 One-Dimensional Illustration

We consider the following Poisson problem

$\displaystyle-\triangle u$ $\displaystyle=$ $\displaystyle f(x)\qquad\forall
x\in\Omega=[0,1],$ (11) $\displaystyle u(x)$ $\displaystyle=$ $\displaystyle
g(x)\qquad\forall x\in\partial\Omega.$

In order to study the performance of the method we choose the exact solution
to be

$\displaystyle u_{\textrm{exact}}(x)=\sin(7x)+\cos(7x)$

and we set $f(x)$ and $g(x)$ accordingly. All local basis function sets are
chosen to be centered and rescaled polynomials, i.e. for each $x_{i}$ and over
each corresponding neighborhood $V_{i}$, we define

$\displaystyle\phi_{i}^{q}(x_{i})=\left(\frac{x-x_{i}}{|V_{i}|}\right)^{q-1},\qquad|V_{i}|=\max_{j\in
V_{i}}|x_{i}-x|,\qquad x\in V_{i},\quad q=\\{1,..,Q=P+1\\}.$ (12)

Each neighborhood $V_{i}$ is first created by choosing the node
$\mathbf{x}_{i}$ itself, its closest left and right neighbors. This ensures
that the neighborhoods overlap with at least another neighborhood in each
direction and therefore that the connectivity condition ii is satisfied.
Additionally, because the cardinality of the basis set is equal to $P+1$, we
need at least $P+1$ different neighbors in $V_{i}$ to satisfy both the system
size i and independence iii conditions. Therefore, if $P>2$, we add the
remaining $P-2$ closest points to the neighborhood. As a result, each
neighborhood will contain $\eta=\max(3,P+1)$ elements. For simplicity, we
choose the evaluation points to be
$\mathbf{x}_{ij}^{PDE}=\mathbf{x}_{ij}^{C}=\mathbf{x}_{j}$.

Once again, we emphasize the fact that advanced numerical analysis techniques
are not required to construct our system, nor to obtain our solution (e.g. no
Taylor expansions, convoluted integrals, or quadrature rules). Our
implementation contains around 100 lines of MATLAB code and the least squares
system is solved using the simple backslash solve command (i.e.
$\alpha=A^{T}A\backslash A^{T}b$). The backslash command for general square
matrices uses LU decomposition; however, for symmetric matrices (which is our
case) the Cholesky decomposition is used instead. The results can be seen in
Figure 3.

To study the h-P convergence of the method, we initialize the points
$x_{i=1..10}$ at random and then recursively add midpoints between adjacent
neighbors to double the local spatial resolution. We will refer to splits as
the number of subdivisions performed from the original points set, and point
out that the local resolution scales like $h^{-\text{splits}}$, where $h$ is
the original local spatial resolution. In parallel, we vary the polynomial
order using $P=2,3,4$. The error is computed by comparing the exact solution
and the local approximation at each point. We will focus on the maximum error
only.

Convergence results are compiled in Figure 3, where we see the error
diminishing both as the resolution and the polynomial order increase. On the
other hand, the condition number inflates with both the resolution and maximum
polynomial order; this phenomenon should be seen as a proxy for the
computational time needed to solve the least squares problem. This is as one
might expect and is the exchange we are willing to make in the name of
simplicity for this method.

With the system being symmetric, we can precondition it using a standard
incomplete Cholesky decomposition, which, for this example, drastically
reduces the condition number. Typically, the condition number of the
preconditioned system is the square root of the unpreconditioned one. As we
will see in section 4.6 the discontinuities between adjacent expansions can be
used to estimate.

Figure 2: 1D Method - solution representation for increasing spatial and
polynomial resolution. We start from a randomly generated set of points
$\mathbf{x}_{i=1..10}$ and then recursively split the grid by adding one grid
point exactly between each existing pair of direct neighbors. The exact
solution is depicted in dashed lines.

Figure 3: 1D Method - (a) h-P convergence analysis for Poisson’s equation with
Dirichlet boundary conditions. $L^{\infty}$ error and corresponding Estimated
Order of Accuracy (EOA). (b) Condition number of the least squares matrix
$A^{T}A$. (c) Condition number of the preconditioned least squares matrix,
using an incomplete Cholesky decomposition.

##  3 Extension to Higher Dimensions & Implementation on Quadtree Grids

Figure 4: Extension to higher dimensions: domain of interest $\Omega$, and
local vicinities $\Omega_{i},\Omega_{j}$. The neighborhoods $V_{i},V_{j}$ are
defined as the set of points contained in the corresponding vicinities.

Despite the one-dimensional example serving as an obligatory proof of concept
in terms of convergence, the full capabilities of our method, such as dealing
with highly adaptive data structure and arbitrary geometries and boundary
conditions, only manifest themselves in higher dimensions.

For our method, the extension to higher dimensions is straightforward from
both the theoretical and practical perspectives. As illustrated in Figure 4,
the domain of interest $\Omega$ is still discretized by a set of points
$\mathbf{x}_{i=1..N}$. The solution is then decomposed over a local basis
$\phi_{i}^{1}(\mathbf{x})..\phi_{i}^{Q}(\mathbf{x})$ around the set of points,
and the coefficients for these expansions are obtained by evaluating the
partial differential equations, boundary, and continuity conditions over local
neighborhoods, as underscored in section 2.1. The construction of the linear
system remains unambiguous and the properties of the final system are
unchanged ($A^{T}A$ is symmetric, positive definite).

The main difference lies in the neighborhood construction: both in their shape
and size. In 1D, finding close neighbors from a pre-sorted list is trivial. In
higher dimensions, for an efficient construction, the pre-sorting must be done
using more complex data structures, such as sorting trees or hash tables.

Also, because the derivative continuity condition will be evaluated in each
direction, the total number of equation written per neighboring pair will
differ, and thus the minimal neighborhood size. Namely, in $d$-dimensions, for
a neighborhood $V_{i}$ composed of $\eta$ neighbors (including
$\mathbf{x}_{i}$), we will evaluate

  * -

$\eta$ PDE equations,

  * -

$\eta-1$ continuity conditions,

  * -

$d(\eta-1)$ first-order derivative continuity conditions ($\eta-1$ for each
partial derivative),

so a total of $\eta\left(d+2\right)-d-1$ equations. Since we want the matrix
$A$ to be full rank, we will systematically pick the number of basis of local
basis function $Q$ such that

$\displaystyle Q\leq\eta\left(d+2\right)-d-1.$ (13)

Note that this strategy does not account for the boundary condition
evaluations. This is an intentional decision to derive a robust and simple
strategy: the total number of boundary condition evaluations is hard to
predict a priori, and these conditions will only increase the system size and
rank.

Since sorting structures will be needed to construct neighborhoods, we will
locate the points at the cell centers of a Quadtree to limit the number of
data structures and facilitate the presentation and implementation of our
method. Quadtree grids [47], a paradigmatic data structure for high-
performance computations (see e.g. [30, 7, 17, 23]). The points
$\mathbf{x}_{i}$ will be located at the tree’s cell centers. We will represent
the domain $\Omega$ through a level-set function and only approximate the
solution at the points located inside $\Omega$ (see Figure 6).

###  3.1 2D Basis Functions

Figure 5: LSQD Approximation of the function $u_{\textrm{exact}}=\cos(x+y)$
for increasing polynomial order on an arbitrary Quadtree grid.

We chose our basis functions as products of one-dimensional rescaled
polynomials. Namely around $\mathbf{x}_{i}=(x_{i},y_{i})$:

$\displaystyle\phi_{i}^{q}\left(\mathbf{x}\right)=\left(\frac{x-x_{i}}{|V^{x}_{i}|}\right)^{p_{x}}\left(\frac{y-y_{i}}{|V^{y}_{i}|}\right)^{p_{y}}\qquad\forall\mathbf{x}=(x,y)\in\Omega_{i},$
(14)

where the rescaling coefficients are

$\displaystyle\qquad|V_{i}^{x}|=\max_{\mathbf{x}\in
V_{i}}|x_{i}-x|,\qquad|V_{i}^{x}|=\max_{\mathbf{y}\in V_{i}}|y_{i}-y|.$ (15)

The rescaling helps keep the coefficients in $A$ – therefore of $A^{T}A$ – of
similar magnitude, consequently keeping the condition number of $A^{T}A$ low.
In this context, this is essential as we are dealing with very disparate
length scales and potentially high-order polynomial functions.

We consider the set of polynomial functions of degree $P$, meaning that
$p_{x}+p_{y}\leq P$ (see Figure 5). The total number of basis functions will
be

$\displaystyle Q=\frac{(P+1)(P+2)}{2},$ (16)

and so conditions (13) dictates that we will need to choose the neighborhood
size $\eta$ such that

$\displaystyle\eta\geq\frac{(P+1)(P+2)+6}{8}.$ (17)

Consider, for example, using a polynomial order $P=7$, from the above
condition we see that $\eta\geq 9$. The neighborhood will grow with the number
of basis functions and, thus, local unknown coefficients. However, this growth
will occur slowly, implying that the neighborhood will remain very compact;
this attribute is highly desirable, especially in higher dimensions, as
keeping “things” as local as possible implies keeping $A$ and $A^{T}A$ as
sparse as possible. The above inequality (17) is not enough to specify $\eta$
as the connectivity (ii) and independence (iii) conditions must also be
satisfied.

###  3.2 Tree-Based Neighborhood Construction

Figure 6: Illustration of neighborhood construction around a reference cell
(in blue) for polynomial basis order $P=4$. The process is carried out one
layer at a time until enough neighbors are identified to satisfy the system
size and independence conditions. When the reference cell lays next to the
interface, the neighbors are chosen within $\Omega$ (bottom row). The solution
is only constructed at the cell centers of the orange cells located inside
$\Omega$. The outside grey cells are only used to locate the interface while
evaluating the boundary condition.

Figure 6 illustrates the construction of the neighborhood $V_{i}$ for any
given point $\mathbf{x}_{i}$. We start by adding the point itself (blue) to
the neighborhood. Next, the algorithm recursively picks a direct neighbor in
each direction for each point in the neighborhood and store them in a
candidate set (pink). In this set, the candidates are sorted in
counterclockwise order starting to the right. In case one of the candidates
lies outside $\Omega$ (see bottom row of Figure 6), we move it to the ghost
set $G_{i}$ that will be used to find evaluation points for the boundary
condition.

On the first pass, we add the entire candidate set to the neighborhood to
ensure that the connectivity condition (ii) is met. On the subsequent passes,
while the other two conditions are not met, we go through the candidates set
and add a candidate $\mathbf{c}$ to the neighborhood if either

  1. (a)

The system size condition (i), reformulated as inequality (17), is not met.

  2. (b)

The independence condition (iii) is not met in at least a direction $d$ and
the coordinate $c_{d}$ of the candidate in that direction is not already
present in the neighborhood, i.e. $\forall\mathbf{v}\in V_{i}\quad v_{d}\neq
c_{d}$.

We point out that using non-uniform point distributions reduces repetition in
the coordinates, which eases the fulfillment of condition (iii) and ultimately
generates smaller and more compact neighborhoods. If the grid is refined
enough, the algorithm mentioned above will always converge, which allows us to
systematically construct neighborhoods. However, in case the grid is under-
resolved, there may not be enough neighbors available, which can lead to
infinite loops. To prevent this, we set a maximum number of passes (usually
$P+4$). In these degenerated situations, the grid is inconsistent with the
polynomial order and geometry, so the problem is completely under-resolved,
and the method will produce an inaccurate solution. In the next section, we
will explain how we use the neighborhood to evaluate the problem and form a
linear system.

###  3.3 Evaluation points

For simplicity, we keep the continuity and pde evaluation points identical
($\mathbf{x}^{PDE}_{ij}=\mathbf{x}^{C}_{ij}=\mathbf{x}_{j}$). For the
evaluation of the boundary condition, we go through the set of ghost neighbors
$G_{i}$ obtained through the neighborhood construction, and for each ghost
$\mathbf{g}$ neighbor we find a point
$\mathbf{x}^{BC}_{ig}\in[\mathbf{x}_{i},\mathbf{g}]$ lying on the interface
and evaluate the boundary condition there. Since we are using an implicit
representation for the interface (level-set), finding $\mathbf{x}^{BC}_{ig}$
boils down to a zero-finding problem on the segment
$[\mathbf{x}_{i},\mathbf{g}]$ which is easily solved using, for example, the
bisection method.

If $\mathbf{x}_{i}$ is located at the center of a cell at the edge of the
computational domain, the boundary condition is evaluated on the corresponding
edges.

###  3.4 Matrix Rescaling

In addition to the basis function rescaling, we rescale the rows of the matrix
$A$ by their respective norms. This operation should be seen as a re-weighting
of each linear equation, i.e. multiplying the entire system by a diagonal
matrix with non-zero diagonal entries. In doing so, we keep the magnitude of
the coefficients in $A$ as close to each other as possible to ensure that the
final condition number is as low as possible. Specifically, we start by
computing the weighting coefficients $w_{k}$ as the inverse of the row norms
of $A=(a_{kl})$

$\displaystyle w_{k}=\left(\sum_{l=1}^{Q\times
N}|a_{kl}|\right)^{-1}\qquad\forall k=1,2,...M$ (18)

and then rescale the entire system

$\displaystyle a_{kl}$ $\displaystyle=$
$\displaystyle{a_{kl}}{w_{k}}\qquad\forall k=1,2,...M,\quad l=1,2,..Q\times N$
(19) $\displaystyle b_{k}$ $\displaystyle=$
$\displaystyle{b_{k}}{w_{k}}\qquad\forall k=1,2,...M.$ (20)

Recall $M$ is the total number of equations in the overdetermined system.

###  3.5 Stabilization of the LSQ System

As discussed earlier, the final least squares system is guaranteed to be
invertible if a small diagonal perturbation $\epsilon$ is added. This addition
may not be necessary, especially as it essentially penalizes large
coefficients, thus affecting the accuracy of the method. In practice, this
stabilization should be seen as a means to keep the condition number of the
system bounded.

Adding the diagonal perturbation effectively shifts the spectrum of the
system. We know that the eigenvalues $\lambda_{k}$ of $A^{T}A$ are positive,
and so the condition number of $A^{T}A$ is the ratio of the largest and
smallest eigenvalues $\lambda_{\max}$ and $\lambda_{\min}$

$\displaystyle\kappa(A^{T}A)=\frac{\lambda_{\max}}{\lambda_{\min}}.$

If $\lambda_{k}$ is an eigenvalue of $A^{T}A$, $\lambda_{k}+\epsilon$ is an
eigenvalue of the stabilized problem $A^{T}A+\epsilon I$, and so the condition
number of the stabilized problem is

$\displaystyle\kappa(A^{T}A+\epsilon
I)=\frac{\lambda_{\max}+\epsilon}{\lambda_{\min}+\epsilon}.$

If we set a practical maximum condition number $\kappa_{M}$ for the stabilized
system, and seek the minimal perturbation needed $\epsilon_{\min}$ to remain
below this limit

$\displaystyle\kappa_{\max}=\frac{\lambda_{\max}+\epsilon_{\min}}{\lambda_{\min}+\epsilon_{\min}},$

we find that

$\displaystyle\epsilon_{\min}=\frac{\lambda_{\max}-\kappa_{\max}\lambda_{\min}}{\kappa_{\max}-1}.$

For the current implementation we set $\kappa_{\max}=10^{40}$, and estimate
the largest and smallest eigenvalues of $A^{T}A$ using Gershgorin circles
theorem

$\displaystyle\lambda_{\max}$ $\displaystyle=$
$\displaystyle\max_{k=1..M}c_{kk}+\sum_{l\neq k}|c_{kl}|,$
$\displaystyle\lambda_{\min}$ $\displaystyle=$
$\displaystyle\max\left(0,\min_{k=1..M}c_{kk}+\sum_{l\neq k}|c_{kl}|\right),$

where the $c_{kl}$ are the coefficients of the matrix $A^{T}A$. In practice,
we measured $\lambda_{\max}\approx 10^{2}$, $\lambda_{\min}=0$, leading to
$\epsilon\approx 10^{-38}$. We note that $\kappa_{\max}$ is the maximum
allowed condition number for the preconditioned system. In practice, we
preconditioned the system with an incomplete Cholesky precondition, and so the
final condition number will be much lower. Based on our one-dimensional
observations, we expect the preconditioned condition number to be bounded
above by $\approx 10^{20}$. Again, this is an upper limit that should be seen
as the absolute worst-case scenario. For all cases presented in the following
section, we find the least squares system to be invertible and construct the
numerical solution using an iterative Conjugate Gradient method, typically
taking $10^{2}-10^{4}$ iterations to converge to residuals
$||(A^{T}A+\epsilon_{\min}I)\alpha-A^{T}b||_{\infty}$ smaller than $10^{-26}$.

Figure 7: Configurations for the convergence analysis: (a) Illustration of
arbitrary tree structures. The grid is randomly generated (first grid) which
we call Split 0, then recursively refined (left to right - Split 1 to Split
3). (b) Grids and domains used - For each boundary condition type (Dirichlet,
Neumann, Mixed Dirichlet Robin), enforced on a square (top rows) or the
octofoil domain (bottom rows), we will consider both uniform (row 1 and 3) and
adaptive grids (rows 2 and 4). This representation will be used to label all
convergence figures hereafter

.

Table 1: Convergence summary - Average estimated order of convergence for various boundary conditions, configurations, and polynomial orders $P=2,..,5$. | Dirichlet | Neumann | Dir.+Robin | Vary $\frac{\alpha}{\mu}$ | Degenerate | Overall Avg
---|---|---|---|---|---|---
$\mathbf{P=2}$ | 1.86 | 2.17 | 2.02 | 1.63 | 1.46 | 1.83
$\mathbf{P=3}$ | 1.96 | 2.20 | 2.07 | 2.16 | 1.91 | 2.06
$\mathbf{P=4}$ | 3.45 | 3.30 | 3.46 | 3.51 | 3.45 | 3.43
$\mathbf{P=5}$ | 3.66 | 4.53 | 3.89 | 4.46 | 3.60 | 4.03

##  4 Results

Presented here are convergence results for various domain types, levels of
mesh refinement, boundary condition types, and numerous formulations of
equation (1) (by changing the coefficients). Furthermore, we present built-in
error estimate. Finally, we discuss the implications of Neumann boundary
conditions with a computational time analysis.

###  4.1 Problem Setup

We consider the proposed elliptic problem (1) with varying $\alpha,\mu$ values
over two domains. The first domain is the unit square
$\Omega=[-1,1]\times[-1,1]$. The second domain is the interior of an octofoil
centered at $(0,0)$ and defined by the zero contour of the following level set
equation (see Figure 7)

$\displaystyle\phi(x,y)=-\Biggl{[}\frac{4}{5}+\frac{4}{25}\sin\biggl{(}8\arctan\left(\frac{y}{x}\right)\biggr{)}\Biggr{]}+\sqrt{x^{2}+y^{2}}.$
(21)

For both of these domains, we use the boundary condition formulation (2) and
adjust the weights $\beta$ and $\gamma$ to achieve the desired condition
types. We also consider two types of mesh refinement: uniform and randomly
generated (see Figure 7). Each of the representations will be used for the
convergence analysis in section 4.

We choose the exact solution to be $u_{\textrm{exact}}(x,y)=e^{x+y}$ and we
set the forcing term $f(\mathbf{x})$ and the boundary condition
$g(\mathbf{x})$ accordingly. For each of the representations, we vary the
polynomial order of our basis from 2 to 5, and recursively split all grid
cells to study the h-P convergence of the $L^{\infty}$-errors (see Figure
7(a)).
The local $L^{1}$ and global $L^{\infty}$-error are defined respectively as

$\displaystyle||e_{u}(\mathbf{x_{i}})||_{1}=\left\lvert
u_{\textrm{exact}}(\mathbf{x}_{i})-\sum_{q=1}^{Q}\alpha_{i}^{q}\phi_{i}^{q}(\mathbf{x}_{i})\right\rvert.$
(22) $\displaystyle||e_{u}||_{\infty}=\max_{i=1..N}\left\lvert
u_{\textrm{exact}}(\mathbf{x}_{i})-\sum_{q=1}^{Q}\alpha_{i}^{q}\phi_{i}^{q}(\mathbf{x}_{i})\right\rvert.$
(23)

The Estimated Order of Convergence (EOC) will be computed from the global
error, using power fitting curves.

Figure 8: Convergence of the local error: Octofoil domain with Dirichlet
boundary conditions and an adaptive grid. We see the error decaying to zero
both as the polynomial order increases or the grid is refined.

###  4.2 Classical Boundary Conditions: Dirichlet, Neumann, and Robin

Consider $\alpha=\mu=1$ for all twelve configurations presented in Figure 7.
The local error convergence for the octofoil domain with a Dirichlet Boundary
condition ($\beta=1$), using an adaptive refinement is depicted on Figure 8.
The h-P convergence of the global error for each set of boundary conditions is
represented in Figure 9. For all twelve configurations, we see the error
decaying to zero as the grid resolution and the polynomial order increase. The
estimated order of convergence (see Table 1 and Figure 9) are fairly constant
across all examples, ranging from approximately $2$ (P=2) to $4$ (P=5).

We point out that the octofoil examples on the adaptive grid, and even more so
with the mixed boundary condition case (labeled ”Dir. + Robin” in Figure 7),
are particularly challenging to deal with traditional discretization
techniques. Again, non-trivial boundary conditions are enforced on complex
geometries, and the grid is randomly generated and not conforming to the
interface. In fact, for these last two reasons, the traditional Finite Element
Method cannot be employed here, as they require conforming mesh, and when
implemented on quadtree grids, requires graded tree (i.e. the size ratio
between neighboring cells cannot exceed 2). There are in theory no
obstructions from using the Finite Differences of Volumes methods, but the
constructions of the discrete differential operators and fluxes rapidly
becomes a practical nightmare, especially in the high-order context. The modal
Discontinuous Galerkin method allows for a systematic construction of the
discrete system, where volume and contour integrals must be performed on the
fractions of the computational cells and faces contained in the domain of
interest, which can become very involved when the mesh is not adapted to the
interface. Again, LSQD requires no convoluted discretization nor quadrature
formula. Even in the presence of complex data structure or geometries, simple
basis functions evaluations at pre-specified locations are enough to yield
above second-order convergence.

Figure 9: Convergence of $L^{\infty}$ global error for several configurations.
Columns denote boundary condition type. Rows correspond to (1) a uniform grid
on a square domain, (2) adaptive square, (3) uniform octofoil, and (4)
adaptive octofoil (top to bottom). See Figure 7 for mesh examples.

###  4.3 Degenerate Cases

Here, we investigate degenerate cases where the PDE or boundary conditions may
lead to a non-coercive weak formulation and thus potentially ill-posed
problems, non-positive or even non-invertible systems. The LSQD method obtains
the polynomial coefficient from an SPD system and does not rely on the weak
formulation. Similarly to the Least Square Finite Element Method (LSFEM [25,
24]), the LSQD method constructs the discrete system by evaluating the PDE and
boundary conditions everywhere (formation of the matrix A), and then defining
the Galerkin matrix taking the scalar product of these evaluations ($A^{T}A$).
As a result, it is not affected by the potential ill-posedness of the weak
form, and the resulting linear system is guaranteed to be invertible.

For the first case, we considered equation 1 with $a=1,\mu=1$ and imposed a
split boundary conditions of Dirichlet and Robin. The degeneration comes from
the Robin boundary condition, for which we set $\beta=1,\gamma=-1$. For the
second case, we keep $a=1$, but set $\mu=-100$ (i.e. Helmholtz’s equation),
and enforce a Dirichlet boundary condition. The convergence results are
depicted in Figure 10. The method achieves h-P convergence for both of the
degenerate cases and all grid configurations. The numerical error dips rapidly
below reasonable error tolerances, and the EOC behaves normally, increasing
with the polynomial order and split number. We note that results using
exclusively traditional Robin boundary conditions ($\beta=1$, $\gamma=1$) on
the entire domain achieved the similar results.

Figure 10: Two instances of a degenerate case: (a) Mixed Robin and Dirichlet
boundary conditions, (b) Dirichlet boundary condition for a highly negative
value of $\mu$.

###  4.4 Sensitivity to Problem Parameters

Despite consistent patterns of convergence for various boundary conditions
tested here (i.e. dirichlet, dirichlet + robin, degenerate, etc.), Neumann
boundary conditions seem to perform poorly – especially on adaptive grids. In
fact, in Figure 9, we observe errors of magnitude upwards of $10^{1}\text{ and
even }10^{2}$.

From (1), we note that when $\mu$ is significantly smaller than $a$, the
problem resembles an interpolation problem, which is relatively
straightforward to solve. On the other hand, if $a$ is considerably smaller
than $\mu$, the equation resembles a – potentially – ill-posed Poisson
problem, e.g. failing the solution uniqueness criterion. In fact, the well-
posedness of the Poisson problem under Neumann boundary conditions hinges on
factors such as the smoothness of the domain and boundary, the characteristics
of the Neumann boundary data $g$, and how well the boundary conditions align
with the given problem. When the Neumann boundary condition fails to offer
sufficient information to determine $u$ uniquely, it can result in non-
uniqueness or instability within the solution, rendering the problem ill-
posed. Therefore, by setting $a=\mu=1$, the well-posedness of the problem
could have impacted the convergence/accuracy in the Neumann boundary condition
setting. To investigate this, we select four different values of
$\frac{a}{\mu}$. The expectation is that as we increase the ratio, the problem
should become easier to solve.

Figure 11: Neumann boundary condition on the octofoil domain: impact of
$\frac{\alpha}{\mu}$ ratios. (a) Global error trend for increasing split
number for a uniform grid for each polynomial basis. (b) Impact of
computational time on global error for both adaptive and uniform grid types.

Presented in Figure 11 (a) are convergence results for Neumann boundary
conditions imposed on the octofoil interface with a uniform grid for
$\frac{a}{\mu}=1,10,100,1000$. We observe that for each polynomial basis
order, $P$, the global error maintains a similar trend among the four tested
ratios $\frac{a}{\mu}$, respectively. This behavior corroborates the congruity
of the Estimated Order of Convergence (EOC) between various $\frac{a}{\mu}$
given order $P$. Clearly, the method performs as expected with the Neumann
constraint across the various ratios; nonetheless, it is noteworthy that
larger values of $\frac{a}{\mu}$ allow for a significant decrease in the
global error ($\approx 2$ orders of magnitude). In fact, when
$\frac{a}{\mu}\geq 100$, the method demonstrates accuracy comparable to, or
surpassing, that which was achieved with Dirichlet boundary conditions. Hence,
with Neumann boundary conditions, smaller values of $\frac{a}{\mu}$ can
compromise the effectiveness of the method.

###  4.5 Computational Time

Previously, it was established that lower $\frac{\alpha}{\mu}$ ratio values
associated with a Neumann boundary condition create a more challenging
computational problem (4.4). For this reason, we showcase the performance and
efficiency of the LSQD method by exploring a variety of $\frac{\alpha}{\mu}$
ratios. Figure 11 (b) illustrates the maximum error obtained by our solver in
comparison to the time it took to run in serial on a Dell precision 3240. It
is apparent that, in general, achieving a lower maximum error requires a
longer runtime and higher polynomial order. However, we also see that as the
ratio $\frac{\alpha}{\mu}$ decreases, not only does it take more time to solve
the system, but the order of accuracy declines.

We also note that the slope of the error versus computational times increases
with the polynomial order, suggesting that in the asymptotic limit of highly
resolved simulations (small error or large CPU time), high-polynomial
computations will always be the most efficient. Though, for the presented
results this asymptotic regime is never reached. Therefore, in practice, it is
imperative to decide which variable (minimizing error or computational cost)
takes precedence in a given setting. Presuming that the modeled system has a
tight tolerance for error, while computational resources are not a primary
concern, the conspicuous objective is achieving a low global error. Figure
11(b) hints at the necessary setup for attaining specific accuracy levels; the
higher the polynomial basis order, the higher the accuracy. However, if the
priority shifts towards saving computational time and/or memory or simply
capturing the general behavior of the solution, opting for a moderately
compact setup and a lower-order polynomial basis might suffice.

###  4.6 Error Estimator

We propose an alternative way to compute local error estimate as:

$\displaystyle||e_{u}(\mathbf{\widetilde{x}})||_{Est}=\max\left\lvert\sum_{q=1}^{Q}\alpha_{l_{1}}^{q}\phi_{l_{1}}^{q}(\mathbf{\widetilde{x}})-\sum_{q=1}^{Q}\alpha_{l_{2}}^{q}\phi_{l_{2}}^{q}(\mathbf{\widetilde{x}})\right\rvert,\quad
l_{1},l_{2}\in L=(\text{neighboring cells}).$ (24)

Obtaining analytical solutions to PDEs is only possible in the simplest cases,
so it should be expected that the exact solution to the PDE (1) may not be
available. This error estimator is built into the algorithm, as evaluating
accuracy predicated on traditional methods tends to be infeasible. For a
comparison between the theoretical prediction and practical results, we
restrict our analysis to the Dirichlet case.

The function of our node-based error estimator is to compare local
approximations from neighboring cells. All corners of the cells located inside
the domain become node locations $\mathbf{\widetilde{x}_{i}}$. For each node,
the touching cells are used to calculate the error. In one dimension, we can
picture the error as the discontinuity between two local expansions around
that node (see Figure 3). To visualize this process in two dimensional space,
refer once again to Figure 12(a). If we want to compute the estimate error
around the reference node (in blue), we use the local approximation functions
constructed for the closest neighbor cells (black), and evaluate at the
considered reference node $\mathbf{\widetilde{x}_{i}}$ (note that the maximum
amount of neighboring cells is four: left, right, top, and bottom). Finally,
the error estimate is computed as the maximum difference between the local
values obtained.

Focusing, in particular, on Figure 12(b-c), we are able to see a comparison
between simulation error and estimated error for multiple polynomial orders
and for increasing grid refinement. These results are encapsulated,
numerically, in 12(b), where calculations (on legend) generalize that the
predicted EOC is consistently higher than the actual EOC. However, when taking
a closer look at the distribution of the colored (simulated) and dashed lines
(predicted), the predicted global error is not free of variance, meaning that
neither the simulated, nor the predicted global error is consistently
superior. It becomes evident (see fig. 12(c)) that the main differences
between our actual and estimated error are transpiring at the boundary of our
octofoil domain. Nonetheless, the most stark differences dissipate as the grid
is refined (this is referred to as ‘Splits’). Overall, we find that this cheap
built-in estimator offers meaningful global error and convergence estimates.

Figure 12: Error Estimator Analysis on adaptive grids: (a) Three distinct node
locations and what is considered their respective surrounding/neighboring
cells. (b) Comparison between computational order of accuracy (colored lines)
and predicted order of accuracy (dashed lines). (c) Illustration of
computational and predicted error for increasing polynomial order and grid
refinement.

##  5 Discussion / Conclusion

In this paper, we proposed a new paradigm for solving elliptic PDEs, one that
would retain all of the advantages of meshless methods, while remaining
accessible in theory and in practice. Least Squares Discretization Method
(LSQD) does not require any quadarature rules, weak formulations, body-fitted
grids, nor does it necessitate advanced mathematical theory. Nevertheless, it
constructs h-P converging numerical solutions on non-trivial data structures
with reasonable computational complexity. LSQD method easily handles complex
geometries and is virtually impervious to ill-posed problems.

Although we began with a one-dimensional proof of concept, we demonstrated the
feasibility in extending our method to a higher dimensional space, and its
applicability to a broad range of problems, boundary conditions, and
computational parameters (e.g. grid resolution, polynomial orders). For all
considered cases, we consistently observe h-P convergence, with Estimated
Order of Convergence (EOC) ranging from 1.83 to 4.53 (see Table 1), and in
practice the global error falls rapidly below any reasonable error tolerance.
Moreover, when investigating degenerate cases (i.e. where the problem is
potentially ill-posed as a consequence of the PDE coefficients or boundary
conditions), the method observes kindred results (see 4.3), exemplifying the
robustness of the LSQD method.

Finally, we proposed an a posteriori node-based error estimator to unveil the
local solution agreement. From Fig 12 (c), we see that as the polynomial order
increases and the grid is refined, the estimated and measured/experimental
errors display the same convergence pattern and similar spatial localization.
This estimator cannot be used as a lower bound for actual error as seen in
12(c), but provides a general blueprint for where error is to be expected and,
thus, can be used for optimal mesh refinement. Coupled with a spatial
adaptivity capability, this error estimator could be used as an error
controller.

From our computational time analysis (see 4.5), we observed steeper efficiency
curves for higher polynomial orders, suggesting that using high-order
representations might be advantageous, at least for very large systems. In
practice, for small test problems like the ones considered here, these
benefits may not manifest, and low polynomial orders seem to be the best
strategy to get an accurate solution rapidly. In the examples presented here,
computational times peaked around $10^{3}$ seconds on a single core, or
$\approx 17$ minutes for the largest systems ($\approx 1.3\times 10^{6}$
degrees of freedom for 4 splits and P=5). Even though we are willing to
sacrifice computational time for ease of implementation, these values are
significant and suggest that the method must be optimized to be of full
practical interest.

In this case, we have two avenues for optimization, either we can prioritize
reducing the CPU time or the numerical error. Since most of the computational
time is dedicated to solving the least square system, accelerating the method
would entail speeding up the linear solver and ultimately reducing the
condition number. The method offers a number of free parameters and tuning
knobs, such as the positions of the evaluation points, the number of
neighbors, the basis functions selection, or the possibility to weight
equations differently in the overdetermined system. We believe the method can
be optimized by investigating the impact of these parameters on the accuracy
of the solution and conditioning of the system.

Despite our presentation and analysis being restricted to low-dimensional
examples and an un-optimized approach, we see that the LSQD method is a robust
and versatile tool for computational scientists. Its minimal theoretical and
programming barriers render it accessible to a broad scientific community, and
we anticipate that further research endeavors could propel it into a
formidable simulation paradigm.

##  6 Acknowledgments

The authors thank Arnold Kim for valuable discussions. This material is based
upon work supported by the National Science Foundation under Grant No.
DMS-1840265.

## Credit author statement

Anna Kucherova: Investigation, Formal Analysis, Software, Visualization,
Writing - Original Draft, Writing - Review $\&$ Editing. Gbocho Terasaki:
Formal Analysis, Investigation, Software, Visualization, Writing - Original
Draft, Writing - Review $\&$ Editing. Selma Strango: Preliminary
Investigation, Formal Analysis, Software, Writing - Review $\&$ Editing.
Maxime Theillard: Conceptualization, Supervision, Formal Analysis,
Methodology, Software, Project Administration, Writing - Original Draft,
Writing - Review $\&$ Editing.

## References

  * [1] S. N. Atluri and T. Zhu.  A new meshless local petrov-galerkin (mlpg) approach in computational mechanics.  Computational Mechanics, 22(2):117–127, Aug 1998.
  * [2] S. N. Atluri and T.-L. Zhu.  The meshless local petrov-galerkin (mlpg) approach for solving problems in elasto-statics.  Computational Mechanics, 25(2):169–179, Mar 2000.
  * [3] I. Babuska and J. M. Melenk.  The partition of unity method.  International Journal for Numerical Methods in Engineering, 40(4):727–758, 1997.
  * [4] T. Belytschko, Y. Krongauz, D. Organ, M. Fleming, and P. Krysl.  Meshless methods: An overview and recent developments.  Computer Methods in Applied Mechanics and Engineering, 139(1):3 – 47, 1996.
  * [5] T. Belytschko, Y. Y. Lu, and L. Gu.  Element-free galerkin methods.  International Journal for Numerical Methods in Engineering, 37(2):229–256, 1994.
  * [6] T. Belytschko, D. Organ, and C. Gerlach.  Element-free galerkin methods for dynamic fracture in concrete.  Computer Methods in Applied Mechanics and Engineering, 187(3):385 – 399, 2000.
  * [7] M. Blomquist, S. R. West, A. L. Binswanger, and M. Theillard.  Stable nodal projection method on octree grids.  Journal of Computational Physics, 499:112695, 2024.
  * [8] K. C. Chung.  A generalized finite-difference method for heat transfer problems of irregular geometries.  Numerical Heat Transfer, 4(3):345–357, 1981.
  * [9] P. C. Curtis.  $n$-parameter families and best approximation.  Pacific J. Math., 9(4):1013–1027, 1959.
  * [10] Q. Du, M. D. Gunzburger, and L. Ju.  Voronoi-based finite volume methods, optimal voronoi meshes, and pdes on the sphere.  Computer Methods in Applied Mechanics and Engineering, 192(35):3933 – 3957, 2003.
  * [11] Y. Duan.  A note on the meshless method using radial basis functions.  Computers & Mathematics with Applications, 55(1):66–75, 2008.
  * [12] C. A. Duarte and J. T. Oden.  H-p clouds—an h-p meshless method.  Numerical Methods for Partial Differential Equations, 12(6):673–705, 1996.
  * [13] A. Eirís, L. Ramírez, I. Couceiro, J. Fernández-Fidalgo, J. París, and X. Nogueira.  Mls-sph-ale: A review of meshless-fv methods and a unifying formulation for particle discretizations.  Archives of Computational Methods in Engineering, 30(8):4959–4981, 2023.
  * [14] B. Fornberg and N. Flyer.  Solving pdes with radial basis functions.  Acta Numerica, 24:215–258, 2015.
  * [15] R. Franke.  Scattered data interpolation: tests of some methods.  Mathematics of Computation, 38:181–200, 1982.
  * [16] E. J. Fuselier and G. B. Wright.  A high-order kernel method for diffusion and reaction-diffusion equations on surfaces.  Journal of Scientific Computing, 56:535–565, 2013.
  * [17] C. Fylling, J. Tamayo, A. Gopinath, and M. Theillard.  Multi-population dissolution in confined active fluids.  Soft Matter, 20:1392–1409, 2024.
  * [18] R. A. Gingold and J. J. Monaghan.  Smoothed particle hydrodynamics: theory and application to non-spherical stars.  Monthly Notices of the Royal Astronomical Society, 181(3):375–389, 12 1977.
  * [19] D. González, E. Cueto, and M. Doblaré.  Towards an Isogeometric Meshless Natural Element Method.  Springer Netherlands, Dordrecht, 2009.
  * [20] A. Guittet, M. Theillard, and F. Gibou.  A stable projection method for the incompressible navier-stokes equations on arbitrary geometries and adaptive quad/octrees.  J. Comput. Phys., 292:215–238, 2015.
  * [21] R. L. Hardy.  Multiquadric equations of topography and other irregular surfaces.  Journal of Geophysical Research (1896-1977), 76(8):1905–1915, 1971.
  * [22] B. He, J. Li, M. Chen, and Y. Huang.  A leap-frog meshless method with radial basis functions for simulating electromagnetic wave splitter and rotator.  Engineering Analysis with Boundary Elements, 118:225–242, 2020.
  * [23] A. A. Heydari, S. S. Sindi, and M. Theillard.  Conservative finite volume method on deforming geometries: The case of protein aggregation in dividing yeast cells.  Journal of Computational Physics, 448:110755, 2022.
  * [24] B.-n. Jiang.  Basis of LSFEM, pages 47–79.  Springer Berlin Heidelberg, Berlin, Heidelberg, 1998.
  * [25] B.-N. Jiang and L. A. Povinelli.  Least-squares finite element method for fluid dynamics.  Computer Methods in Applied Mechanics and Engineering, 81(1):13–37, 1990.
  * [26] G. Joldes, G. Bourantas, B. Zwick, H. Chowdhury, A. Wittek, S. Agrawal, K. Mountris, D. Hyde, S. K. Warfield, and K. Miller.  Suite of meshless algorithms for accurate computation of soft tissue deformation for surgical simulation.  Medical Image Analysis, 56:152–171, 2019.
  * [27] E. Kansa.  Multiquadrics—a scattered data approximation scheme with applications to computational fluid-dynamics—ii solutions to parabolic, hyperbolic and elliptic partial differential equations.  Computers Mathematics with Applications, 19(8):147 – 161, 1990.
  * [28] E. Kansa and Y. Hon.  Circumventing the ill-conditioning problem with multiquadric radial basis functions: Applications to elliptic partial differential equations.  Computers & Mathematics with Applications, 39(7):123–137, 2000.
  * [29] I. Kovacevic, A. Poredos, and B. Sarler.  Solving the stefan problem with the radial basis function collocation method.  Numerical Heat Transfer: Part B: Fundamentals, 44:575–599, 12 2003.
  * [30] A. Kucherova, S. Strango, S. Sukenik, and M. Theillard.  Computational modeling of protein conformational changes - application to the opening sars-cov-2 spike.  Journal of Computational Physics, 444:110591, 2021.
  * [31] M. Kumar, A. K. Jha, Y. Bhagoria, and P. Gupta.  A review to explore different meshless methods in various structural problems.  IOP Conference Series: Materials Science and Engineering, 1116(1):012119, apr 2021.
  * [32] S. Kumar, K. Danas, and D. M. Kochmann.  Enhanced local maximum-entropy approximation for stable meshfree simulations.  Computer Methods in Applied Mechanics and Engineering, 344:858–886, 2019.
  * [33] T. Liszka and J. Orkisz.  The finite difference method at arbitrary irregular grids and its application in applied mechanics.  Computers & Structures, 11(1):83 – 95, 1980.  Special Issue-Computational Methods in Nonlinear Mechanics.
  * [34] W. K. Liu, S. Jun, and Y. F. Zhang.  Reproducing kernel particle methods.  International Journal for Numerical Methods in Fluids, 20(8‐9):1081–1106, 1995.
  * [35] L. B. Lucy.  A numerical approach to the testing of the fission hypothesis.  82:1013–1024, 1977.
  * [36] N. Mai-Duy and T. Tran-Cong.  Numerical solution of navier–stokes equations using multiquadric radial basis function networks.  International Journal for Numerical Methods in Fluids, 37(1):65–86, 2001.
  * [37] N. Mai-Duy and T. Tran-Cong.  Mesh-free radial basis function network methods with domain decomposition for approximation of functions and numerical solution of poisson’s equations.  Engineering Analysis with Boundary Elements, 26(2):133 – 156, 2002.
  * [38] J. C. Mairhuber.  On haar’s theorem concerning chebychev approximation problems having unique solutions.  Proceedings of the American Mathematical Society, 7(4):609–615, 1956.
  * [39] F. Mazhar, A. Javed, J. T. Xing, A. Shahzad, M. Mansoor, A. Maqsood, S. I. A. Shah, and K. Asim.  On the meshfree particle methods for fluid-structure interaction problems.  Engineering Analysis with Boundary Elements, 124:14–40, 2021.
  * [40] C. A. Micchelli.  Interpolation of scattered data: Distance matrices and conditionally positive definite functions.  Constructive Approximation, 2(1):11–22, 1986.
  * [41] V. P. Nguyen, T. Rabczuk, S. Bordas, and M. Duflot.  Meshless methods: A review and computer implementation aspects.  Mathematics and Computers in Simulation, 79(3):763–813, 2008.
  * [42] E. Onate, S. Idelsohn, O. C. Zienkiwicz, and R. L. TAYLOR.  A finite point method in computational mechanics. applications to convective transport and fluid flow.  International Journal for Numerical Methods in Engineering, 39(22):3839–3866, 1996.
  * [43] C. C. Paige and M. A. Saunders.  Lsqr: An algorithm for sparse linear equations and sparse least squares.  ACM Trans. Math. Softw., 8(1):43–71, Mar. 1982.
  * [44] V. G. Patel and N. V. Rachchh.  Meshless method – review on recent developments.  Materials Today: Proceedings, 26:1598–1603, 2020.  10th International Conference of Materials Processing and Characterization.
  * [45] C. Piret.  The orthogonal gradients method: A radial basis functions method for solving partial differential equations on arbitrary surfaces.  Journal of Computational Physics, 231(14):4662 – 4675, 2012.
  * [46] M. J. Powell.  Five lectures on radial basis functions.  Technical report, Informatics and Mathematical Modelling, Technical University of Denmark, DTU, Richard Petersens Plads, Building 321, DK-2800 Kgs. Lyngby, 2005.
  * [47] H. Samet.  An overview of quadtrees, octrees, and related hierarchical data structures.  1988\.
  * [48] B. Sarler.  A radial basis function collocation approach in computational fluid dynamics.  Computer Modeling in Engineering & Sciences, 7(2):185–194, 2005.
  * [49] B. Sarler, J. Perko, and C.-S. Chen.  Radial basis function collocation method solution of natural convection in porous media.  International Journal of Numerical Methods for Heat & Fluid Flow, 14:187–212, 03 2004.
  * [50] B. Sarler and R. Vertnik.  Meshfree explicit local radial basis function collocation method for diffusion problems.  Computers Mathematics with Applications, 51(8):1269 – 1282, 2006.  Radial Basis Functions and Related Multivariate Meshfree Approximation Methods: Theory and Applications.
  * [51] M. A. Saunders.  Solution of sparse rectangular systems using lsqr and craig.  BIT Numerical Mathematics, 35(4):588–604, 1995.
  * [52] C. Shu, H. Ding, H. Chen, and T. Wang.  An upwind local rbf-dq method for simulation of inviscid compressible flows.  Computer Methods in Applied Mechanics and Engineering, 194(18):2001 – 2017, 2005.
  * [53] S. Soleimani, M. Jalaal, H. Bararnia, E. Ghasemi, D. Ganji, and F. Mohammadi.  Local rbf-dq method for two-dimensional transient heat conduction problems.  International Communications in Heat and Mass Transfer, 37(9):1411 – 1418, 2010.
  * [54] N. Sukumar, B. Moran, and T. Belytschko.  The natural element method in solid mechanics.  International Journal for Numerical Methods in Engineering, 43(5):839–887, 1998.
  * [55] G. B. Wright, N. Flyer, and D. A. Yuen.  A hybrid radial basis function–pseudospectral method for thermal convection in a 3-d spherical shell.  Geochemistry, Geophysics, Geosystems, 11(7), 2010.
  * [56] J. Yvonnet, D. Ryckelynck, P. Lorong, and F. Chinesta.  A new extension of the natural element method for non-convex and discontinuous problems: the constrained natural element method (c-nem).  International Journal for Numerical Methods in Engineering, 60(8):1451–1474, 2004.
  * [57] T. Zhu, J.-D. Zhang, and S. N. Atluri.  A local boundary integral equation (lbie) method in computational mechanics, and a meshless discretization approach.  Computational Mechanics, 21(3):223–235, Apr 1998.

@articleNGAI2011559, title = ”The application of data mining techniques in
financial fraud detection: A classification framework and an academic review
of literature”, journal = ”Decision Support Systems”, volume = ”50”, number =
”3”, pages = ”559 - 569”, year = ”2011”, note = ”On quantitative methods for
detection of financial fraud”, issn = ”0167-9236”, doi =
”https://doi.org/10.1016/j.dss.2010.08.006”, url =
”http://www.sciencedirect.com/science/article/pii/S0167923610001302”, author =
”E.W.T. Ngai and Yong Hu and Y.H. Wong and Yijun Chen and Xin Sun”, keywords =
”Financial fraud, Fraud detection, Literature review, Data mining, Business
intelligence”, abstract = ”This paper presents a review of — and
classification scheme for — the literature on the application of data mining
techniques for the detection of financial fraud. Although financial fraud
detection (FFD) is an emerging topic of great importance, a comprehensive
literature review of the subject has yet to be carried out. This paper thus
represents the first systematic, identifiable and comprehensive academic
literature review of the data mining techniques that have been applied to FFD.
49 journal articles on the subject published between 1997 and 2008 was
analyzed and classified into four categories of financial fraud (bank fraud,
insurance fraud, securities and commodities fraud, and other related financial
fraud) and six classes of data mining techniques (classification, regression,
clustering, prediction, outlier detection, and visualization). The findings of
this review clearly show that data mining techniques have been applied most
extensively to the detection of insurance fraud, although corporate fraud and
credit card fraud have also attracted a great deal of attention in recent
years. In contrast, we find a distinct lack of research on mortgage fraud,
money laundering, and securities and commodities fraud. The main data mining
techniques used for FFD are logistic models, neural networks, the Bayesian
belief network, and decision trees, all of which provide primary solutions to
the problems inherent in the detection and classification of fraudulent data.
This paper also addresses the gaps between FFD and the needs of the industry
to encourage additional research on neglected topics, and concludes with
several suggestions for further FFD research.”

@articleBLOMQUIST2024112695, title = Stable nodal projection method on octree
grids, journal = Journal of Computational Physics, volume = 499, pages =
112695, year = 2024, issn = 0021-9991, doi =
https://doi.org/10.1016/j.jcp.2023.112695, url =
https://www.sciencedirect.com/science/article/pii/S0021999123007908, author =
Matthew Blomquist and Scott R. West and Adam L. Binswanger and Maxime
Theillard, keywords = Incompressible Navier-Stokes, Collocated node-based,
Projection, Stability, Octree grids, Sharp interface, abstract = We propose a
novel collocated projection method for solving the incompressible Navier-
Stokes equations with arbitrary boundaries. Our approach employs non-graded
octree grids, where all variables are stored at the nodes. To discretize the
viscosity and projection steps, we utilize supra-convergent finite difference
approximations with sharp boundary treatments. We demonstrate the stability of
our projection on uniform grids, identify a sufficient stability condition on
adaptive grids, and validate these findings numerically. We further
demonstrate the accuracy and capabilities of our solver with several canonical
two- and three-dimensional simulations of incompressible fluid flows. Overall,
our method is second-order accurate, allows for dynamic grid adaptivity with
arbitrary geometries, and reduces the overhead in code development through
data collocation.

@articleNGUYEN2008763, title = Meshless methods: A review and computer
implementation aspects, journal = Mathematics and Computers in Simulation,
volume = 79, number = 3, pages = 763-813, year = 2008, issn = 0378-4754, doi =
https://doi.org/10.1016/j.matcom.2008.01.003, url =
https://www.sciencedirect.com/science/article/pii/S0378475408000062, author =
Vinh Phu Nguyen and Timon Rabczuk and Stéphane Bordas and Marc Duflot,
keywords = Meshless methods, Intrinsic enrichment, Extrinsic discontinuities,
Computer implementation, MATLAB, abstract = The aim of this manuscript is to
give a practical overview of meshless methods (for solid mechanics) based on
global weak forms through a simple and well-structured MATLAB code, to
illustrate our discourse. The source code is available for download on our
website and should help students and researchers get started with some of the
basic meshless methods; it includes intrinsic and extrinsic enrichment, point
collocation methods, several boundary condition enforcement schemes and
corresponding test cases. Several one and two-dimensional examples in
elastostatics are given including weak and strong discontinuities and testing
different ways of enforcing essential boundary conditions.

@inproceedingsSamet1988AnOO, title=An Overview of Quadtrees, Octrees, and
Related Hierarchical Data Structures, author=Hanan Samet, year=1988

@articleSARSCOV2nucleocapsid, author = Chang, C., Michalska, K., Jedrzejczak,
R., Maltseva, N., Endres, M., Godzik, A., Kim, Y., Joachimiak, A., title =
Crystal structure of RNA binding domain of nucleocapsid phosphoprotein from
SARS coronavirus 2, journal = (to be published), year = 2020, doi =
10.2210/pdb6vyo/pdb,

@miscFoldinghome, title = Folding @ home , howpublished =
http:/foldingathome.org,

@articlePolitzer1985, author = P Politzer and P R Laurence and K Jayasuriya ,
title = Molecular electrostatic potentials: an effective tool for the
elucidation of biochemical phenomena., journal = Environmental Health
Perspectives, volume = 61, number = , pages = 191-202, year = 1985, doi =
10.1289/ehp.8561191, URL =
https://ehp.niehs.nih.gov/doi/abs/10.1289/ehp.8561191, eprint =
https://ehp.niehs.nih.gov/doi/pdf/10.1289/ehp.8561191, abstract = The
electrostatic potential V(r) that is created in the space around a molecule by
its nuclei and electrons (treated as static distributions of charge) is a very
useful property for analyzing and predicting molecular reactive behavior. It
is rigorously defined and can be determined experimentally as well as
computationally. The potential has been particularly useful as an indicator of
the sites or regions of a molecule to which an approaching electrophile is
initially attracted, and it has also been applied successfully to the study of
interactions that involve a certain optimum relative orientation of the
reactants, such as between a drug and its cellular receptor. A variety of
methods for calculating V(r) is available, at different levels of rigor. For
large biologically active molecules, multipole expansions and superposition of
potentials computed for subunits have been found to be effective. A large
number of chemical and biochemical systems and processes have now been studied
in terms of electrostatic potentials. Three examples of such applications are
surveyed in this paper. These deal with: (a) reactive properties of nucleic
acids, including their component bases; (b) biological recognition processes,
including drug-receptors and enzyme-substrate interactions; and (c) chemical
carcinogenesis, referring specifically to the polycyclic aromatic hydrocarbons
and halogenated olefins and their epoxides. For each of these areas, examples
of the use of electrostatic potentials in elucidating structure-activity
patterns are given.

@ARTICLEShashikala2019,

AUTHOR=Shashikala, H. B. Mihiri and Chakravorty, Arghya and Alexov, Emil,

TITLE=Modeling Electrostatic Force in Protein-Protein Recognition,

JOURNAL=Frontiers in Molecular Biosciences,

VOLUME=6,

PAGES=94,

YEAR=2019,

URL=https://www.frontiersin.org/article/10.3389/fmolb.2019.00094,

DOI=10.3389/fmolb.2019.00094,

ISSN=2296-889X,

ABSTRACT=Electrostatic interactions are important for understanding molecular
interactions, since they are long-range interactions and can guide binding
partners to their correct binding positions. To investigate the role of
electrostatic forces in molecular recognition, we calculated electrostatic
forces between binding partners separated at various distances. The
investigation was done on a large set of 275 protein complexes using recently
developed DelPhiForce tool and in parallel, evaluating the total electrostatic
force via electrostatic association energy. To accomplish the goal, we
developed a method to find an appropriate direction to move one chain of
protein complex away from its bound position and then calculate the
corresponding electrostatic force as a function of separation distance. It is
demonstrated that at large distances between the partners, the electrostatic
force (magnitude and direction) is consistent among the protocols used and the
main factors contributing to it are the net charge of the partners and their
interfaces. However, at short distances, where partners form specific pair-
wise interactions or de-solvation penalty becomes significant, the outcome
depends on the precise balance of these factors. Based on the electrostatic
force profile (force as a function of distance), we group the cases into four
distinctive categories, among which the most intriguing is the case termed
“soft landing.” In this case, the electrostatic force at large distances is
favorable assisting the partners to come together, while at short distance it
opposes binding, and thus slows down the approach of the partners toward their
physical binding.

@articleAZN_LANCET, journal = The Lancet, year = 2020, author = Pedro M
Folegatti and Katie J Ewer and Parvinder K Aley, Brian Angus and Prof Stephan
Becker and Sandra Belij-Rammerstorfer and Duncan Bellamy and Sagida Bibi and
Mustapha Bittaye and Elizabeth A Clutterbuck and Christina Dold and Prof Saul
N Faust and Prof Adam Finn and Amy L Flaxman and Bassam Hallis and Prof Paul
Heath and Daniel Jenkin and Rajeka Lazarus and Rebecca Makinson and Angela M
Minassian and Katrina M Pollock and Maheshi Ramasamy and Hannah Robinson and
Matthew Snape and Richard Tarrant and Merryn Voysey and Catherine Green and
Alexander D Douglas and Prof Adrian V S Hill and Teresa Lambe and Sarah C
Gilbert and Andrew J Pollard, title = Safety and immunogenicity of the ChAdOx1
nCoV-19 vaccine against SARS-CoV-2: a preliminary report of a phase 1/2,
single-blind, randomised controlled trial,

@articleAZN_LANCET1, journal = The Lancet, year = 2020, author = Pedro M
Folegatti and Katie J Ewer and Parvinder K Aley, title = Safety and
immunogenicity of the ChAdOx1 nCoV-19 vaccine against SARS-CoV-2: a
preliminary report of a phase 1/2, single-blind, randomised controlled trial,

@INPROCEEDINGS6406271,

author=Y. Chen and Y. Zhou and S. Zhu and H. Xu,

booktitle=2012 International Conference on Privacy, Security, Risk and Trust
and 2012 International Confernece on Social Computing,

title=Detecting Offensive Language in Social Media to Protect Adolescent
Online Safety,

year=2012,

volume=,

number=,

pages=71-80,

@ArticleStilgoe2019, author=Stilgoe, Jack, title=Self-driving cars will take a
while to get right, journal=Nature Machine Intelligence, year=2019, month=May,
day=01, volume=1, number=5, pages=202-203, abstract=If we are to realize the
potential of self-driving cars, we need to recognize the limits of machine
learning. We should not pretend self-driving cars are around the corner: it
will still take substantial time and effort to integrate the technology safely
and fairly into our societies., issn=2522-5839, doi=10.1038/s42256-019-0046-z,
url=https://doi.org/10.1038/s42256-019-0046-z

@articleHEYDARI2022110755, title = Conservative finite volume method on
deforming geometries: The case of protein aggregation in dividing yeast cells,
journal = Journal of Computational Physics, volume = 448, pages = 110755, year
= 2022, issn = 0021-9991, doi = https://doi.org/10.1016/j.jcp.2021.110755, url
= https://www.sciencedirect.com/science/article/pii/S0021999121006501, author
= A. Ali Heydari and Suzanne S. Sindi and Maxime Theillard, keywords =
Conservative, Finite volume, Level set, Octree, Protein aggregation, Dividing
cell, abstract = Motivated by experimental observations of asymmetric protein
aggregates distributions in dividing yeast cells, we present a conservative
finite volume approach for reaction-diffusion systems defined over deforming
geometries. The key idea of our approach is to use spatio-temporal control
volumes instead of integrating the time-discretized equations in space, as it
is common practice. Both our theoretical and computational results demonstrate
the convergence of our method and highlight how traditional approaches can
lead to inaccurate solutions. We employ this novel approach to investigate the
partitioning of protein aggregates in dividing yeast cells, leveraging the
flexibility of the level set method to construct realistic biological
geometries. Using a simple reaction-diffusion model, we find that spatial
heterogeneity in yeast cells during division can alone create asymmetries in
the concentration of protein aggregates. Moreover, we find that obstructing
intracellular entities, such as nuclei or insoluble protein compartments,
amplify these asymmetries, suggesting that they may play an essential role in
regulating molecular partitioning. Beyond these findings, our results
illustrate the flexibility of our approach and its potential to design
realistic predictive tools to explore intracellular bio-mechanisms.

@ArticleD3SM01196H, author =”Fylling, Cayce and Tamayo, Joshua and Gopinath,
Arvind and Theillard, Maxime”, title =”Multi-population dissolution in
confined active fluids”, journal =”Soft Matter”, year =”2024”, volume =”20”,
issue =”7”, pages =”1392-1409”, publisher =”The Royal Society of Chemistry”,
doi =”10.1039/D3SM01196H”, url =”http://dx.doi.org/10.1039/D3SM01196H”,
abstract =”Autonomous out-of-equilibrium agents or cells in suspension are
ubiquitous in biology and engineering. Turning chemical energy into mechanical
stress, they generate activity in their environment, which may trigger
spontaneous large-scale dynamics. Often, these systems are composed of
multiple populations that may reflect the coexistence of multiple species,
differing phenotypes, or chemically varying agents in engineered settings.
Here, we present a new method for modeling such multi-population active fluids
subject to confinement. We use a continuum multi-scale mean-field approach to
represent each phase by its first three orientational moments and couple their
evolution with those of the suspending fluid. The resulting coupled system is
solved using a parallel adaptive level-set-based solver for high computational
efficiency and maximal flexibility in the confinement geometry. Motivated by
recent experimental work, we employ our method to study the spatiotemporal
dynamics of confined bacterial suspensions and swarms dominated by fluid
hydrodynamic effects. Our in silico explorations reproduce observed emergent
collective patterns, including features of active dissolution in two-
population active-passive swarms, with results clearly suggesting that
hydrodynamic effects dominate dissolution dynamics. Our work lays the
foundation for a systematic characterization and study of collective phenomena
in natural or synthetic multi-population systems such as bacteria colonies,
bird flocks, fish schools, colloid swimmers, or programmable active matter.”

@articleJIANG199013, title = Least-squares finite element method for fluid
dynamics, journal = Computer Methods in Applied Mechanics and Engineering,
volume = 81, number = 1, pages = 13-37, year = 1990, issn = 0045-7825, doi =
https://doi.org/10.1016/0045-7825(90)90139-D, url =
https://www.sciencedirect.com/science/article/pii/004578259090139D, author =
Bo-Nan Jiang and Louis A. Povinelli, abstract = This paper gives an overview
of new developments of the least-squares finite element method (LSFEM) in
fluid dynamics. Special emphasis is placed on the universality of LSFEM, the
symmetry and positiveness of the algebraic systems obtained from LSFEM, the
accommodation of LSFEM to equal-order interpolations for incompressible
viscous flows, and the natural numerical dissipation of LSFEM for convective
transport problems and high-speed compressible flows. The performance of LSFEM
is illustrated by numerical examples.

@InbookJiang1998, author=”Jiang, Bo-nan”, title=”Basis of LSFEM”,
bookTitle=”The Least-Squares Finite Element Method: Theory and Applications in
Computational Fluid Dynamics and Electromagnetics”, year=”1998”,
publisher=”Springer Berlin Heidelberg”, address=”Berlin, Heidelberg”,
pages=”47–79”, abstract=”In this chapter an attempt is made to present a
unified theory and formulation of the least-squares finite element method for
multi-dimensional problems so that in principle, it is not necessary to repeat
the same argument in different instances. In the following chapters, these
theorems, inequalities and formulations will be extensively utilized. The
mathematics is kept as simple as possible, since our purpose is to become
familiar with certain modern concepts that will help readers understand the
basic theory, formulation and properties of the least-squares method.”,
isbn=”978-3-662-03740-9”, doi=”10.1007/978-3-662-03740-9_4”,
url=”https://doi.org/10.1007/978-3-662-03740-9_4”

@articleKUCHEROVA2021110591, title = Computational modeling of protein
conformational changes - Application to the opening SARS-CoV-2 spike, journal
= Journal of Computational Physics, volume = 444, pages = 110591, year = 2021,
issn = 0021-9991, doi = https://doi.org/10.1016/j.jcp.2021.110591, url =
https://www.sciencedirect.com/science/article/pii/S0021999121004861, author =
Anna Kucherova and Selma Strango and Shahar Sukenik and Maxime Theillard,
keywords = SARS-CoV-2, Covid-19, Spike protein, Molecular trajectory, Poisson-
Boltzmann, Multiscale modeling, abstract = We present a new approach to
compute and analyze the dynamical electro-geometric properties of proteins
undergoing conformational changes. The molecular trajectory is obtained from
Markov state models, and the electrostatic potential is calculated using the
continuum Poisson-Boltzmann equation. The numerical electric potential is
constructed using a parallel sharp numerical solver implemented on adaptive
Octree grids. We introduce novel a posteriori error estimates to quantify the
solution’s accuracy on the molecular surface. To illustrate the approach, we
consider the opening of the SARS-CoV-2 spike protein using the recent
molecular trajectory simulated through the Folding@home initiative. We analyze
our results, focusing on the characteristics of the receptor-binding domain
and its vicinity. This work lays the foundation for a new class of hybrid
computational approaches, producing high-fidelity dynamical computational
measurements serving as a basis for protein bio-mechanism investigations.

@articlegooglelanguage2016, author = Yonghui Wu and Mike Schuster and Zhifeng
Chen and Quoc V. Le and Mohammad Norouzi and Wolfgang Macherey and Maxim
Krikun and Yuan Cao and Qin Gao and Klaus Macherey and Jeff Klingner and
Apurva Shah and Melvin Johnson and Xiaobing Liu and Lukasz Kaiser and Stephan
Gouws and Yoshikiyo Kato and Taku Kudo and Hideto Kazawa and Keith Stevens and
George Kurian and Nishant Patil and Wei Wang and Cliff Young and Jason Smith
and Jason Riesa and Alex Rudnick and Oriol Vinyals and Greg Corrado and
Macduff Hughes and Jeffrey Dean, title = Google’s Neural Machine Translation
System: Bridging the Gap between Human and Machine Translation, journal =
CoRR, volume = abs/1609.08144, year = 2016, url =
http://arxiv.org/abs/1609.08144, archivePrefix = arXiv, eprint = 1609.08144,
timestamp = Thu, 14 Mar 2019 09:34:18 +0100, biburl =
https://dblp.org/rec/journals/corr/WuSCLNMKCGMKSJL16.bib, bibsource = dblp
computer science bibliography, https://dblp.org

@articleWeiss2018, author = Ron J. Weiss and Jan Chorowski and Navdeep Jaitly
and Yonghui Wu and Zhifeng Chen, title = Sequence-to-Sequence Models Can
Directly Transcribe Foreign Speech, journal = CoRR, volume = abs/1703.08581,
year = 2017, url = http://arxiv.org/abs/1703.08581, archivePrefix = arXiv,
eprint = 1703.08581, timestamp = Mon, 13 Aug 2018 16:46:35 +0200, biburl =
https://dblp.org/rec/journals/corr/WeissCJWC17.bib, bibsource = dblp computer
science bibliography, https://dblp.org

@inproceedingsNobata2016, author = Nobata, Chikashi and Tetreault, Joel and
Thomas, Achint and Mehdad, Yashar and Chang, Yi, title = Abusive Language
Detection in Online User Content, year = 2016, isbn = 9781450341431, publisher
= International World Wide Web Conferences Steering Committee, address =
Republic and Canton of Geneva, CHE, url =
https://doi.org/10.1145/2872427.2883062, doi = 10.1145/2872427.2883062,
booktitle = Proceedings of the 25th International Conference on World Wide
Web, pages = 145–153, numpages = 9, keywords = hate speech, stylistic
classification, natural language processing, abusive language, nlp, discourse
classification, location = Montréal, Québec, Canada, series = WWW ’16

@articleBaskin2016, author = Igor I. Baskin and David Winkler and Igor V.
Tetko, title = A renaissance of neural networks in drug discovery, journal =
Expert Opinion on Drug Discovery, volume = 11, number = 8, pages = 785-795,
year = 2016, publisher = Taylor & Francis, doi =
10.1080/17460441.2016.1201262, note =PMID: 27295548,

URL = https://doi.org/10.1080/17460441.2016.1201262

, eprint = https://doi.org/10.1080/17460441.2016.1201262

@ArticleEkins2016, author=Ekins, Sean, title=The Next Era: Deep Learning in
Pharmaceutical Research, journal=Pharmaceutical Research, year=2016,
month=Nov, day=01, volume=33, number=11, pages=2594-2603, abstract=Over the
past decade we have witnessed the increasing sophistication of machine
learning algorithms applied in daily use from internet searches, voice
recognition, social network software to machine vision software in cameras,
phones, robots and self-driving cars. Pharmaceutical research has also seen
its fair share of machine learning developments. For example, applying such
methods to mine the growing datasets that are created in drug discovery not
only enables us to learn from the past but to predict a molecule’s properties
and behavior in future. The latest machine learning algorithm garnering
significant attention is deep learning, which is an artificial neural network
with multiple hidden layers. Publications over the last 3 years suggest that
this algorithm may have advantages over previous machine learning methods and
offer a slight but discernable edge in predictive performance. The time has
come for a balanced review of this technique but also to apply machine
learning methods such as deep learning across a wider array of endpoints
relevant to pharmaceutical research for which the datasets are growing such as
physicochemical property prediction, formulation prediction, absorption,
distribution, metabolism, excretion and toxicity (ADME/Tox), target prediction
and skin permeation, etc. We also show that there are many potential
applications of deep learning beyond cheminformatics. It will be important to
perform prospective testing (which has been carried out rarely to date) in
order to convince skeptics that there will be benefits from investing in this
technique., issn=1573-904X, doi=10.1007/s11095-016-2029-7,
url=https://doi.org/10.1007/s11095-016-2029-7

@articleGawehn2016, author = Gawehn, Erik and Hiss, Jan A. and Schneider,
Gisbert, title = Deep Learning in Drug Discovery, journal = Molecular
Informatics, volume = 35, number = 1, pages = 3-14, keywords = bioinformatics,
cheminformatics, drug design, machine-learning, neural network, virtual
screening, doi = 10.1002/minf.201501008, url =
https://onlinelibrary.wiley.com/doi/abs/10.1002/minf.201501008, eprint =
https://onlinelibrary.wiley.com/doi/pdf/10.1002/minf.201501008, abstract =
Abstract Artificial neural networks had their first heyday in molecular
informatics and drug discovery approximately two decades ago. Currently, we
are witnessing renewed interest in adapting advanced neural network
architectures for pharmaceutical research by borrowing from the field of “deep
learning”. Compared with some of the other life sciences, their application in
drug discovery is still limited. Here, we provide an overview of this emerging
field of molecular informatics, present the basic concepts of prominent deep
learning methods and offer motivation to explore these techniques for their
usefulness in computer-assisted drug discovery and design. We specifically
emphasize deep neural networks, restricted Boltzmann machine networks and
convolutional networks., year = 2016

@articleLAVECCHIA20192017, title = ”Deep learning in drug discovery:
opportunities, challenges and future prospects”, journal = ”Drug Discovery
Today”, volume = ”24”, number = ”10”, pages = ”2017 - 2032”, year = ”2019”,
issn = ”1359-6446”, doi = ”https://doi.org/10.1016/j.drudis.2019.07.006”, url
= ”http://www.sciencedirect.com/science/article/pii/S135964461930282X”, author
= ”Antonio Lavecchia”, abstract = ”Artificial Intelligence (AI) is an area of
computer science that simulates the structures and operating principles of the
human brain. Machine learning (ML) belongs to the area of AI and endeavors to
develop models from exposure to training data. Deep Learning (DL) is another
subset of AI, where models represent geometric transformations over many
different layers. This technology has shown tremendous potential in areas such
as computer vision, speech recognition and natural language processing. More
recently, DL has also been successfully applied in drug discovery. Here, I
analyze several relevant DL applications and case studies, providing a
detailed view of the current state-of-the-art in drug discovery and
highlighting not only the problematic issues, but also the successes and
opportunities for further advances.”

@articleCHEN20181241, title = ”The rise of deep learning in drug discovery”,
journal = ”Drug Discovery Today”, volume = ”23”, number = ”6”, pages = ”1241 -
1250”, year = ”2018”, issn = ”1359-6446”, doi =
”https://doi.org/10.1016/j.drudis.2018.01.039”, url =
”http://www.sciencedirect.com/science/article/pii/S1359644617303598”, author =
”Hongming Chen and Ola Engkvist and Yinhai Wang and Marcus Olivecrona and
Thomas Blaschke”, abstract = ”Over the past decade, deep learning has achieved
remarkable success in various artificial intelligence research areas. Evolved
from the previous research on artificial neural networks, this technology has
shown superior performance to other machine learning algorithms in areas such
as image and voice recognition, natural language processing, among others. The
first wave of applications of deep learning in pharmaceutical research has
emerged in recent years, and its utility has gone beyond bioactivity
predictions and has shown promise in addressing diverse problems in drug
discovery. Examples will be discussed covering bioactivity prediction, de novo
molecular design, synthesis prediction and biological image analysis.”

@ARTICLE8661629,

author=Y. Zhao and Q. Chen and W. Cao and J. Yang and J. Xiong and G. Gui,

journal=IEEE Access,

title=Deep Learning for Risk Detection and Trajectory Tracking at Construction
Sites,

year=2019,

volume=7,

number=,

pages=30905-30912,

@articlePALTRINIERI2019475, title = ”Learning about risk: Machine learning for
risk assessment”, journal = ”Safety Science”, volume = ”118”, pages = ”475 -
486”, year = ”2019”, issn = ”0925-7535”, doi =
”https://doi.org/10.1016/j.ssci.2019.06.001”, url =
”http://www.sciencedirect.com/science/article/pii/S0925753518311184”, author =
”Nicola Paltrinieri and Louise Comfort and Genserik Reniers”, keywords = ”Risk
assessment, Dynamic risk analysis, Machine learning, Deep learning”, abstract
= ”Risk assessment has a primary role in safety-critical industries. However,
it faces a series of overall challenges, partially related to technology
advancements and increasing needs. There is currently a call for continuous
risk assessment, improvement in learning past lessons and definition of
techniques to process relevant data, which are to be coupled with adequate
capability to deal with unexpected events and provide the right support to
enable risk management. Through this work, we suggest a risk assessment
approach based on machine learning. In particular, a deep neural network (DNN)
model is developed and tested for a drive-off scenario involving an Oil & Gas
drilling rig. Results show reasonable accuracy for DNN predictions and general
suitability to (partially) overcome risk assessment challenges. Nevertheless,
intrinsic model limitations should be taken into account and appropriate model
selection and customization should be carefully carried out to deliver
appropriate support for safety-related decision-making.”

@ARTICLE8316850,

author=C. Jiang and J. Song and G. Liu and L. Zheng and W. Luan,

journal=IEEE Internet of Things Journal,

title=Credit Card Fraud Detection: A Novel Approach Using Aggregation Strategy
and Feedback Mechanism,

year=2018,

volume=5,

number=5,

pages=3637-3647,

@articleDORNADULA2019631, title = ”Credit Card Fraud Detection using Machine
Learning Algorithms”, journal = ”Procedia Computer Science”, volume = ”165”,
pages = ”631 - 641”, year = ”2019”, note = ”2nd International Conference on
Recent Trends in Advanced Computing ICRTAC -DISRUP - TIV INNOVATION , 2019
November 11-12, 2019”, issn = ”1877-0509”, doi =
”https://doi.org/10.1016/j.procs.2020.01.057”, url =
”http://www.sciencedirect.com/science/article/pii/S187705092030065X”, author =
”Vaishnavi Nath Dornadula and S Geetha”, keywords = ”Card-Not-Present frauds,
Card-Present-Frauds, Concept Drift”, abstract = ”Credit card frauds are easy
and friendly targets. E-commerce and many other online sites have increased
the online payment modes, increasing the risk for online frauds. Increase in
fraud rates, researchers started using different machine learning methods to
detect and analyse frauds in online transactions. The main aim of the paper is
to design and develop a novel fraud detection method for Streaming Transaction
Data, with an objective, to analyse the past transaction details of the
customers and extract the behavioural patterns. Where cardholders are
clustered into different groups based on their transaction amount. Then using
sliding window strategy [1], to aggregate the transaction made by the
cardholders from different groups so that the behavioural pattern of the
groups can be extracted respectively. Later different classifiers
[3],[5],[6],[8] are trained over the groups separately. And then the
classifier with better rating score can be chosen to be one of the best
methods to predict frauds. Thus, followed by a feedback mechanism to solve the
problem of concept drift [1]. In this paper, we worked with European credit
card fraud dataset.”

@articleRAISSI2018125, title = ”Hidden physics models: Machine learning of
nonlinear partial differential equations”, journal = ”Journal of Computational
Physics”, volume = ”357”, pages = ”125 - 141”, year = ”2018”, issn =
”0021-9991”, doi = ”https://doi.org/10.1016/j.jcp.2017.11.039”, url =
”http://www.sciencedirect.com/science/article/pii/S0021999117309014”, author =
”Maziar Raissi and George Em Karniadakis”, keywords = ”Probabilistic machine
learning, System identification, Bayesian modeling, Uncertainty
quantification, Fractional equations, Small data”, abstract = ”While there is
currently a lot of enthusiasm about “big data”, useful data is usually “small”
and expensive to acquire. In this paper, we present a new paradigm of learning
partial differential equations from small data. In particular, we introduce
hidden physics models, which are essentially data-efficient learning machines
capable of leveraging the underlying laws of physics, expressed by time
dependent and nonlinear partial differential equations, to extract patterns
from high-dimensional data generated from experiments. The proposed
methodology may be applied to the problem of learning, system identification,
or data-driven discovery of partial differential equations. Our framework
relies on Gaussian processes, a powerful tool for probabilistic inference over
functions, that enables us to strike a balance between model complexity and
data fitting. The effectiveness of the proposed approach is demonstrated
through a variety of canonical problems, spanning a number of scientific
domains, including the Navier–Stokes, Schrödinger, Kuramoto–Sivashinsky, and
time dependent linear fractional equations. The methodology provides a
promising new direction for harnessing the long-standing developments of
classical methods in applied mathematics and mathematical physics to design
learning machines with the ability to operate in complex domains without
requiring large quantities of data.” @articleRAISSI2019686, title = ”Physics-
informed neural networks: A deep learning framework for solving forward and
inverse problems involving nonlinear partial differential equations”, journal
= ”Journal of Computational Physics”, volume = ”378”, pages = ”686 - 707”,
year = ”2019”, issn = ”0021-9991”, doi =
”https://doi.org/10.1016/j.jcp.2018.10.045”, url =
”http://www.sciencedirect.com/science/article/pii/S0021999118307125”, author =
”M. Raissi and P. Perdikaris and G.E. Karniadakis”, keywords = ”Data-driven
scientific computing, Machine learning, Predictive modeling, Runge–Kutta
methods, Nonlinear dynamics”, abstract = ”We introduce physics-informed neural
networks – neural networks that are trained to solve supervised learning tasks
while respecting any given laws of physics described by general nonlinear
partial differential equations. In this work, we present our developments in
the context of solving two main classes of problems: data-driven solution and
data-driven discovery of partial differential equations. Depending on the
nature and arrangement of the available data, we devise two distinct types of
algorithms, namely continuous time and discrete time models. The first type of
models forms a new family of data-efficient spatio-temporal function
approximators, while the latter type allows the use of arbitrarily accurate
implicit Runge–Kutta time stepping schemes with unlimited number of stages.
The effectiveness of the proposed framework is demonstrated through a
collection of classical problems in fluids, quantum mechanics,
reaction–diffusion systems, and the propagation of nonlinear shallow-water
waves.”

@articleBotu2015, author = Botu, Venkatesh and Ramprasad, Rampi, title =
Adaptive machine learning framework to accelerate ab initio molecular
dynamics, journal = International Journal of Quantum Chemistry, volume = 115,
number = 16, pages = 1074-1083, keywords = ab initio molecular dynamics,
accelerate, machine learning, fingerprint, adaptive, doi = 10.1002/qua.24836,
url = https://onlinelibrary.wiley.com/doi/abs/10.1002/qua.24836, eprint =
https://onlinelibrary.wiley.com/doi/pdf/10.1002/qua.24836, abstract = Quantum
mechanics-based ab initio molecular dynamics (MD) simulation schemes offer an
accurate and direct means to monitor the time evolution of materials.
Nevertheless, the expensive and repetitive energy and force computations
required in such simulations lead to significant bottlenecks. Here, we lay the
foundations for an accelerated ab initio MD approach integrated with a machine
learning framework. The proposed algorithm learns from previously visited
configurations in a continuous and adaptive manner on-the-fly, and predicts
(with chemical accuracy) the energies and atomic forces of a new configuration
at a minuscule fraction of the time taken by conventional ab initio methods.
Key elements of this new accelerated ab initio MD paradigm include
representations of atomic configurations by numerical fingerprints, a learning
algorithm to map the fingerprints to the properties, a decision engine that
guides the choice of the prediction scheme, and requisite amount of ab initio
data. The performance of each aspect of the proposed scheme is critically
evaluated for Al in several different chemical environments. This work has
enormous implications beyond ab initio MD acceleration. It can also lead to
accelerated structure and property prediction schemes, and accurate force
fields. © 2014 Wiley Periodicals, Inc., year = 2015

@ArticleSmith2019, author=Smith, Justin S. and Nebgen, Benjamin T. and
Zubatyuk, Roman and Lubbers, Nicholas and Devereux, Christian and Barros,
Kipton and Tretiak, Sergei and Isayev, Olexandr and Roitberg, Adrian E.,
title=Approaching coupled cluster accuracy with a general-purpose neural
network potential through transfer learning, journal=Nature Communications,
year=2019, month=Jul, day=01, volume=10, number=1, pages=2903,
abstract=Computational modeling of chemical and biological systems at atomic
resolution is a crucial tool in the chemist’s toolset. The use of computer
simulations requires a balance between cost and accuracy: quantum-mechanical
methods provide high accuracy but are computationally expensive and scale
poorly to large systems, while classical force fields are cheap and scalable,
but lack transferability to new systems. Machine learning can be used to
achieve the best of both approaches. Here we train a general-purpose neural
network potential (ANI-1ccx) that approaches CCSD(T)/CBS accuracy on
benchmarks for reaction thermochemistry, isomerization, and drug-like
molecular torsions. This is achieved by training a network to DFT data then
using transfer learning techniques to retrain on a dataset of gold standard QM
calculations (CCSD(T)/CBS) that optimally spans chemical space. The resulting
potential is broadly applicable to materials science, biology, and chemistry,
and billions of times faster than CCSD(T)/CBS calculations., issn=2041-1723,
doi=10.1038/s41467-019-10827-4, url=https://doi.org/10.1038/s41467-019-10827-4

@articleMoreau2006, author = Moreau,A. and Teytaud,O. and Bertoglio,J. P. ,
title = Optimal estimation for large-eddy simulation of turbulence and
application to the analysis of subgrid models, journal = Physics of Fluids,
volume = 18, number = 10, pages = 105101, year = 2006, doi =
10.1063/1.2357974,

URL = https://doi.org/10.1063/1.2357974

, eprint = https://doi.org/10.1063/1.2357974

@articlefukami_fukagata_taira_2019, title=Super-resolution reconstruction of
turbulent flows with machine learning, volume=870, DOI=10.1017/jfm.2019.238,
journal=Journal of Fluid Mechanics, publisher=Cambridge University Press,
author=Fukami, Kai and Fukagata, Koji and Taira, Kunihiko, year=2019,
pages=106–120

@articleUNGER20081994, title = ”Coupling of scales in a multiscale simulation
using neural networks”, journal = ”Computers & Structures”, volume = ”86”,
number = ”21”, pages = ”1994 - 2003”, year = ”2008”, issn = ”0045-7949”, doi =
”https://doi.org/10.1016/j.compstruc.2008.05.004”, url =
”http://www.sciencedirect.com/science/article/pii/S0045794908001430”, author =
”Jörg F. Unger and Carsten Könke”, keywords = ”Neural network, Homogenization,
Multilayer perceptron, Support vector machines, Constitutive modelling”,
abstract = ”Multiscale approaches require the coupling of models on different
spatial scales. In this paper, a coupling using neural networks is proposed.
Based on a set of mesoscale simulations of concrete a system of neural
networks is trained to approximate the response. A macroscale constitutive
model is obtained by homogenizing the mesoscale response. Special focus is put
on the mesh sensitivity, since the mesoscale model includes softening and
consequently the width of the localization zone compared to the dimension of
the mesoscale model is an additional parameter in the model.”

@articleMingTryggvason2015, author = Ma,Ming and Lu,Jiacai and
Tryggvason,Gretar , title = Using statistical learning to close two-fluid
multiphase flow equations for a simple bubbly system, journal = Physics of
Fluids, volume = 27, number = 9, pages = 092101, year = 2015, doi =
10.1063/1.4930004,

URL = https://aip.scitation.org/doi/abs/10.1063/1.4930004

, eprint = https://aip.scitation.org/doi/pdf/10.1063/1.4930004

@articleGIBOU2019442, title = ”Sharp interface approaches and deep learning
techniques for multiphase flows”, journal = ”Journal of Computational
Physics”, volume = ”380”, pages = ”442 - 463”, year = ”2019”, issn =
”0021-9991”, doi = ”https://doi.org/10.1016/j.jcp.2018.05.031”, url =
”http://www.sciencedirect.com/science/article/pii/S0021999118303371”, author =
”Frederic Gibou and David Hyde and Ron Fedkiw”, keywords = ”Multiphase flows,
Ghost-fluid method, Voronoi interface method, Surface tension, Quadtree and
octree, Parallel”, abstract = ”We present a review on numerical methods for
simulating multiphase and free surface flows. We focus in particular on
numerical methods that seek to preserve the discontinuous nature of the
solutions across the interface between phases. We provide a discussion on the
Ghost-Fluid and Voronoi Interface methods, on the treatment of surface tension
forces that avoid stringent time step restrictions, on adaptive grid
refinement techniques for improved efficiency and on parallel computing
approaches. We present the results of some simulations obtained with these
treatments in two and three spatial dimensions. We also provide a discussion
of Machine Learning and Deep Learning techniques in the context of multiphase
flows and propose several future potential research thrusts for using deep
learning to enhance the study and simulation of multiphase flows.”

@articleChang6591, author = Chang, Fi-John and Yang, Han-Chung and Lu, Jau-Yau
and Hong, Jian-Hao, title = Neural network modelling for mean velocity and
turbulence intensities of steep channel flows, journal = Hydrological
Processes, volume = 22, number = 2, pages = 265-274, keywords = artificial
neural network (ANN), velocity profile, turbulent open channel flow, log-law
model, Reynolds stress model, doi = 10.1002/hyp.6591, url =
https://onlinelibrary.wiley.com/doi/abs/10.1002/hyp.6591, eprint =
https://onlinelibrary.wiley.com/doi/pdf/10.1002/hyp.6591, abstract = Abstract
The main purpose of this study is to evaluate the potential of simulating the
profiles of the mean velocity and turbulence intensities for the steep open
channel flows over a smooth boundary using artificial neural networks. In a
laboratory flume, turbulent flow conditions were measured using a fibre-optic
laser doppler velocimeter (FLDV). One thousand and sixty-four data sets were
collected for different slopes and aspect ratios at different locations. These
data sets were randomly split into two subsets, i.e. training and validation
sets. The multi-layer functional link network (MFLN) was used to construct the
simulation model based on the training data. The constructed MFLN models can
almost perfectly simulate the velocity profile and turbulence intensity. The
values of correlation coefficient (γ) are close to one and the values of root
mean square error (RMSE) are close to zero in all conditions. The results
demonstrate that the MFLN can precisely simulate the velocity profiles, while
the log law and Reynolds stress model (RSM) are less effective when used to
simulate the velocity profiles close to the side wall. The simulated
longitudinal turbulence intensities yielded by the MFLN were also fairly
consistent with the measured data, while the simulated vertical turbulence
intensities by the RSM were not consistent with the measured data. Copyright ©
2007 John Wiley & Sons, Ltd., year = 2008

@articleMILANO20021, title = ”Neural Network Modeling for Near Wall Turbulent
Flow”, journal = ”Journal of Computational Physics”, volume = ”182”, number =
”1”, pages = ”1 - 26”, year = ”2002”, issn = ”0021-9991”, doi =
”https://doi.org/10.1006/jcph.2002.7146”, url =
”http://www.sciencedirect.com/science/article/pii/S0021999102971469”, author =
”Michele Milano and Petros Koumoutsakos”, abstract = ”A neural network
methodology is developed in order to reconstruct the near wall field in a
turbulent flow by exploiting flow fields provided by direct numerical
simulations. The results obtained from the neural network methodology are
compared with the results obtained from prediction and reconstruction using
proper orthogonal decomposition (POD). Using the property that the POD is
equivalent to a specific linear neural network, a nonlinear neural network
extension is presented. It is shown that for a relatively small additional
computational cost nonlinear neural networks provide us with improved
reconstruction and prediction capabilities for the near wall velocity fields.
Based on these results advantages and drawbacks of both approaches are
discussed with an outlook toward the development of near wall models for
turbulence modeling and control.”

@articleGiralt2000, author = Giralt,Francesc and Arenas,A. and Ferre-Giné,J.
and Rallo,R. and Kopp,G. A. , title = The simulation and interpretation of
free turbulence with a cognitive neural system, journal = Physics of Fluids,
volume = 12, number = 7, pages = 1826-1835, year = 2000, doi =
10.1063/1.870430,

URL = https://doi.org/10.1063/1.870430

, eprint = https://doi.org/10.1063/1.870430

@ArticleAsproulis2013, author=Asproulis, Nikolaos and Drikakis, Dimitris,
title=An artificial neural network-based multiscale method for hybrid
atomistic-continuum simulations, journal=Microfluidics and Nanofluidics,
year=2013, month=Oct, day=01, volume=15, number=4, pages=559-574,
abstract=This paper presents an artificial neural network-based multiscale
method for coupling continuum and molecular simulations. Molecular dynamics
modelling is employed as a local “high resolution” refinement of computational
data required by the continuum computational fluid dynamics solver. The
coupling between atomistic and continuum simulations is obtained by an
artificial neural network (ANN) methodology. The ANN aims to optimise the
transfer of information through minimisation of (1) the computational cost by
avoiding repetitive atomistic simulations of nearly identical states, and (2)
the fluctuation strength of the atomistic outputs that are fed back to the
continuum solver. Results are presented for prototype flows such as the
isothermal Couette flow with slip boundary conditions and the slip Couette
flow with heat transfer., issn=1613-4990, doi=10.1007/s10404-013-1154-4,
url=https://doi.org/10.1007/s10404-013-1154-4

@articleZhenwei2015, title = Molecular Dynamics with On-the-Fly Machine
Learning of Quantum-Mechanical Forces, author = Li, Zhenwei and Kermode, James
R. and De Vita, Alessandro, journal = Phys. Rev. Lett., volume = 114, issue =
9, pages = 096405, numpages = 5, year = 2015, month = Mar, publisher =
American Physical Society, doi = 10.1103/PhysRevLett.114.096405, url =
https://link.aps.org/doi/10.1103/PhysRevLett.114.096405

@article Chmielae1603015, author = Chmiela, Stefan and Tkatchenko, Alexandre
and Sauceda, Huziel E. and Poltavsky, Igor and Schütt, Kristof T. and Müller,
Klaus-Robert, title = Machine learning of accurate energy-conserving molecular
force fields, volume = 3, number = 5, elocation-id = e1603015, year = 2017,
doi = 10.1126/sciadv.1603015, publisher = American Association for the
Advancement of Science, abstract = Using conservation of energy—a fundamental
property of closed classical and quantum mechanical systems—we develop an
efficient gradient-domain machine learning (GDML) approach to construct
accurate molecular force fields using a restricted number of samples from ab
initio molecular dynamics (AIMD) trajectories. The GDML implementation is able
to reproduce global potential energy surfaces of intermediate-sized molecules
with an accuracy of 0.3 kcal mol-1 for energies and 1 kcal mol-1 Å˝̊-1 for
atomic forces using only 1000 conformational geometries for training. We
demonstrate this accuracy for AIMD trajectories of molecules, including
benzene, toluene, naphthalene, ethanol, uracil, and aspirin. The challenge of
constructing conservative force fields is accomplished in our work by learning
in a Hilbert space of vector-valued functions that obey the law of energy
conservation. The GDML approach enables quantitative molecular dynamics
simulations for molecules at a fraction of cost of explicit AIMD calculations,
thereby allowing the construction of efficient force fields with the accuracy
and transferability of high-level ab initio methods., URL =
https://advances.sciencemag.org/content/3/5/e1603015, eprint =
https://advances.sciencemag.org/content/3/5/e1603015.full.pdf, journal =
Science Advances

@Articlecomputation8010015, AUTHOR = Frank, Michael and Drikakis, Dimitris and
Charissis, Vassilis, TITLE = Machine-Learning Methods for Computational
Science and Engineering, JOURNAL = Computation, VOLUME = 8, YEAR = 2020,
NUMBER = 1, ARTICLE-NUMBER = 15, URL = https://www.mdpi.com/2079-3197/8/1/15,
ISSN = 2079-3197, ABSTRACT = The re-kindled fascination in machine learning
(ML), observed over the last few decades, has also percolated into natural
sciences and engineering. ML algorithms are now used in scientific computing,
as well as in data-mining and processing. In this paper, we provide a review
of the state-of-the-art in ML for computational science and engineering. We
discuss ways of using ML to speed up or improve the quality of simulation
techniques such as computational fluid dynamics, molecular dynamics, and
structural analysis. We explore the ability of ML to produce computationally
efficient surrogate models of physical applications that circumvent the need
for the more expensive simulation techniques entirely. We also discuss how ML
can be used to process large amounts of data, using as examples many different
scientific fields, such as engineering, medicine, astronomy and computing.
Finally, we review how ML has been used to create more realistic and
responsive virtual reality applications., DOI = 10.3390/computation8010015

@articleDESPRES2020109275, title = ”Machine Learning design of Volume of Fluid
schemes for compressible flows”, journal = ”Journal of Computational Physics”,
volume = ”408”, pages = ”109275”, year = ”2020”, issn = ”0021-9991”, doi =
”https://doi.org/10.1016/j.jcp.2020.109275”, url =
”http://www.sciencedirect.com/science/article/pii/S0021999120300498”, author =
”Bruno Després and Hervé Jourdren”, keywords = ”VOF, CFD, ML”, abstract = ”Our
aim is to establish the feasibility of Machine-Learning-designed Volume of
Fluid algorithms for compressible flows. We detail the incremental steps of
the construction of a new family of Volume of Fluid-Machine Learning (VOF-ML)
schemes adapted to bi-material compressible Euler calculations on Cartesian
grids. An additivity principle is formulated for the Machine Learning
datasets. We explain a key feature of this approach which is how to adapt the
compressible solver to the preservation of natural symmetries. The VOF-ML
schemes show good accuracy for advection of a variety of interfaces, including
regular interfaces (straight lines and arcs of circle), Lipschitz interfaces
(corners) and non Lipschitz triple point (the Trifolium test problem). Basic
comparisons with a SLIC/Downwind scheme are presented together with elementary
bi-material calculations with shocks.”

@article Brunton3932, author = Brunton, Steven L. and Proctor, Joshua L. and
Kutz, J. Nathan, title = Discovering governing equations from data by sparse
identification of nonlinear dynamical systems, volume = 113, number = 15,
pages = 3932–3937, year = 2016, doi = 10.1073/pnas.1517384113, publisher =
National Academy of Sciences, abstract = Understanding dynamic constraints and
balances in nature has facilitated rapid development of knowledge and enabled
technology, including aircraft, combustion engines, satellites, and electrical
power. This work develops a novel framework to discover governing equations
underlying a dynamical system simply from data measurements, leveraging
advances in sparsity techniques and machine learning. The resulting models are
parsimonious, balancing model complexity with descriptive ability while
avoiding overfitting. There are many critical data-driven problems, such as
understanding cognition from neural recordings, inferring climate patterns,
determining stability of financial markets, predicting and suppressing the
spread of disease, and controlling turbulence for greener transportation and
energy. With abundant data and elusive laws, data-driven discovery of dynamics
will continue to play an important role in these efforts.Extracting governing
equations from data is a central challenge in many diverse areas of science
and engineering. Data are abundant whereas models often remain elusive, as in
climate science, neuroscience, ecology, finance, and epidemiology, to name
only a few examples. In this work, we combine sparsity-promoting techniques
and machine learning with nonlinear dynamical systems to discover governing
equations from noisy measurement data. The only assumption about the structure
of the model is that there are only a few important terms that govern the
dynamics, so that the equations are sparse in the space of possible functions;
this assumption holds for many physical systems in an appropriate basis. In
particular, we use sparse regression to determine the fewest terms in the
dynamic governing equations required to accurately represent the data. This
results in parsimonious models that balance accuracy with model complexity to
avoid overfitting. We demonstrate the algorithm on a wide range of problems,
from simple canonical systems, including linear and nonlinear oscillators and
the chaotic Lorenz system, to the fluid vortex shedding behind an obstacle.
The fluid example illustrates the ability of this method to discover the
underlying dynamics of a system that took experts in the community nearly 30
years to resolve. We also show that this method generalizes to parameterized
systems and systems that are time-varying or have external forcing., issn =
0027-8424, URL = https://www.pnas.org/content/113/15/3932, eprint =
https://www.pnas.org/content/113/15/3932.full.pdf, journal = Proceedings of
the National Academy of Sciences

@ArticleYan2020, author=Yan, Li and Zhang, Hai-Tao and Goncalves, Jorge and
Xiao, Yang and Wang, Maolin and Guo, Yuqi and Sun, Chuan and Tang, Xiuchuan
and Jing, Liang and Zhang, Mingyang and Huang, Xiang and Xiao, Ying and Cao,
Haosen and Chen, Yanyan and Ren, Tongxin and Wang, Fang and Xiao, Yaru and
Huang, Sufang and Tan, Xi and Huang, Niannian and Jiao, Bo and Cheng, Cheng
and Zhang, Yong and Luo, Ailin and Mombaerts, Laurent and Jin, Junyang and
Cao, Zhiguo and Li, Shusheng and Xu, Hui and Yuan, Ye, title=An interpretable
mortality prediction model for COVID-19 patients, journal=Nature Machine
Intelligence, year=2020, month=May, day=01, volume=2, number=5, pages=283-288,
abstract=The sudden increase in COVID-19 cases is putting high pressure on
healthcare services worldwide. At this stage, fast, accurate and early
clinical assessment of the disease severity is vital. To support decision
making and logistical planning in healthcare systems, this study leverages a
database of blood samples from 485 infected patients in the region of Wuhan,
China, to identify crucial predictive biomarkers of disease mortality. For
this purpose, machine learning tools selected three biomarkers that predict
the mortality of individual patients more than 10 days in advance with more
than 90% accuracy: lactic dehydrogenase (LDH), lymphocyte and high-sensitivity
C-reactive protein (hs-CRP). In particular, relatively high levels of LDH
alone seem to play a crucial role in distinguishing the vast majority of cases
that require immediate medical attention. This finding is consistent with
current medical knowledge that high LDH levels are associated with tissue
breakdown occurring in various diseases, including pulmonary disorders such as
pneumonia. Overall, this Article suggests a simple and operable decision rule
to quickly predict patients at the highest risk, allowing them to be
prioritized and potentially reducing the mortality rate., issn=2522-5839,
doi=10.1038/s42256-020-0180-7, url=https://doi.org/10.1038/s42256-020-0180-7

@ArticleVinyals2019, author=Vinyals, Oriol and Babuschkin, Igor and Czarnecki,
Wojciech M. and Mathieu, Michaël and Dudzik, Andrew and Chung, Junyoung and
Choi, David H. and Powell, Richard and Ewalds, Timo and Georgiev, Petko and
Oh, Junhyuk and Horgan, Dan and Kroiss, Manuel and Danihelka, Ivo and Huang,
Aja and Sifre, Laurent and Cai, Trevor and Agapiou, John P. and Jaderberg, Max
and Vezhnevets, Alexander S. and Leblond, Rémi and Pohlen, Tobias and
Dalibard, Valentin and Budden, David and Sulsky, Yury and Molloy, James and
Paine, Tom L. and Gulcehre, Caglar and Wang, Ziyu and Pfaff, Tobias and Wu,
Yuhuai and Ring, Roman and Yogatama, Dani and Wünsch, Dario and McKinney,
Katrina and Smith, Oliver and Schaul, Tom and Lillicrap, Timothy and
Kavukcuoglu, Koray and Hassabis, Demis and Apps, Chris and Silver, David,
title=Grandmaster level in StarCraft II using multi-agent reinforcement
learning, journal=Nature, year=2019, month=Nov, day=01, volume=575,
number=7782, pages=350-354, abstract=Many real-world applications require
artificial agents to compete and coordinate with other agents in complex
environments. As a stepping stone to this goal, the domain of StarCraft has
emerged as an important challenge for artificial intelligence research, owing
to its iconic and enduring status among the most difficult professional
esports and its relevance to the real world in terms of its raw complexity and
multi-agent challenges. Over the course of a decade and numerous
competitions1–3, the strongest agents have simplified important aspects of the
game, utilized superhuman capabilities, or employed hand-crafted sub-systems4.