Higher order nonlocal operator method

TL;DR

This paper introduces a higher order nonlocal operator method that enhances the original approach by achieving higher convergence rates and simultaneously obtaining multiple derivatives without shape functions, demonstrated through various numerical examples.

Contribution

The paper extends the nonlocal operator method to higher order schemes using Taylor series, improving convergence and derivative computation without shape functions.

Findings

Achieves higher order convergence compared to the original method.

Simultaneously obtains multiple derivatives up to a specified order.

Demonstrates effectiveness through numerical examples in strong and weak forms.

Abstract

We extend the nonlocal operator method to higher order scheme by using a higher order Taylor series expansion of the unknown field. Such a higher order scheme improves the original nonlocal operator method proposed by the authors in [A nonlocal operator method for solving partial differential equations], which can only achieve one-order convergence. The higher order nonlocal operator method obtains all partial derivatives with specified maximal order simultaneously without resorting to shape functions. The functional based on the nonlocal operators converts the construction of residual and stiffness matrix into a series of matrix multiplication on the nonlocal operator matrix. Several numerical examples solved by strong form or weak form are presented to show the capabilities of this method.

Tables

Table 1. Table 1: Statistical results for 2 dimensional Poisson equation

Nnode	$\Delta x$	$L_2$ norm	$\frac{u_{max}}{u_{max}^{exact}}-1$	$p$-order	$p^{hg}$
1681	0.025	0.0485	-0.0281	1	1
1681	0.025	0.0262	0.01	2	1
1681	0.025	0.0139	-0.00256	3	1
1681	0.025	0.0175	-0.00308	4	1
6561	0.0125	0.0379	0.033	1	0
6561	0.0125	0.0179	0.0714	1	1
6561	0.0125	0.011	0.00505	2	1
25921	0.00625	0.0202	0.0221	1	1
25921	0.00625	0.00501	0.00266	2	1
25921	0.00625	0.00191	-0.000417	3	1
40401	0.005	0.00777	-0.00263	1	1
160801	0.0025	0.00291	0.0007	1	1

Table 2. Table 2: Statistical results for 3 dimensional Poisson equation

Nnode	$\Delta x$	$L_2$ norm	$\frac{u_{max}}{u_{max}^{exact}}-1$	$p$-order	$p^{hg}$
10648	0.04763	0.0907	-0.0406	1	1
29791	0.03333	0.0604	-0.0248	1	1
68921	0.025	0.0485	-0.02	1	1

Table 3. Table 3: Statistical results for 4 dimensional Poisson equation

Nnode	$\Delta x$	$L_2$ norm	$\frac{u_{max}}{u_{max}^{exact}}-1$	$p$-order	$p^{hg}$
14641	0.1	0.169	-0.0514	1	1
65536	0.0667	0.118	-0.0171	1	1
160000	0.0526	0.0983	-0.0203	1	1
810000	0.0345	0.0579	0.00304	1	1
2560000	0.0256	0.0454	0.00152	1	1

Table 4. Table 4: Statistical results for 5 dimensional Poisson equation

Nnode	$\Delta x$	$L_2$ norm	$\frac{u_{max}}{u_{max}^{exact}}-1$	$p$-order	$p^{hg}$
7776	0.2	0.229	-0.114	1	1
100000	0.111	0.181	-0.0944	1	1
1048576	0.0667	0.13	-0.0485	1	1
4084101	0.05	0.0985	-0.0352	1	1

Table 5. Table 5: Nearly incompressible block: displacement w m ​ a ​ x subscript 𝑤 𝑚 𝑎 𝑥 w_{max} (mm)

H1 element	13.17 ($8^3$ mesh)	19.52 (${32}^3$ mesh)
H2 element	19.54 ($8^3$ mesh)	20.01 (${32}^3$ mesh)
NOM	19.14 (${11}^3$ nodes)	20.43 (${21}^3$ nodes)

Table 6. Table 6: Minimal number of neighbors in support. d 𝑑 d =number of spatial dimensions, n 𝑛 n =maximal order of derivatives

$N_d^n$	n=1	n=2	n=3	n=4	n=5	n=6
d=1	1	2	3	4	5	6
d=2	2	5	9	14	20	27
d=3	3	9	19	34	55	83
d=4	4	14	34	69	125	209
d=5	5	20	55	125	251	461
d=6	6	27	83	209	461	923

Equations

$${\mathcal{S}}_i^{\prime }=\{{\mathbf{x}}_j|{\mathbf{x}}_i\in {\mathcal{S}}_j\}.$$

$$\displaystyle\tilde{\nabla}\mathbf{v}_{i}:=\int_{\mathcal{S}_{i}}w(\mathbf{r})\mathbf{v}_{ij}\otimes\mathbf{r}\text{d}V_{j}\cdot\Big{(}\int_{\mathcal{S}_{i}}w(\mathbf{r})\mathbf{r}\otimes\mathbf{r}\text{d}V_{j}\Big{)}^{-1}.$$

$$\tilde{\nabla }{\mathbf{v}}_i$$

$$\tilde{\nabla }\delta {\mathbf{v}}_i$$

$${\mathcal{F}}_i^{hg}$$

$$u_j=u_i+\sum _{(n_1,...,n_d)\in {\alpha }_d^n}\frac{r_1^{n_1}...r_d^{n_d}}{n_1!...n_d!}{u}_{i,n_1...n_d}+O(r^{|\alpha |+1})$$

$$\mathbf{r}$$

$$u_{i,n_1...n_d}$$

$$|\alpha |$$

$${\alpha }_d^n$$

$$\text{mbox}or{\alpha }_d^n$$

$$u_{i,n_1...n_d},\frac{r_1^{n_1}...r_d^{n_d}}{n_1!...n_d!},\forall (n_1,...,n_d)\in {\alpha }_d^n.$$

$$r_1^{n_1}...r_d^{n_d}\propto h^{n_1+...+n_d}$$

$$u_j$$

$$=u_i+\sum _{(n_1,...,n_d)\in {\alpha }_d^n}\frac{r_1^{n_1}...r_d^{n_d}}{h_i^{n_1+...+n_d}}{u}_{i,n_1...n_d}^h+O(r^{n+1})$$

$$u_{i,n_1...n_d}^h=\frac{h_i^{n_1+...+n_d}}{n_1!...n_d!}{u}_{i,n_1...n_d}$$

$${\mathbf{p}}_j^h$$

$${\partial }_{\alpha }^h{u}_i$$

$${\partial }_{\alpha }{u}_i$$

$${\partial }_{\alpha }{u}_i={\mathbf{H}}_i^{-1}{\partial }_{\alpha }^h{u}_i$$

$$\displaystyle\mathbf{H}_{i}=\mbox{diag}\big{[}h_{i},...,\frac{h_{i}^{n_{1}+...+n_{d}}}{n_{1}!...n_{d}!},...,\frac{h_{i}^{n}}{n!}\big{]}$$

$$u_{ij}=({\partial }_{\alpha }^h{u}_i{)}^T{\mathbf{p}}_j^h,\forall j\in {\mathcal{S}}_i$$

$${\int }_{{\mathcal{S}}_i}w(\mathbf{r})u_{ij}({\mathbf{p}}_j^h{)}^T\text{d}{V}_j$$

$$=({\partial }_{\alpha }{u}_i{)}^T{\mathbf{H}}_i{\int }_{{\mathcal{S}}_i}w(\mathbf{r}){\mathbf{p}}_j^h\otimes ({\mathbf{p}}_j^h{)}^T\text{d}{V}_j$$

$$\displaystyle\tilde{\partial}_{\alpha}u_{i}:=\mathbf{H}_{i}^{-1}\Big{(}\int_{\mathcal{S}_{i}}w(\mathbf{r})\mathbf{p}^{h}_{j}\otimes(\mathbf{p}^{h}_{j})^{T}\text{d}V_{j}\Big{)}^{-1}\int_{\mathcal{S}_{i}}w(\mathbf{r})u_{ij}\mathbf{p}^{h}_{j}\text{d}V_{j}=\mathbf{K}_{i}\cdot\int_{\mathcal{S}_{i}}w(\mathbf{r})\mathbf{p}^{h}_{j}u_{ij}\text{d}V_{j}$$

$$\displaystyle\mathbf{K}_{i}:=\mathbf{H}_{i}^{-1}\Big{(}\int_{\mathcal{S}_{i}}w(\mathbf{r})\mathbf{p}^{h}_{j}\otimes(\mathbf{p}^{h}_{j})^{T}\text{d}V_{j}\Big{)}^{-1}.$$

$${\tilde{\partial }}_{\alpha }\delta u_i:={\mathbf{K}}_i\cdot {\int }_{{\mathcal{S}}_i}w(\mathbf{r}){\mathbf{p}}_j^h(\delta u_j-\delta u_i)\text{d}{V}_j$$

$$\Omega =\sum _{i=1}^N\Delta V_i$$

$${\mathcal{S}}_i$$

$${\tilde{\partial }}_{\alpha }{u}_i$$

$${\tilde{\partial }}_{\alpha }\delta u_i$$

$${\mathbf{K}}_i$$

$${\mathbf{p}}_{wi}^h$$

$$\Delta {\mathbf{u}}_i$$

$${\tilde{\partial }}_{\alpha }{u}_i={\mathbf{K}}_i{\mathbf{p}}_{wi}^h\Delta {\mathbf{u}}_i={\mathbf{B}}_{\alpha i}{\mathbf{u}}_i$$

$${\mathbf{B}}_{\alpha i}$$

Full text

1 Introduction

In the field of solving Partial Differential Equations (PDEs), methods can be generally divided into (semi-)analytical methods and numerical methods. Analytical methods include the method of separation of variables [1], integral transforms [2], Homotopy Analysis Method (HAM) [3], Variational Iteration Method (VIM) [4] and so on. Analytical methods have advantages in finding the approximate/exact solutions but are often restricted to regular geometry domain. The numerical methods contain Rayleigh-Ritz method, Finite Difference Method(FDM), Finite Element Methods(FEMs), Meshless Methods(MMs), isogeometric analysis [5], to just name a few. In finite element methods, the computation domain is meshed into discrete elements and the shape function defined on the element is used to interpolate the field value within the element. Meshless methods comprise many different formulations [6, 7], for example, Smoothed Particle Hydrodynamics (SPH) [8, 9], Element-Free Galerkin method (EFG)[10], Reproducing Kernel Particle Method (RKPM) [11] and so on. Finite element methods and most of the meshless methods interpolate the field value in the domain by means of shape functions, and the derivatives in PDEs are constructed from the derivatives of the shape functions. Different from the methods by interpolation technique, finite difference method expresses the partial derivatives with finite difference. However, the finite difference method is only applicable for domain with regular geometry. For a more complete review of the PDEs by numerical methods, we refer to [12].

When it comes to higher order PDEs in higher dimensional space, finite element method, meshless methods and finite difference method confront some problems. For finite element methods, the topology of element in higher dimensional space is complicated. Though the simplex element is valid in any dimensions, the representation of the topology and calculation of the shape functions and their partial derivatives are cumbersome. Other difficulties involve the numerical integration and the continuity required on the interface between adjoint elements. Nevertheless, some finite element schemes are developed for arbitrary order of derivative (i.e. [13, 14, 15]). For meshless methods based on the shape functions, there is no problem for the mesh construction in the higher dimensional space. However, the numerical integration in meshless methods requires a background mesh, which is the same as the finite element methods. What’s worse, the calculation of higher order derivative of the shape function is very expensive. One method to circumvent the numerical integration in background mesh is the nodal integration, which however surfers the rank-deficiency problem. The finite difference method can construct higher order finite difference to replace the higher order partial derivatives, but the stencil becomes more complicated. Other problems with higher order PDEs in high dimensional space involve the complicated boundary conditions at different orders of derivatives, the proof on uniqueness, robust, stability of the solution.

The fundamental elements in PDEs are various partial differential operators of different orders. How to deal with these operators is the central topic of various numerical methods. FEMs and most meshless methods start from the shape function for interpolation, while the derivatives of shape function are used to represent the differential operators. Such process is expensive for higher order differential derivatives in higher dimensions. The difficulties to numerically describe the differential operators arise from the locality of the operator, where the locality denotes the operator being defined at a point. To circumvent the difficulties arising from locality, Nonlocal Operator Method (NOM) was proposed by the authors [16]. NOM starts from the common differential operators such as gradient, curl, divergence and Hessian operators, to define the nonlocal gradient, nonlocal curl, nonlocal divergence and nonlocal Hessian operators by introducing the support with finite characteristic length. These nonlocal operators can be viewed as the generalization of the local operators. When the support degenerates to one point, the nonlocal operators recover the local operators. Unlike FEMs, meshless methods or finite difference method, NOM is a “true” meshless method and only requires the neighbor list in the support in order to construct the nonlocal derivatives. The low order nonlocal operators [16] can solve low order (not more than 4th-order) PDEs, but not higher order PDEs.

The purpose of the paper is to develop a higher order nonlocal operator method for solving higher order PDEs of multiple fields in multiple spatial dimensions. The nonlocal operator method obtains a set of partial derivatives of different orders at once. Combining with weighed residual method and variational principles, nonlocal operator method establishes the residual and tangent stiffness matrix for PDEs by some matrix operation on common terms, operator matrix. In contrast with finite element method or meshless method with shape functions, the nonlocal operator method leads to the differential operators directly and adopts the nodal integration method. The remainder of the paper is outlined as follows. In section 2, the basic concepts such as support and dual-support, and the low order nonlocal operators are reviewed and then the higher order nonlocal operator method based on Taylor series expansion of multiple variables is developed. We define a special quadratic functional to derive the nonlocal strong form for a $2n$-order PDEs based on the nonlocal operators in section 3. We give several numerical examples to demonstrate the capabilities of this method in solving PDEs by strong form in section 4 and by weak form in section 5. Finally, we conclude in section 6.

2 Nonlocal operator method

2.1 Basic concepts

Consider a domain as shown in Fig.1, let $\mathbf{x}_{i}$ be spatial coordinates in the domain $\mathbf{\Omega}$; $\mathbf{r}:=\mathbf{x}_{j}-\mathbf{x}_{i}$ is a spatial vector starts from $\mathbf{x}_{i}$ to $\mathbf{x}_{j}$; $\mathbf{v}_{i}:=\mathbf{v}(\mathbf{x}_{i},t)$ and $\mathbf{v}_{j}:=\mathbf{v}(\mathbf{x}_{j},t)$ are the field value for $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$, respectively; $\mathbf{v}_{ij}:=\mathbf{v}_{j}-\mathbf{v}_{i}$ is the relative field vector for spatial vector $\mathbf{r}$.

Support $\mathcal{S}_{i}$ of point $\mathbf{x}_{i}$ is the domain where any spatial point $\mathbf{x}_{j}$ forms spatial vector $\mathbf{r}(=\mathbf{x}_{j}-\mathbf{x}_{i})$ from $\mathbf{x}_{i}$ to $\mathbf{x}_{j}$. The support serves as the basis for the nonlocal operators. There is no restriction on the support shapes, which can be spherical domain, cube, semi-spherical domain, triangle and so on.

Dual-support is defined as a union of the points whose supports include $\mathbf{x}$, denoted by

[TABLE]

Point $\mathbf{x}_{j}$ forms dual-vector $\mathbf{r}^{\prime}(=\mathbf{x}_{i}-\mathbf{x}_{j}=-\mathbf{r})$ in $\mathcal{S}_{i}^{\prime}$. On the other hand, $\mathbf{r}^{\prime}$ is the spatial vector formed in $\mathcal{S}_{j}$. One example to illustrate the support and dual-support is shown in Fig.1.

The nonlocal operator method uses the basic nonlocal operators to replace the local operator in calculus such as the gradient, divergence, curl and Hessian operators. The functional formulated by the local differential operator can be used to construct the residual or tangent stiffness matrix by replacing the local operator with the corresponding nonlocal operator. However, convergence rate of the original nonlocal operator is limited to 1 since the basic nonlocal operator is one-order.

The nonlocal gradient of a vector field $\mathbf{v}$ for point $\mathbf{x}_{i}$ in support $\mathcal{S}_{i}$ is defined as

[TABLE]

The nonlocal gradient operator and its variation in discrete form are

[TABLE]

The operator energy functional for vector field at point $\mathbf{x}_{i}$ is

[TABLE]

where $p^{hg}$ is the penalty coefficient. The residual and tangent stiffness matrix of $\mathcal{F}_{i}^{hg}$ can be obtained with ease, we refer to [16] for more details.

2.2 Higher order nonlocal operator method

Several formulations of the Taylor series expansion of a function of multiple variables are available in A. A scalar field $u_{j}$ at a point $j\in\mathcal{S}_{i}$ can be obtained by the Taylor series expansion at $u_{i}$ in $d$ dimensions with maximal derivative order not more than $n$,

[TABLE]

where

[TABLE]

$\alpha_{d}^{n}$ is the list of flattened multi-indexes, where $d$ denotes the number of spatial dimensions and $n$ is the maximal order of partial derivative for one index. Two special multi-index can be written as

[TABLE]

where $\mathbb{N}^{0}=\{0,1,2,3,...\}$. Eq.11 gives a multi-index with $(1+n)^{d}-1$ elements and $|\alpha|=nd$, while the multi-index by Eq.10 has $\frac{(n+d)!}{n!d!}-1$ elements according to Combinatorics. In this paper, we adopt the multi-index by Eq.10 since it avoids the mixed higher order terms and has some benefit for numerical computation. The way to obtain all elements in $\alpha_{d}^{n}$ of Eq.10 by the Mathematica sees B.

For any multi-index $(n_{1},...,n_{d})\in\alpha_{d}^{n}$, the partial derivative and the polynomial are

[TABLE]

However, the original form of Taylor series expansion is very sensitive to the round-off error. For example,

[TABLE]

where $h$ is the characteristic length scale of the support. The higher order terms reduce to 0 quickly when $h<1$, or explode as $h>1$. It is expected to have the length scale $h$ approaching $1$. When length scale of support $\mathcal{S}_{i}$ at $u_{i}$ is taken into account, Taylor series expansion by Eq.6 can be written as

[TABLE]

where $h_{i}$ is the characteristic length of $\mathcal{S}_{i}$, and

[TABLE]

Let $\mathbf{p}^{h}_{j}$, $\partial^{h}_{\alpha}u_{i}$ and $\partial_{\alpha}u_{i}$ be the list of the flattened polynomials, scaled partial derivatives, partial derivatives, respectively, based on multi-index notation $\alpha_{d}^{n}$ in Eq.10,

[TABLE]

Introducing $h$ in Eq.15 enables the terms in Eq.15 being in the “same” characteristic length scale. The actual partial derivatives can be recovered by

[TABLE]

where

[TABLE]

where diag[$a_{1},...,a_{n}$] denotes a diagonal matrix whose diagonal entries starting in the upper left corner are $a_{1},...,a_{n}$.

Therefore, Taylor series expansion with $u_{i}$ being moved to left side of the equation can be written as

[TABLE]

where $u_{ij}=u_{j}-u_{i}$.

Integrate $u_{ij}$ with weighted coefficient $w(\mathbf{r})(\mathbf{p}^{h}_{j})^{T}$ in support $\mathcal{S}_{i}$, we obtain

[TABLE]

where $w(\mathbf{r})$ is the weight function.

Therefore, the nonlocal operator $\tilde{\partial}_{\alpha}u_{i}$ can be obtained as

[TABLE]

where

[TABLE]

The reason to call Eq.22 nonlocal operator is that it is defined in the support, in contrast with the local operator defined at a point. The nonlocal operator approximates the local operator with order up to $|\alpha|$. Traditional local operator is suitable for theoretical derivation but not for numerical analysis since its definition is limited to infinitesimal. The nonlocal operator can be viewed as a generalization of the conventional local operator.

The variation of $\tilde{\partial}_{\alpha}u_{i}$ is

[TABLE]

In the continuous form, the number of dimensions of $\partial\delta u_{i}$ is infinite and discretization is required. After discretization of the domain by particles, the whole domain is represented by

[TABLE]

where $i$ is the global index of volume $\Delta V_{i}$, $N$ is the number of particles in $\Omega$.

Particles in $\mathcal{S}_{i}$ are represented by

[TABLE]

where $j_{1},...,j_{k},...,j_{n_{i}}$ are the global indexes of neighbors of particle $i$, $n_{i}$ is the number of neighbors of $i$ in $\mathcal{S}_{i}$.

The discrete form of Eq.22 and its variation are

[TABLE]

where

[TABLE]

When the weight function $w(\mathbf{r})$ is selected as the reciprocal of the volume, Eq.29 and Eq.30 can be simplified further. The nonlocal operator provides all the partial derivatives with maximal order for single index up to $n$. The set of derivatives in PDEs of real application is a subset of the nonlocal operator. It should be noted that when the number of points in support is the same as the length of multi-index $\alpha_{d}^{n}$ and the coefficient matrix from Eq.20 for all points in support is well conditioned, the nonlocal operator can be obtained directly by the inverse of the coefficient matrix. In this case, the nonlocal operator serves as an efficient way to obtain the higher order finite difference scheme.

Each term in $\tilde{\partial}_{\alpha}u_{i}$ corresponds to the row of $\mathbf{K}_{i}\mathbf{p}^{h}_{wi}$ multiplying $\Delta\mathbf{u}_{i}$. Eq.27 can be used to replace the differential operators in PDEs to form the algebraic equations. This way is through strong form of the PDEs. The other ways to solve the linear (nonlinear) PDEs are through the weak formulations (weighted residual method) or the variational formulations (i.e. [16]). In these cases, the variation of $\partial_{\alpha}u_{i}$ in Eq.28 is required.

Eq.27 can be written more concisely as

[TABLE]

with $\mathbf{B}_{\alpha i}$ being the operator matrix for point $i$ based on multi-index $\alpha_{d}^{n}$

[TABLE]

where $(1,\cdots,1)_{n_{p}}\mathbf{K}_{i}\mathbf{p}_{wi}^{h}$ is the column sum of $\mathbf{K}_{i}\mathbf{p}_{wi}^{h}$, $n_{p}$ is the length of $\alpha_{d}^{n}$. The operator matrix obtains all the partial derivatives of maximal order less than $|\alpha|+1$ by the nodal values in support. For real applications, one can select the specific rows in the operator matrix based on the partial derivatives contained in the specific PDEs. The template acts as

The traditional differential operator and their combination of one order or higher order and the corresponding variations can be constructed from Eq.27 and Eq.28, respectively. For example, the multi-index, polynomials and partial derivatives in two dimensions with maximal second-order derivatives are

[TABLE]

For the case of Poisson equation in 2D, $\nabla^{2}u=f$. In the strong form, the operator $\nabla^{2}u=\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}}$ is required, one can select the $\partial_{\alpha}u_{i}[2]$ in Eq.35 for $\frac{\partial^{2}u}{\partial y^{2}}$ and $\partial_{\alpha}u_{i}[5]$ in Eq.35 for $\frac{\partial^{2}u}{\partial x^{2}}$. When solved in weak form, one can select the $\partial_{\alpha}u_{i}[1]$ in Eq.35 for $\frac{\partial u}{\partial y}$ and $\partial_{\alpha}u_{i}[3]$ in Eq.35 for $\frac{\partial u}{\partial x}$ to construct the tangent stiffness matrix.

In fact, the nonlocal operator $\partial_{\alpha}u_{i}$ in discrete form can be obtained by least squares. Consider the weighted square sum of the Taylor series expansion in $\mathcal{S}_{i}$,

[TABLE]

$\frac{\partial\mathcal{F}_{i}(\mathbf{u})}{\partial(\tilde{\partial}_{\alpha}u_{i})}=0$ leads to

[TABLE]

which is the same as Eq.27.

Meanwhile, Eq.36 represents the operator energy functional in nonlocal operator method, and can be used to construct the tangent stiffness matrix of operator energy functional. The operator energy functional is the quadratic functional of the Taylor series expansion. Through Eq.38, Eq.37 can be simplified into

[TABLE]

where

[TABLE]

The first and second variation of $\mathcal{F}_{i}(\mathbf{u})$ read

[TABLE]

Let $\mathbf{v}_{i}(j)=\sum_{k=1}^{n_{i}}\mathbf{M}_{i}(j,k)$ be the sum of row of matrix $\mathbf{M}_{i}$, then the tangent stiffness matrix can be extracted from Eq.43,

[TABLE]

where the first row(column) denotes the entries for point $i$, while the neighbors start from the second row(column), $p_{hg}$ is the penalty coefficient and $m_{i}$ the normalization coefficient

[TABLE]

where $\mathbf{r}$ varies for each $j$.

Let $n_{i}$ be the number of neighbors in $\mathcal{S}_{i}$ and $n_{p}$ be the length of $\tilde{\partial}^{h}_{\alpha}u_{i}$. The dimensions of terms in $\mathbf{M}_{i}$ are

[TABLE]

When $n_{i}<n_{p}$, $\sum_{j\in\mathcal{S}_{i}}w(\mathbf{r})\mathbf{p}^{h}_{j}(\mathbf{p}_{j}^{h})^{T}\Delta V_{j}$ is singular. It is required that $n_{i}\geq n_{p}$, so that ${\mathbf{M}_{i}}$ in Eq.41 is well defined. The number of neighbors is selected as $5p+n_{p}$, where $p$ denotes the order of the nonlocal operator. These extra nodes are used to overcome the rank deficiency in the nodal integration.

The operator energy functional $\mathcal{F}(\mathbf{u})$ represents the topology of the nonlocal operator method. Any field derived from $\mathcal{F}(\mathbf{u})$ should try to satisfy $\mathcal{F}(\mathbf{u})=0$ at the first step, which is independent with the actual physical model to be solved.

3 Quadratic functional

A very special functional has the form

[TABLE]

where $\mathbf{D}$ is an arbitrary symmetric matrix, $\tilde{\partial}u\subset\tilde{\partial}_{\alpha}u$ in Eq.27. The operator matrix $\mathbf{B}$ is constructed from $\mathbf{B}_{\alpha}$ based on the index of terms $\partial u$ in $\partial_{\alpha}u$. Some examples of Eq.46 are given in §3.1.

When $\mathbf{D}$ is independent with the unknown functions $u$, the functional $\mathcal{F}(\tilde{\partial}u)$ is pure quadratic, the first and second variation of $\mathcal{F}(\tilde{\partial}u)$ at a point are

[TABLE]

and the residual and tangent stiffness matrix at a point can be written as

[TABLE]

When $\mathbf{D}:=\mathbf{D}(u,\partial u)$ is nonlinear tensor, the functional can be converted into quadratic functional by linearization and the Newton-Raphson can be employed to find the solution.

According to $\tilde{\partial}u_{i}\subset\tilde{\partial}_{\alpha}u_{i}$ in Eq.22, we write $\tilde{\partial}u_{i}=\mathbf{K}_{i}^{\prime}\cdot\int_{\mathcal{S}_{i}}w(\mathbf{r})\mathbf{p}_{j}^{h}u_{ij}\text{d}V_{j}$, where $\mathbf{K}_{i}^{\prime}\subset\mathbf{K}_{i}$ in Eq.22. Let $\bm{\sigma}_{i}=\mathbf{D}\tilde{\partial}u_{i}$ and consider the variation in domain

[TABLE]

Note that $\mathbf{p}_{i}^{h}$ in $\mathcal{S}^{\prime}_{i}$ varies for different $j$ since $\mathbf{p}_{i}^{h}$ is computed in $j$’s support $\mathcal{S}_{j}$.

The terms with $\delta u_{i}$ in the first order variation $\delta\mathcal{F}=0$ are

[TABLE]

with “equivalent” higher order partial differential term $-\partial^{T}(\mathbf{D}\partial u)$, where $\partial:=(...,\frac{\partial^{n_{1}+...+n_{d}}}{\partial x_{1}^{n_{1}}...\partial x_{d}^{n_{d}}},...)^{T}$ is the differential operator based on subset of multi-index $\alpha_{d}^{n}$ in Eq.10. PDE given by $-\partial^{T}(\mathbf{D}\partial u)$ has a maximal differential order of $2n$. The nonlocal strong form by Eq.51 can be solved directly by explicit integration algorithm. It should be noted that $\bm{\sigma},\tilde{\partial}u$ are in form of column vector for a scalar field $u$. The generalization of $u$ to vector field is straightforward.

Eq.51 alone may suffer numerical instabilities (zero-energy mode), and therefore the operator energy functional by Eq.36 is required. Eq.51 with correction terms can be written as

[TABLE]

3.1 Elastic solid materials

In this section, we give some examples on how to express the linear/nonlinear elastic mechanics by the form of nonlocal operator method. The maximal derivative order in linear elastic mechanics is 2 and the corresponding weak form only requires first order partial derivative. The internal energy functional for plane stress, plane strain and 3D linear elastic solid at a point are

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The tangent stiffness matrix of that point can be extracted by performing the first or second order variation of the above functionals.

For nonlinear elastic material, the strain energy density is a function of the deformation gradient, i.e.

[TABLE]

while $\mathbf{F}$ is consisted with the nonlocal operators in $\tilde{\partial}\mathbf{u}_{3d}$

[TABLE]

Within the framework of total Lagrangian formulation, the first Piola-Kirchhoff stress is the direct derivative of the strain energy over the deformation gradient,

[TABLE]

Furthermore, the material tensor (stress-strain relation) which is required in the implicit analysis can be obtained with the derivative of the first Piola-Kirchhoff stress,

[TABLE]

The 4th order material tensor $\mathbf{D}_{4}$ can be expressed in matrix form when the deformation gradient is flattened.

[TABLE]

where the flattened deformation gradient and first Piola-Kirchhoff stress are

[TABLE]

and

[TABLE]

For the case of nearly incompressible Neo-Hooke material [17], the strain energy can be expressed as

[TABLE]

where $J=\det{\mathbf{F}}$.

The first Piola-Kirchoff stress is

[TABLE]

With some derivation, the material tensor in matrix form can be written as

[TABLE]

where $J_{,F}$ is the vector form of $J_{,\mathbf{F}}$, and

[TABLE]

The numerical example based on material model Eq.85 is given in section 5.4.

4 Numerical examples by strong form

The nonlocal operator defined in Eq.27 can be used to replace the partial derivatives of different orders in the partial differential equation. In other word, we can use the nonlocal operator to solve the PDE by its strong form. In this sense, the nonlocal operator is similar to the finite difference method. However, finite difference scheme of different order is constructed on the regular grid, where the extension to higher dimensions or higher order derivative require special treatment, while the nonlocal operator is established simply based on the neighbor list in the support. In this section, we test the accuracy of nonlocal operator in solving second order ordinary differential equation (ODE) or PDE by strong form. Note that the operator energy functional is not required in solving PDE by strong form.

The first three numerical examples demonstrate the capabilities of nonlocal operator method in obtaining high order finite difference scheme.

4.1 Second-order ODE

The ODE with boundary condition is given by

[TABLE]

with analytic solution

[TABLE]

Since the highest order derivative in the ODE is two, the order of derivative in the nonlocal operator list should be $p\geq 2$. We test the nonlocal operator with $p=2,3,4,5,6$ in solving the second-order ODE. The minimal number of neighbors in the support is selected as the number of terms in the nonlocal operator. The difference between numerical result and theoretical solution is measured by the L2-norm, which is calculated by

[TABLE]

The convergence of the L2-norm for $u$ is shown in Fig.2.

It can be seen that with the increase of order in the nonlocal operator, the convergence rate increases greatly. $p=2,3$ have the same convergence rate.

4.2 1D Schrödinger equation

This section tests the accuracy of the eigenvalue problem in 1D. The Schrödinger equation written in adimensional units for a one-dimensional harmonic oscillator is

[TABLE]

For simplicity, we use $\omega=1$. The particles are uniformly distributed with constant spacing $\Delta x$ on the region [-10,10].

The exact wave functions and eigenvalues can be expressed as

[TABLE]

where $n$ is a non-negative integer. $H_{n}(x)$ is the $n$-order Hermite polynomial. We calculate the lowest eigenvalue and compare the numerical result with $\lambda_{0}=0.5$. The convergence plot of the error is shown in Fig.3.

4.3 Poisson equation

In this section, we test the Poisson equation

[TABLE]

with the boundary conditions

[TABLE]

The analytic solution is

[TABLE]

The number of neighbors for each point is selected as the number of terms in the nonlocal operator. We test the convergence of the L2 error for the $u$ field under uniform discretizations and non-uniform discretization in Fig.4. The convergent plot is given in Fig.5.

5 Numerical examples by weak form

The fourth example aims at solving the Poisson equation in higher dimensional space by both the “equivalent” integral form and operator energy functional. The fifth example is about the biharmornic equation. The sixth example solves the Von-Karmon plate with simply support.

5.1 Poisson equation in higher dimensional space

In this section, we solve the Poisson equation in $n$ dimensional space by nonlocal operator method. The $n$ dimensional Poisson equation is

[TABLE]

with analytic solution

[TABLE]

under the boundary conditions

[TABLE]

where $\mathbf{x}=(x_{1},...,x_{n})$, f(\mathbf{x})=\nabla^{2}\big{(}\exp{(\sum_{i=1}^{n}(-1)^{i-1}x_{i})}\Pi_{i=1}^{n}x_{i}(1-x_{i})\big{)}.

The equivalent integral functional for Eq.104 is

[TABLE]

The tangent stiffness is constructed from the operator matrix $\mathbf{B}$ for $\nabla u$, e.g.

[TABLE]

The Dirichlet boundary condition are applied by penalty method.

The Poisson equation with dimensional number $n=(2,3,4,5)$ under different discretization and order of nonlocal operator are tested, where the statistical results are shown in Tables. (1,2,3,4).

Table.1 gives the statistical results for 2D Poisson equation under different discretizations. When the order of nonlocal operator increased from 1 to 3, the L2 norm and error for $u_{max}$ decrease gradually as shown in several cases. However, for 4-order nonlocal operator, the result is not better than 3-order scheme. The 3-order scheme with 25921 nodes can achieve better result than 1-order scheme with 160801 nodes. The comparison between 5,6 rows shows that the operator energy functional has positive effect in improving the accuracy. In contrast with the scheme by strong form, the convergence property of weak form is slightly affected by the operator energy functional.

For the 3D Poisson equation, we tested three cases with discretization ranged from 22, 31, 41 nodes in each direction. The statistical results are given in Table.3. The L2 norm and error for $u_{max}$ decrease with the point grid space. When 41 nodes used for each direction, the L2 norm is approximately 5%.

For the 4 dimensional Poisson equation, we tested four cases with discretization ranged from 11,16,20,30,40 nodes for each direction. The statistical results are given in Table.3. The L2 norm and error for $u_{max}$ decrease with the point grid space. When 40 nodes used for each direction, the L2 norm is approximately 5%.

For 5 dimensional Poisson equation, when 16 nodes are assigned in each direction, the number of nodes reaches 1,048,576. More nodes in each direction will lead to the dimension disaster. The statistical results for different discretization are given in Table 4. The L2 norm and error for maximal $u$ decrease with the node spacing. We tested maximal 21 nodes in each direction (the computational scale is restricted by the computational power of a desktop PC), the L2 norm is approximately 9.85% and the error for $u_{max}$ with respect to the theoretical solution is less than 4%.

5.2 Square plate with simple support

The plate equation reads

[TABLE]

where $D_{0}=\frac{Et^{3}}{12(1-\nu^{2})}$, with Dirichlet boundary conditions

[TABLE]

The analytic solution for the simply support square plate subjected to uniform load is denoted by [18]

[TABLE]

where $\alpha_{m}=\frac{m\pi}{2}$.

The “equivalent” integral form for Eq.109 is

[TABLE]

where

[TABLE]

[TABLE]

The parameters for the plate include length $a=1$, thickness $t=0.01$ m and uniform pressure $q_{0}$=-100 N, Poisson ratio $\nu=0.3$, elastic modulus $E=30$ GPa and $D_{0}=\frac{Et^{3}}{12(1-\nu^{2})}$

With the aid of nonlocal operator $\tilde{\partial}w$ and its operator matrix $\mathbf{B}$, the first and second variation of the energy functional are,

[TABLE]

where $\delta\mathbf{w}_{i}$ is the vector for all unknowns in support $\mathcal{S}_{i}$.

The plate is discretized uniformly and the number of neighbors for each point is selected as $n=5p+length(\tilde{\partial}u)$, where $p$ is the order of nonlocal operator. The deflection curves for several discretizations are compared with the analytic solution in Fig.7. The contour of the deflection field for discretization of $40\times 40$ is shown in Fig.7. Compared with the original nonlocal operator method, the higher order NOM obtains the nonlocal operator in a simper way.

5.3 Von Kármá equations for a thin plate

The Von Kármá equations [19] are a set of nonlinear partial differential equations describing the large deflections of thin flat plates. The equations are based on Kirchhoff hypothesis : the surface normals to the plane of the plate remain perpendicular to the plate after deformation and in-plane (membrane) displacements are small and the change in thickness of the plate is negligible. These assumptions imply that the displacement field $v$ in the plate can be expressed as [20],

[TABLE]

For a plate of a thickness $h$ defined on the mid-surface $(x_{1},x_{2})$, Von Kármá energy is given by [21, 22]

[TABLE]

where Laplace operator $\Delta w=\frac{\partial^{2}w}{\partial x_{1}^{2}}+\frac{\partial^{2}w}{\partial x_{2}^{2}}$ and

[TABLE]

$\varepsilon_{ij}$ is the strain tensor with nonlinear terms in the deformations $u_{1}=u_{1}(x_{1},x_{2}),u_{2}=u_{2}(x_{1},x_{2}),w=w(x_{1},x_{2})$:

[TABLE]

where $(u_{1},u_{2})$ is the lateral displacement field due to membrane effect, $w$ is the deflection, $\sigma_{ij}$ is the stress tensor, linearly proportional to $\varepsilon_{ij}$, $E=30\times 10^{6}$ Pa and $\nu=0.3$ are the Young modulus and the Poisson ratio, respectively, $q=1000$Pa is the external normal force per unit area of the plate. The dimensions of the plate are $1.0\times 1.0\times 0.01\mbox{ m}^{3}$. The energy functional in Eq.118 leads to the governing equations

[TABLE]

The Cauchy stress tensor in mid-plane can be written as

[TABLE]

In this paper, we write the moment and curvature by tensor form. The conventional vectorial form can be recovered with ease. The moment tenor and curvature tensor are

[TABLE]

[TABLE]

where $D_{0}=\frac{Eh^{3}}{12(1-\nu^{2})}$. The moment tensor is similar to the stress tensor in the plane stress conditions.

The rotation in direction $\bm{n}$ is

[TABLE]

The curvature in direction $\bm{n}$ is

[TABLE]

The momentum in direction $\bm{n}$ is

[TABLE]

The nonlocal differential operators in Eq.118 can be written as

[TABLE]

The gradient of energy functional on $\tilde{\partial}u$ is

[TABLE]

The Hessian matrix of $\mathcal{F}$ can be obtained with ease by computing $\frac{\partial^{2}\mathcal{F}}{\partial{\tilde{\partial}}\mathbf{u}^{2}}$. The solution can be obtained when using the Newton-Raphson method in C. For simplicity, we only consider the simple support boundary conditions.

The plate solved by NOM is discretized by $50\times 50$ nodes. The reference results are calculated by S4R plate/shell element in ABAQUS [23]. S4R element is a 4-node doubly curved thin or thick shell element with reduced integration, hourglass control, finite membrane strains. In ABAQUS, the flat thin plate with the same material parameters are discretized into $100\times 100$ elements.

Displacement in membrane and deflection out-of-plane for nodes on $y=0.5$ under different load levels are depicted in Figs.8,9, respectively, where the lines represent the results by ABAQUS while the discrete symbols are the results by NOM. The displacement results agree well with that by ABAQUS.

Maximal central deflection is plotted in Fig.10, which shows the non-linearity increase with load level significantly. It can be seen that the result by NOM matches well with by finite element method.

5.4 Nearly incompressible block

In this section, we model the nearly incompressible block of material consitituion in Eq.85 by nonlocal operator method with Newton-Raphson iteration method. The nearly incompressible block of height $h=50$ mm, length $2h$ and width $2h$ is loaded by an equally distributed pressure $p=3$ MPa at its top center of area $h\times h$ mm2, as shown in Fig.11. For symmetry reason, only a quarter of the block is modeled. The bottom face is fixed in $z$-direction, while the nodes on plane $y=0$ are fixed in $y$-direction and the nodes on plane $x=0$ are fixed in $x$-direction, as similarly presented in reference [17]. The material parameters are $\kappa=499.92568$ MPa, $\mu=1.61148$ MPa.

The deformed block at final load level is depicted in Fig.12, where good agreement is obtained between the finite element method and the NOM. The maximal displacements in $z$-direction by linear hexahedral element (H1), quadratic hexahedral elements (H2) and nonlocal operator method are given in Table.5.

6 Concluding Remarks

We have proposed a higher order nonlocal operator method for solving higher order PDEs based on the strong form or the equivalent integral formulated by weighed residual method or variational principles.

The relation of nonlocal operator and local operator is that the local operator is defined on a point, while the nonlocal operator method is defined on the support with finite characteristic length scale. When the support decreases to a point, the nonlocal operator degenerates to the local operator. Nonlocal operator is constructed from the Taylor series expansion and approximates the local derivative with orders up to $n$. In order to establish the nonlocal operator, only finite points in support is required. Nonlocal operator can be viewed as a generalization of the local operator. Most rules applied to the local operator can be adopted directly by the nonlocal operator method.

In certain cases such as the regular grid, the nonlocal operator method is similar to the finite difference. One difference with finite difference method is that finite difference method requires a regular grid. When handling multiple fields, the finite difference method should adopt staggered grid for the reason of numerical stability, which complicates the numerical implementation. For nonlocal operator method, all the nodes have the same functions, in contrast with the finite difference method with staggered grid, where different nodes represent different fields. In terms of numerical stability, the nonlocal operator method monitors and enhances the robust of the derivative estimation by the operator energy functional, the quadratic functional of the Taylor series expansion. When adding the quadratic functional of the Taylor series expansion to the functional of physical problem, the numerical stability can be enhanced and traced.

Taylor series expansion of multiple variables based on the mutli-index notation is powerful in deriving various partial derivatives of different orders. Multi-index notation $\alpha_{d}^{n}$ in B can obtain automatically all the partial derivatives with order up to $n$ in $d$ spatial dimensions. In addition, the characteristic length scale is introduced for high precision of the derivative estimation. With all the partial derivatives available, all linear PDEs up to $2n$ orders can be described with ease. By replacing the differential operator with the nonlocal one, nonlocal operator method converts the PDEs into algebraic equations directly. The nonlocal operator method can be viewed as a tool to study the higher order PDEs.

Acknowledgments

The first author acknowledge the supports from the COMBAT Program (Computational Modeling and Design of Lithium-ion Batteries, Grant No.615132). The supports from National Basic Research Program of China (973 Program: 2011CB013800) and NSFC (51474157), the Ministry of Science and Technology of China (Grant No.SLDRCE14-B-28, SLDRCE14-B-31) are acknowledged.

Appendix A Taylor series expansion

There are several formulations for the Taylor series expansion of function of multiple variables. The conventional Taylor series of a function at origin can be written as [25]

[TABLE]

By using the generalization of inner product, the Taylor series expansion is

[TABLE]

where $\mathbf{r}=\mathbf{x}_{j}-\mathbf{x}_{i}$, $\mathbf{r}^{n}=\mathbf{r}\otimes...\otimes\mathbf{r}$, and $\cdot^{(n)}$ is the generalization of inner product, where two special cases are $\cdot^{(1)}=\cdot,\,\cdot^{(2)}=:$.

Or using the $d$-dimensional multi-index, the Taylor series expansion is

[TABLE]

where $\alpha=\{(\alpha_{1},...,\alpha_{d})|\alpha_{i}\in\mathbb{N}^{0},1\leq i\leq d\}$. In this paper, Eq.134 is adopt for Taylor series expansion.

Appendix B Mathematica code for multi-index

For the multi-indexes in

[TABLE]

the Mathematica code with high efficiency is

MultiIndexList[d_,n_]:=Module[{a,b,c},a=Subsets[Range[d+n],{d}]; Do[c=a[[i]];b=c-1;b[[2;;]]-=c[[1;;-2]];a[[i]]=b,{i,Length[a]}];a[[2;;]]]; (note: d=number of spatial dimensions, n=maximal order of derivative)

The number of elements in $\alpha_{d}^{n}$ can be determined by counting the combination of positive integer $k$ as a sum of $d$ non-negative integers up to non-commutativity. Imagine a line of $d+k-1$ positions, where each position can contain either a cat or a divider. If one has $k$ (nameless) cats and $d-1$ dividers, he can split the cats into $d$ groups by choosing positions for the dividers: $C_{k+d-1}^{d-1}=C_{k+d-1}^{k}$, where $C_{n}^{d}$ is binomial coefficient and can be written as $C_{n}^{d}=\binom{n}{d}=\frac{n!}{d!(n-d)!}$. The size of each group of cats corresponds to one of the non-negative integers in the sum.

Therefore, in $d$-dimension space, the number of $k$-order derivatives by Eq.133 is

[TABLE]

The number of all derivatives with maximal order $n$ in $d$ dimensional space is

[TABLE]

In order to obtain the $n$-order derivatives, the number of neighbors in $\mathcal{S}_{i}$ must be not less than $C_{n+d}^{n}-1$, so that the coefficient matrix for all derivatives is invertible. If more neighbors are in the support, the least square method can be used to find the approximation. The minimal number of neighbors in support is listed in Table.6.

Appendix C Newton-Raphson method for nonlinear functional

The core of NOM is the functional, which comprises with physical functional and the operator energy functional. The physical functional may contains the functional on the domain and other functional on the boundaries. In all,

[TABLE]

The first and second derivative on all unknowns lead to the residual and the tangent stiffness matrix, respectively

[TABLE]

When any term in $\mathcal{F}$ is nonlinear functional, the Newton-Raphson is required. The solution is updated by iteration in each step. In the $n$ step, the residual $\mathbf{R}(\mathbf{u}_{n})=0$ is satisfied, $\mathbf{R}(\mathbf{u}_{n+1})$ in the next step can be approximated by Taylor series expansion

[TABLE]

The solution in $n+1$ step can be obtained by the iterations

[TABLE]

where $k$ denotes the iteration number in $n+1$ step, $\mathbf{u}^{0}=\mathbf{u}_{n}$, $\mathbf{u}^{k+1}=\mathbf{u}^{k}+\Delta\mathbf{u}^{k+1}$. When

[TABLE]

the iteration converges.