Truncation error estimates of approximate operators in a generalized particle method

TL;DR

This paper derives and verifies truncation error estimates for approximate operators in a generalized particle method, enhancing the understanding of its numerical accuracy and convergence behavior.

Contribution

It introduces a new regularity framework and error estimates for operators in a generalized particle method, unifying analysis across various particle techniques.

Findings

Error estimates align with theoretical convergence rates.

Numerical results confirm the accuracy of the derived estimates.

The framework applies to multiple particle methods, including SPH.

Abstract

To facilitate the numerical analysis of particle methods, we derive truncation error estimates for the approximate operators in a generalized particle method. Here, a generalized particle method is defined as a meshfree numerical method that typically includes other conventional particle methods, such as smoothed particle hydrodynamics or moving particle semi-implicit methods. A new regularity of discrete parameters is proposed via two new indicators based on the Voronoi decomposition of the domain along with two hypotheses of reference weight functions. Then, truncation error estimates are derived for an interpolant, approximate gradient operator, and approximate Laplace operator in the generalized particle method. The convergence rates for these estimates are determined based on the frequency with which they appear in the regularity and hypotheses. Finally, the estimates are computed numerically and the results are shown to be in good agreement with the theoretical results.

Tables

Table 1. Table 1. Numerical and theoretical convergence rates of (a) the interpolant, (b) approximate gradient operator, and (c) approximate Laplace operator with regular orders m = 1 , 3 , 5 𝑚 1 3 5 m=1,3,5 . The numerical convergence rates were obtained for the cases of Δ ​ x = 2 − 11 Δ 𝑥 superscript 2 11 \Delta x=2^{-11} and 2 − 12 superscript 2 12 2^{-12} .

	$m=1$		$m=3$		$m=5$
	Numer.	Theor.	Numer.	Theor.	Numer.	Theor.
(a) Interpolant
($\Pi 1$)	-0.10	N/A	2.02	2	2.00	2
($\Pi 2$)	0.05	N/A	2.13	2	2.01	2
($\Pi 3$)	0.00	N/A	4.20	2	7.41	4
(b) Gradient
($\nabla 1$)	-0.05	N/A	2.11	2	2.02	2
($\nabla 2$)	-0.02	N/A	2.08	2	2.03	2
($\nabla 3$)	-0.06	N/A	3.56	2	7.69	4
(c) Laplacian
($\Delta 1$)	-1.14	N/A	2.91	1	2.05	2
($\Delta 2$)	-1.09	N/A	2.39	1	2.23	2
($\Delta 3$)	-1.02	N/A	2.76	1	7.50	3

Equations

$${\Vert v\Vert }_{C(\overline{S})}$$

$${\left|v\right|}_{C^{\ell }(\overline{S})}$$

$${\Vert v\Vert }_{C^{\ell }(\overline{S})}$$

$${\Omega }_H:=\left\{x\in {\mathbb{R}}^d|\exists y\in \Omega \text{mbox}s.t.|x-y|<H\right\}.$$

$${\mathcal{X}}_N:=\left\{x_i\in {\Omega }_H;i=1,2,\dots ,N,x_i\ne x_j(i\ne j)\right\},$$

$${\mathcal{V}}_N:=\left\{V_i\in {\mathbb{R}}^{+};i=1,2,\dots ,N,\sum _{i=1}^N{V}_i=\left|{\Omega }_H\right|\right\},$$

$$\mathcal{W}:=\left\{w\in C({\mathbb{R}}_0^{+});\text{mbox}supp(w)=[0,1],{\int }_{{\mathbb{R}}^d}w(|x|)\mathrm{d}x=1,\text{mbox}absolutelycontinuous\right\},$$

$$w_h(r):=\frac{1}{h^d}w\left(\frac{r}{h}\right).$$

$$\text{mbox}supp(w_h)=[0,h],{\int }_{{\mathbb{R}}^d}{w}_h(|x|)\mathrm{d}x=1,$$

$${\Pi }_{h}v(x)$$

$${\nabla }_{h}v(x)$$

$${\Delta }_{h}v(x)$$

$${\Lambda }_0(x,r)$$

$$\Lambda (x,r)$$

$${\sigma }_i:=\left\{x\in {\Omega }_H;|x_i-x|<|x_j-x|,\forall x_j\in {\mathcal{X}}_N(j\ne i)\right\},i=1,2,\dots ,N.$$

$$\Xi :=\left\{{\xi }_i\subset {\Omega }_H;\left|{\xi }_i\right|=V_i,\bigcup _{i=1}^N{\overline{\xi }}_i={\overline{\Omega }}_H,{\xi }_i\cap {\xi }_j=\varnothing (i\ne j)\right\}.$$

$$r_N:=\underset{i=1,2,\dots ,N}{\max }\underset{x\in {\sigma }_i}{\sup }|x_i-x|.$$

$$d_N:=\underset{\Xi }{\inf }{d}_{\Xi }$$

$$d_{\Xi }:=\underset{i=1,2,\dots ,N}{\max }\left\{\sum _{j=1}^N\frac{\left|{\sigma }_i\cap {\xi }_j\right|+\left|{\xi }_i\cap {\sigma }_j\right|}{\left|{\sigma }_i\right|}|x_i-x_j|\right\}.$$

$$h_N^m\ge c_0(r_N+d_N),\forall N\in \mathbb{N}.$$

$${\int }_{{\mathbb{R}}^d}{x}^{\alpha }w(|x|)\mathrm{d}x=0.$$

$$\max \left\{\underset{r\in (0,1)}{\sup }|w^{(k+1)}(r)|,\underset{r\in (0,1)}{\sup }\left|\frac{d}{dr}{w}^{(k)}(r)\right|\right\}<\infty ,$$

$$w^{(j)}(r):=\{\begin{array}{ll}\underset{s\downarrow 0}{\lim }\frac{w(s)}{s^j},& r=0,\\ \frac{w(r)}{r^j},& r>0.\end{array}$$

$${\Vert v-{\Pi }_{h}v\Vert }_{C(\overline{\Omega })}$$

$${\Vert \nabla v-{\nabla }_{h}v\Vert }_{C(\overline{\Omega })}\le c{h}^{\min \{m-1,n+1\}}{\Vert v\Vert }_{C^{n+2}({\overline{\Omega }}_H)},v\in C^{n+2}({\overline{\Omega }}_H),$$

$${\Vert \Delta v-{\Delta }_{h}v\Vert }_{C(\overline{\Omega })}\le c{h}^{\min \{m-2,n+1\}}{\Vert v\Vert }_{C^{n+3}({\overline{\Omega }}_H)},v\in C^{n+3}({\overline{\Omega }}_H).$$

$$I_{\alpha }(x):=\sum _{i\in {\Lambda }_0(x,h)}{V}_i(x_i-x{)}^{\alpha }{w}_h(|x_i-x|)-{\int }_{{\mathbb{R}}^d}{y}^{\alpha }{w}_h(|y|)\mathrm{d}y,x\in \overline{\Omega }.$$

$$I_{\alpha ,\ell }(x):=\sum _{i\in \Lambda (x,h)}{V}_i\frac{(x_i-x{)}^{\alpha }}{|x_i-x{|}^{\ell }}{w}_h(|x_i-x|)-{\int }_{{\mathbb{R}}^d}\frac{y^{\alpha }}{|y{|}^{\ell }}{w}_h(|y|)\mathrm{d}y,x\in \overline{\Omega }.$$

$$J_{\ell }(x):=\sum _{i\in {\Lambda }_0(x,h)}{V}_i|x_i-x{|}^{\ell }|w_h(|x_i-x|)|,x\in \overline{\Omega }.$$

$${\Vert v-{\Pi }_{h}v\Vert }_{C(\overline{\Omega })}\le c\left(\sum _{1\le |\alpha |\le n}{\Vert I_{\alpha }\Vert }_{C(\overline{\Omega })}+{\Vert J_{n+1}\Vert }_{C(\overline{\Omega })}\right){\Vert v\Vert }_{C^{n+1}({\overline{\Omega }}_H)},v\in C^{n+1}({\overline{\Omega }}_H),$$

$${\Vert \nabla v-{\nabla }_{h}v\Vert }_{C(\overline{\Omega })}\le c\left(\sum _{2\le |\alpha |\le n+2}{\Vert I_{\alpha ,2}\Vert }_{C(\overline{\Omega })}+{\Vert J_{n+1}\Vert }_{C(\overline{\Omega })}\right){\Vert v\Vert }_{C^{n+2}({\overline{\Omega }}_H)},v\in C^{n+2}({\overline{\Omega }}_H),$$

$${\Vert \Delta v-{\Delta }_{h}v\Vert }_{C(\overline{\Omega })}\le c\left(\sum _{1\le |\alpha |\le n+3}{\Vert I_{\alpha ,2}\Vert }_{C(\overline{\Omega })}+{\Vert J_{n+1}\Vert }_{C(\overline{\Omega })}\right){\Vert v\Vert }_{C^{n+3}({\overline{\Omega }}_H)},v\in C^{n+3}({\overline{\Omega }}_H).$$

Full text

1. Introduction

Particle methods, such as the smoothed particle hydrodynamics (SPH) [10, 18, 19] and moving particle semi-implicit (MPS) [15, 16, 29] methods, are numerical methods for solving partial differential equations that are based on points called particles distributed in a domain. In such methods, an interpolant and several approximate differential operators are defined in terms of linear combinations of weighted interactions between neighboring particles. When such methods are applied to partial differential equations, the equations are effectively discretized in space. As the discretization procedure does not require mesh generation in the domain, particle methods can be applied to moving boundary problems, such as the deformation and destruction of structures [5, 22] and flow problems associated with free surfaces [21, 23].

The accuracy of particle methods has been widely researched. From an engineering perspective, many studies have been conducted into the convergence of such methods in practical applications, such as Amicarelli [1, 2], Fulk [9], and Quinlan et al. [25]. On the other hand, few studies in the literature have presented numerical analyses of these methods from a mathematical perspective. In the 1980s, Mas-Gallic and Ravirt [20] and Ravirt [26] provided error estimates for particle methods when applied to parabolic and hyperbolic partial differential equations on unbounded domains. In the 2000s, Moussa and Via [4] and Moussa [3] provided error estimates of nonlinear conservation lows on bounded domains. In their work, the time integrations of the particle positions and volumes were obtained by solving the differential equations with respect to advection fields. However, as their method is only applicable to problems described by solvable differential equations, it cannot be used with other problems, such as those involving the Navier–Stokes equations.

Sometime later, Ishijima and Kimura [13] developed a truncation error estimate for an approximate gradient operator in MPS. By introducing a regularity for particle distributions based on an indicator called the equivolume partition radius, they determined the conditions that depend solely on the space distributions of the particles. However, a practical limitation is that the indicator cannot be computed.

In previous works, we established truncation error estimates for an interpolant, approximate gradient operator, and approximate Laplace operator of a generalized particle method in which the particle volumes were given as Voronoi volumes [11, 12]. A generalized particle method is a numerical method that typically includes conventional particle methods, such as SPH and MPS. In previous studies, we derived truncation error estimates by introducing a regularity using an indicator known as the covering radius, which is used in the numerical analysis of meshfree methods based on moving least-square methods and radial basis functions [17, 27, 30]. Although the formulations and conditions in those works are computable, they are difficult to deploy in practical computations as the computational costs associated with particle volumes based on Voronoi decomposition are high.

The focus of the current work was to analyze particle methods under more practical conditions by extending our results to cases with commonly used particle volumes. We also introduce another indicator of particle volumes, which we refer to as a Voronoi deviation, that represents the deviation between particle volumes and Voronoi volumes. Then, utilizing the Voronoi deviation, we extend the regularity and introduce two hypotheses of reference weight functions. Using the regularity and hypotheses, we derive truncation error estimates of the interpolant, approximate gradient operator, and approximate gradient operator of the generalized particle method. Finally, we numerically analyze our estimates and compare the results to those from the theory.

The remainder of this paper is organized as follows. The interpolant and approximate operators of the generalized particle method are introduced in Section 2. A regularity describing the family of discrete parameters is discussed in Section 3, after which we propose our primary theorem with respect to the truncation error estimates and provide some corollaries. Then, the primary theorem is proven in Section 4, numerical results are detailed in Section 5, and some concluding remarks are outlined in Section 6.

In the remainder of this section, we describe some notation and define some relevant function spaces. Let $\mathbb{R}^{+}$, $\mathbb{R}^{+}_{0}$, and $\mathbb{N}_{0}$ be the set of positive real numbers, the set of nonnegative real numbers, and the set of nonnegative integers, respectively. Let $d$ be the dimension of a space. Let $\mathbb{A}^{d}$ be the set of all $d$-dimensional multi-indices. If there is no ambiguity, the symbol $|\cdot|$ is used to denote the following: $|x|$ denotes the Euclidean norm for a vector $x\in\mathbb{R}^{d}$; $|S|$ denotes the volume of $S$ for open set $S\subset\mathbb{R}^{d}$; $|\alpha|$ denotes $|\alpha|\mathrel{\mathop{:}}=\alpha_{1}+\alpha_{2}+\dots+\alpha_{d}$ for $\alpha\in\mathbb{A}^{d}$. For $S\subset\mathbb{R}^{d}$, let $\textrm{diam}(S)$ be $\textrm{diam}(S)\mathrel{\mathop{:}}=\sup\left\{|x-y|;\,x,y\in S\right\}$. For $S\subset\mathbb{R}^{d}$, let $C(\overline{S})$ be the space of real continuous functions defined in $\overline{S}$ with the norm $\left\|\,\cdot\,\right\|_{C(\overline{S})}$ defined as

[TABLE]

For $S\subset\mathbb{R}^{d}$ and $\ell\in\mathbb{N}$, let $C^{\ell}(\overline{S})$ be the space of functions in $C(\overline{S})$ with derivatives up to the $\ell$th order with its seminorm $\left|\,\cdot\,\right|_{C^{\ell}(\overline{S})}$ and norm $\left\|\,\cdot\,\right\|_{C^{\ell}(\overline{S})}$ defined as

[TABLE]

respectively. Here $D^{\alpha}v\mathrel{\mathop{:}}=\partial_{1}^{\alpha_{1}}\partial_{2}^{\alpha_{2}}\dots\partial_{d}^{\alpha_{d}}v$ with multi-index $\alpha=(\alpha_{1},\alpha_{2},\dots,\alpha_{d})$.

2. Approximate operators in a generalized particle method

Let $\Omega$ be a bounded domain in $\mathbb{R}^{d}$. Let $H$ be a fixed positive number. For $\Omega$ and $H$, we define extended domain $\Omega_{H}$ as

[TABLE]

For $N\in\mathbb{N}$, we define a particle distribution $\mathcal{X}_{N}$ and particle volume set $\mathcal{V}_{N}$ as

[TABLE]

[TABLE]

respectively. We refer to $x_{i}\in\mathcal{X}_{N}$ and $V_{i}\in\mathcal{V}_{N}$ as a particle and particle volume, respectively. An example of the particle distribution $\mathcal{X}_{N}$ in $\Omega_{H}\,(\subset\mathbb{R}^{2})$ is shown in Figure 1.

We define an admissible reference weight function set $\mathcal{W}$ as

[TABLE]

we refer to $w\in\mathcal{W}$ as a reference weight function, and we define the influence radius $h_{N}\in\mathbb{R}$ as satisfying $0<h_{N}<H$ and $h_{N}\rightarrow 0\,(N\rightarrow\infty)$. If there is no ambiguity, we denote $h_{N}$ as $h$. For reference weight function $w$ and influence radius $h$, we define the weight function $w_{h}\in C(\mathbb{R}^{+}_{0})$ as

[TABLE]

Note that the weight function $w_{h}$ satisfies

[TABLE]

and is absolutely continuous.

For $v\in C(\overline{\Omega}_{H})$, we define interpolant $\Pi_{h}{}$, approximate gradient operator $\nabla_{h}{}$, and approximate Laplace operator $\Delta_{h}{}$ as

[TABLE]

respectively. Here, for $x\subset\mathbb{R}^{d}$ and $r\in\mathbb{R}^{+}$, $\Lambda_{0}(x,r)$ and $\Lambda(x,r)$ are index sets of particles defined as

[TABLE]

respectively.

As discussed later in Appendix A, the approximate operators (10), (11), and (12) indicate a wider class of approximate operators of particle methods that those in the SPH and MPS methods. Therefore, we refer to the approximate operators (10), (11), and (12) as generalized approximate operators and to a particle method that uses them as a generalized particle method.

3. Truncation error estimates of approximate operators

We first introduce a regularity of discrete parameters. Let $\{\sigma_{i}\}$ be the Voronoi decomposition of $\Omega_{H}$ associated with the particle distribution $\mathcal{X}_{N}$, where $\sigma_{i}$ is the Voronoi region defined as

[TABLE]

We define a particle volume decomposition $\Xi$ as a decomposition of $\Omega_{H}$ associated with particle volume set $\mathcal{V}_{N}$ defined as

[TABLE]

An example of the Voronoi decomposition of $\Omega_{H}$ associated with the particle distribution $\mathcal{X}_{N}$ is shown in Figure 2.

We define a covering radius $r_{N}$ for particle distribution $\mathcal{X}_{N}$ as

[TABLE]

Moreover, we define a Voronoi deviation $d_{N}$ for the particle distribution $\mathcal{X}_{N}$ and the particle volume set $\mathcal{V}_{N}$ as

[TABLE]

with

[TABLE]

Then, we define a regularity for a family consisting of a particle distribution $\mathcal{X}_{N}$, particle volume set $\mathcal{V}_{N}$, and influence radius $h$ as follows:

Definition 1**.**

Family $\{(\mathcal{X}_{N},\mathcal{V}_{N},h_{N})\}_{N\rightarrow\infty}$ is said to be regular with order $m\,(m\geq 1)$ if there exists a positive constant $c_{0}$ such that

[TABLE]

Remark 2**.**

As shown in Figure 3, the covering radius $r_{N}$ becomes large in the case of a particle distribution with both dense and sparse regions. Therefore, the covering radius $r_{N}$ can be considered as an indicator representing the uniformness of particle distribution $\mathcal{X}_{N}$.

Remark 3**.**

A Voronoi deviation $d_{N}$ equals zero if and only if the particle volumes are given as the Voronoi volume ($V_{i}=\left|\sigma_{i}\right|$). Moreover, the Voronoi deviation $d_{N}$ becomes large if the particle volumes are given as values far from the Voronoi volumes. Therefore, the Voronoi deviation $d_{N}$ can be regarded as an indicator of the deviation between the particle volume set and the Voronoi volume set.

Remark 4**.**

For a given family $\{(\mathcal{X}_{N},\mathcal{V}_{N},h_{N})\}_{N\rightarrow\infty}$ and given constant $m$ $(m\geq 1)$, it is possible to determine whether or not the family is regular with order $m$ as the covering radius $r_{N}$ and Voronoi deviation $d_{N}$ are absolutely computable, as shown in Appendix B.

Next, we introduce two hypotheses of reference weight function $w$:

Hypothesis 5**.**

For $n\in\mathbb{N}$, the reference weight function $w$ satisfies for all $\alpha\in\mathbb{A}^{d}$ with $1\leq|\alpha|\leq n$,

[TABLE]

Hypothesis 6**.**

For $k\in\mathbb{N}_{0}$, the reference weight function $w$ satisfies

[TABLE]

where for $j\in\mathbb{N}_{0}$, $w^{(j)}(r):(0,\infty)\rightarrow\mathbb{R}$ is defined as

[TABLE]

Remark 7**.**

All reference functions $w\in\mathcal{W}$ satisfy Hypothesis 5 with $n=1$. Moreover, for all $n\in\mathbb{N}$ and $k\in\mathbb{N}$, reference weight functions satisfying Hypothesis 5 with $n$ and Hypothesis 6 with $k$ can be constructed as shown in Appendix C.

We now state a theorem that defines truncation error estimates of approximate operators in the generalized particle method with a continuous norm:

Theorem 8**.**

Suppose that family $\{(\mathcal{X}_{N},\mathcal{V}_{N},h_{N})\}_{N\rightarrow\infty}$ is regular with order $m\,(m\geq 1)$ and that reference weight function $w$ satisfies Hypothesis 5 with $n$. Then, there exists a positive constant $c$ independent of $N$ such that

[TABLE]

In addition, if $w\in\mathcal{W}$ satisfies Hypothesis 6 with $k=0$, then we have

[TABLE]

and if $w\in\mathcal{W}$ satisfies Hypothesis 6 with $k=1$, then we have

[TABLE]

The proof of Theorem 8 is presented in the next section. As shown in the corollaries in Appendix A, the approximate operators commonly used in the SPH and MPS methods are valid for Theorem 8 under appropriate settings.

4. Proof of truncation error estimates

The following notation will be used in the subsequent proof of Theorem 8. Hereafter, let $c$ be a generic positive constant independent of $N$. For $\alpha\in\mathbb{A}^{d}$, set $I_{\alpha}$ as

[TABLE]

For $\alpha\in\mathbb{A}^{d}$ and $\ell\in\mathbb{N}$, set $I_{\alpha,\ell}$ as

[TABLE]

For $\ell\in\mathbb{N}$, set $J_{\ell}$ as

[TABLE]

We now present the following lemma.

Lemma 9**.**

Suppose that $w\in\mathcal{W}$ satisfies Hypothesis 5 with $n$. Then, there exists a positive constant $c$ independent of $N$ such that

[TABLE]

[TABLE]

[TABLE]

Proof.

First, we prove (30). We fix $x\in\overline{\Omega}$. Then, let $B(x,r)$ be the open ball in $\mathbb{R}^{d}$ with center $x$ and radius $r$, i.e.,

[TABLE]

From $h<H$, we have $B(x,h)\subset\Omega_{H}$. Then, for all $v\in C^{\ell+1}(\overline{\Omega}_{H})$ $(\ell\in\mathbb{N})$ and $y\in B(x,h)$, we obtain the Taylor expansion of $v$ as

[TABLE]

From (10) and (34) with $\ell=n+1$, we have

[TABLE]

Moreover, by Hypothesis 5, we have

[TABLE]

Because

[TABLE]

we have

[TABLE]

Moreover, we have

[TABLE]

Therefore, from (37), (39), and (40), we obtain (30).

Next, we prove (31). From (11) and (34) with $\ell=n+2$, we have

[TABLE]

Because for $\beta\in\mathbb{A}^{d}$ with $|\beta|=2$,

[TABLE]

we have

[TABLE]

Hypothesis 5 with $n$ yields

[TABLE]

From (43) and (44), we have

[TABLE]

From (38), we have

[TABLE]

Moreover, we have

[TABLE]

Therefore, from (47), (48), and (49), we obtain (31).

Finally, we prove (32). From (12) and (34) with $\ell=n+3$, we have

[TABLE]

From (42), we have

[TABLE]

Hypothesis 5 with $n$ yields

[TABLE]

Therefore, we have

[TABLE]

From (38), we have

[TABLE]

Moreover, we have

[TABLE]

Therefore, from (54), (55), and (56), we obtain (32). ∎

Next, we show estimates of $I_{\alpha}$, $I_{\alpha,\ell}$, and $J_{\ell}$.

Lemma 10**.**

There exists a positive constant $c$ independent of $N$ such that

[TABLE]

Proof.

We arbitrarily fix $x\in\overline{\Omega}$, $\alpha\in\mathbb{A}^{d}$, and particle volume decomposition $\Xi=\{\xi_{i}\mid i=1,2,\dots,N\}$ and split $I_{\alpha}$ into

[TABLE]

with

[TABLE]

Then, we estimate $E_{1}$, $E_{2}$, and $E_{3}$.

First, we estimate $E_{1}$. Because

[TABLE]

we can rewrite $E_{1}$ as

[TABLE]

From

[TABLE]

we obtain

[TABLE]

From

[TABLE]

we have

[TABLE]

Because $w_{h}$ is absolutely continuous, we have

[TABLE]

for all $y,z\in\mathbb{R}^{d}$. Moreover, we have

[TABLE]

From (67), (68), and (69), we have

[TABLE]

Therefore, from (65) and (70), we obtain

[TABLE]

Next, we estimate $E_{2}$. Because $\mbox{supp}(w_{h})=[0,h]$ and $\sigma_{j}\subset B(x_{j},r_{N})$, we have

[TABLE]

From (72), we have

[TABLE]

Moreover, from (68) and (69), we have

[TABLE]

Therefore, from (64) and (74), we obtain

[TABLE]

Finally, we estimate $E_{3}$. Because

[TABLE]

we can rewrite $E_{3}$ as

[TABLE]

Because $E_{3}=0$ when $|\alpha|=0$, we estimate when $|\alpha|\geq 1$. Let $\beta_{k}\,(k=1,2,\dots,|\alpha|)$ be $d$-dimensional multi-indices with satisfying

[TABLE]

Then, we have, for all $y,z\in\mathbb{R}^{d}$,

[TABLE]

From (64) and (79), we obtain

[TABLE]

By $\mbox{supp}(w_{h})=[0,h]$ and $\sigma_{j}\subset B(x_{j},r_{N})$, if $j\not\in\Lambda_{0}(x,h+r_{N})$, then

[TABLE]

Moreover, from $w\in\mathcal{W}\subset C(\mathbb{R}^{+}_{0})$, we have

[TABLE]

From (69), (81), and (82), we have

[TABLE]

Therefore, from (80), (83), and $h\leq H$, we obtain

[TABLE]

From the estimates of $E_{1}$, $E_{2}$, and $E_{3}$, we obtain

[TABLE]

Because $\Xi$ is arbitrary, we establish (57). ∎

Lemma 11**.**

Suppose that reference weight function $w$ satisfies Hypothesis 6 with $k$. Then, there exists a positive constant $c$ independent of $N$ such that for all $\alpha\in\mathbb{A}^{d}$ and $\ell\in\mathbb{N}$ with $1\leq\ell-k\leq|\alpha|$,

[TABLE]

Proof.

We arbitrarily fix $x\in\overline{\Omega}$, $\alpha\in\mathbb{A}^{d}$, particle volume decomposition $\Xi=\{\xi_{i}\mid i=1,2,\dots,N\}$, and $\ell\in\mathbb{N}$ with $1\leq\ell-k\leq|\alpha|$ and split $I_{\alpha,\ell}$ into

[TABLE]

with

[TABLE]

Then, we estimate $E_{4}$, $E_{5}$, and $E_{6}$.

First, we estimate $E_{4}$ and set $w^{(k)}$ as (23) and $w^{(k)}_{h}$ as

[TABLE]

Then, from (62), we can rewrite $E_{4}$ as

[TABLE]

Because

[TABLE]

we obtain

[TABLE]

From $\mbox{supp}(w^{(k)}_{h})=[0,h]$, we have

[TABLE]

Thus, we obtain

[TABLE]

Using an argument similar to (68), if $w$ satisfies Hypothesis 6 with $k$, then for all $y,z\in\mathbb{R}^{d}$,

[TABLE]

From (96) and (97), we obtain

[TABLE]

Next, we estimate $E_{5}$. By using $w^{(k)}_{h}$, we can rewrite $E_{5}$ as

[TABLE]

From (93), we obtain

[TABLE]

By $\mbox{supp}(w^{(k)}_{h})=[0,h]$ and $\sigma_{j}\subset B(x_{j},r_{N})$, we have

[TABLE]

From (97) and (101), we obtain

[TABLE]

Finally, we estimate $E_{6}$. Using $w^{(k)}_{h}$, we can rewrite $E_{6}$ as

[TABLE]

where $\xi^{\ast}_{i}(x)$ is

[TABLE]

For $\alpha\in\mathbb{A}^{d}$, let $\beta_{j}\,(j=1,2,\dots,|\alpha|)$ be $d$-dimensional multi-indices satisfying

[TABLE]

Let $\beta_{j}^{\ast}\,(j=0,1,\dots,|\alpha|)$ be $d$-dimensional multi-indices defined as

[TABLE]

For all $y,z\in\mathbb{R}^{d}\setminus\{0\}$, when $|\alpha|=\ell-k$, we have

[TABLE]

Moreover, from (79) and (107), when $|\alpha|>\ell-k$, we have

[TABLE]

Therefore, when $|\alpha|\geq\ell-k$, we have for all $y\in\Omega_{H}\setminus\{x\}$ and $i\in\Lambda(x,\infty)$,

[TABLE]

From (109), we obtain

[TABLE]

Because $|\alpha|\geq\ell-k$, we have

[TABLE]

Therefore, we have

[TABLE]

Because for all $y\in\Omega_{H}$,

[TABLE]

by the same procedure as (83), we have

[TABLE]

Therefore, we obtain

[TABLE]

From the estimates of $E_{4}$, $E_{5}$, and $E_{6}$, we obtain

[TABLE]

Because $\Xi$ is arbitrary, we establish (86). ∎

Lemma 12**.**

There exists a positive constant $c$ independent of $N$ such that

[TABLE]

Proof.

We arbitrarily fix $x\in\overline{\Omega}$ and particle volume decomposition $\Xi=\{\xi_{i}\mid i=1,2,\dots,N\}$, and split $J_{\ell}$ into

[TABLE]

with

[TABLE]

Then, we estimate $E_{7}$, $E_{8}$, $E_{9}$, and $E_{10}$.

From (62), we can rewrite $E_{7}$ as

[TABLE]

For all $y\in\Omega_{H}$, we have

[TABLE]

From (70) and (124), we obtain

[TABLE]

From (74) and (124), we obtain

[TABLE]

For all $x_{i}\in\mathcal{X}_{N}$ and $y\in\Omega_{H}$, we have

[TABLE]

From (83), (127), and $h<H$, we obtain

[TABLE]

From (8), we obtain

[TABLE]

From the estimates of $E_{7}$, $E_{8}$, $E_{9}$, and $E_{10}$, we obtain

[TABLE]

Because $\Xi$ is arbitrary, we establish (117). ∎

Using the lemmas defined above, we now prove Theorem 8.

Theorem 8.

By Lemmas 9, 10, and 12, we have for all $v\in C^{n+1}(\overline{\Omega}_{H})$

[TABLE]

Moreover, by Lemmas 9, 11, and 12, when $w$ satisfies Hypothesis 6 with $k=0$, we have for all $v\in C^{n+2}(\overline{\Omega}_{H})$

[TABLE]

and when $w$ satisfies Hypothesis 6 with $k=1$ for all $v\in C^{n+3}(\overline{\Omega}_{H})$,

[TABLE]

Because the family $\{(\mathcal{X}_{N},\mathcal{V}_{N},h_{N})\}_{N\rightarrow\infty}$ is regular, by applying (20) to (131), (132), and (133), we obtain (24), (25), and (26), respectively. We now conclude the proof of Theorem 8. ∎

5. Numerical results

Set $\Omega=(0,1)^{2}$ and $H=0.1$. Then, $\Omega_{H}=(-0.1,1.1)^{2}$. We now compute the truncation errors of $v:\Omega_{H}\rightarrow\mathbb{R}$, which are defined as $v(x,y)=\sin(2\pi(x+y))$. Particle distribution $\mathcal{X}_{N}$ is set as

[TABLE]

where $\Delta x$ is taken by $2^{-5},2^{-6},\dots,2^{-12}$ and $\eta_{ij}^{(k)}\,(i,j\in\mathbb{Z},k=1,2)$ are random numbers satisfying $|\eta_{ij}^{(k)}|<1/4$. Particle distribution $\mathcal{X}_{N}$ with $\Delta x=2^{-5}$ is shown in Figure 4.

Particle volume set $\mathcal{V}_{N}$ is defined as

[TABLE]

For $m=1,3,5$, the influence radius $h_{N}$ is set as

[TABLE]

Note that if $\Delta x=2^{-5}$, then $h=2.6\times 2^{-5}$ for all $m$. Using the discrete parameters above, the covering radius $r_{N}$ satisfies $r_{N}\leq\sqrt{2}(1+1/4)\Delta x/2$. Moreover, the Voronoi deviation $d_{N}$ satisfies $d_{N}\leq 64(1+\sqrt{2})\Delta x/\pi$. Therefore, the family $\{(\mathcal{X}_{N},\mathcal{V}_{N},h_{N})\}$ is regular with order $m$.

For the interpolant, we consider the following three cases of reference weight functions:

[TABLE]

$({\rm\Pi 1})$ is the lowest order polynomial function belonging to $\mathcal{W}$. $({\rm\Pi 2})$ is the cubic B-spline commonly used in SPH and belonging to $\mathcal{W}$. $({\rm\Pi 3})$ is the lowest order polynomial function belonging to $\mathcal{W}$ that satisfies Hypothesis 5 with $n=3$.

For the approximate gradient operator, we consider the following three cases of reference weight functions:

[TABLE]

$({\rm\nabla 1})$ is the lowest order polynomial function belonging to $\mathcal{W}$ and Hypothesis 6 with $k=0$. $({\rm\nabla 2})$ is chosen so that approximate gradient operator (11) with $({\rm\nabla 2})$ coincides with that in SPH with the cubic B-spline (see Appendix A). $({\rm\nabla 3})$ is the lowest order polynomial function belonging to $\mathcal{W}$, Hypothesis 5 with $n=3$, and Hypothesis 6 with $k=0$.

For the approximate Laplace operator, we consider the following three cases of reference weight functions:

[TABLE]

$({\rm\Delta 1})$ is the lowest order polynomial function belonging to $\mathcal{W}$ and Hypothesis 6 with $k=1$. $({\rm\Delta 2})$ is chosen so that approximate Laplace operator (12) with $({\rm\Delta 2})$ coincides with that in SPH with the cubic B-spline (see Appendix A). $({\rm\Delta 3})$ is the lowest order polynomial function belonging to $\mathcal{W}$, Hypothesis 5 with $n=3$, and Hypothesis 6 with $k=1$.

The above settings were used in the computation of the following relative errors

[TABLE]

Here, the discrete norm $\left\|\cdot\right\|_{\ell^{\infty}(\Omega)}$ is defined as

[TABLE]

Figure 5 shows graphs of the relative errors of (a) interpolant $\Pi_{h}{}$, (b) approximate gradient operator $\nabla_{h}{}$, and (c) approximate Laplace operator $\Delta_{h}{}$ versus the influence radius $h_{N}$ with regular orders $m=1,3,5$. In Figure 5, the slopes of the triangles show the theoretical convergence rates obtained via Theorem 8. Table 1 lists the numerical and theoretical convergence rates obtained from the cases of $\Delta x=2^{-11}$ and $2^{-12}$, where the theoretical convergence rates correspond to Theorem 8. In the case of $m=1$, as the settings could not be applied to Theorem 8, only numerical results without convergence were obtained. In contrast, the settings in cases $m=3$ and $5$ could be applied Theorem 8; thus, the numerical results with convergence were obtained. Moreover, the approximate operators with reference weight functions satisfying Hypothesis 5 with $n=3$ became higher convergence orders in the cases where $m=5$ as per Theorem 8.

6. Conclusions

We analyzed truncation errors in a generalized particle method, which is a wider class of particle methods that includes commonly used methods such as SPH and MPS. In our analysis, we introduced two indicators: the first was the covering radius, which represents the maximum radius of the Voronoi region associated with the particle distribution, while the second was the Voronoi deviation, which indicates the deviation between particle volumes and Voronoi volumes. With the covering radius and Voronoi deviation, we introduced a regularity of a family of discrete parameters, which includes the particle distribution, particle volume set, and influence radius associated with the number of particles. Moreover, we introduced two hypotheses of reference weight functions. With the regularity and hypotheses of reference weight functions, we established truncation error estimates for the continuous norm. The convergence rates are dependent on the regular order and order of the reference weight functions appearing in a hypothesis. Moreover, as it was possible to validate the conditions by calculation, we showed the numerical convergence orders were in good agreement with the theoretical ones.

In a forthcoming paper, we plan to establish error estimates of the generalized particle method for the Poisson and heat equations.

Acknowledgments

We would like to thank Drs. Daisuke Tagami and Hayato Waki for their helpful comments over the course of this work. This research was supported by JSPS KAKENHI Grant Number 17K17585 and the JSPS A3 Foresight Program.

Appendix A Description of conventional particle methods by the generalized particle method

This appendix provides a description of conventional particle methods, such as the smoothed particle hydrodynamics (SPH) [18, 24] and the moving particle semi-implicit (MPS) [16], in the context of the generalized particle method. In SPH, by using the reference weight function $w^{\mbox{\rm\tiny SPH}}\in\mathcal{W}$ and parameters $m_{i},\rho_{i}\in\mathbb{R}^{+}\,(i=1,2,\dots,N)$, for $v\in C(\overline{\Omega}_{H})$, the approximate operators are defined as

[TABLE]

By setting $w=w^{\mbox{\rm\tiny SPH}}$ and $\mathcal{V}_{N}=\{V_{i}=m_{i}/\rho_{i};\,i=1,2,\dots,N\}$, the generalized interpolant (10) coincides with (139). Moreover, because

[TABLE]

by setting

[TABLE]

and $\mathcal{V}_{N}=\{V_{i}=m_{i}/\rho_{i};\,i=1,2,\dots,N\}$, (11) and (12) coincide with (140) and (141), respectively.

From Theorem 8, we obtain the following corollary that is a truncation error estimate of approximate operators (140) and (141).

Corollary 13**.**

Suppose that parameters $\rho_{i},m_{i}$ satisfy

[TABLE]

and that $\{(\mathcal{X}_{N},\mathcal{V}_{N},h_{N})\}_{N\rightarrow\infty}$ is regular with order $m$, where $\mathcal{V}_{N}=\{\rho_{i}/m_{i};\,i=1,2,\dots,N\}$. Moreover, suppose that $w^{\mbox{\rm\tiny SPH}}$ satisfies the following conditions;**

[TABLE]

Then, there exists a positive constant $c$ independent of $N$ such that

[TABLE]

Remark 14**.**

Note that representative reference weight functions employed in SPH, such as the cubic B-spline, quintic B-spline, and Wendland function (5-order positive definite function) [8, 18], satisfy (145).

In MPS [14], by using reference weight function $w^{\mbox{\rm\tiny MPS}}\in\mathcal{W}$ and parameters $\widehat{n},\widehat{\lambda}\in\mathbb{R}^{+}$ for $v\in C^{1}(\overline{\Omega}_{H})$, approximate differential operators can be defined as

[TABLE]

Note that an interpolant is not defined in MPS. By setting $w=w^{\mbox{\rm\tiny MPS}}$ and $\mathcal{V}_{N}=\{V_{i}=\widehat{n}^{-1};\,i=1,2,\dots,N\}$, the approximate gradient operator (11) coincides with (149). Moreover, by setting $w(r)=\widehat{\lambda}^{-1}r^{2}w^{\mbox{\rm\tiny MPS}}(r)$ and $\mathcal{V}_{N}=\{V_{i}=\widehat{n}^{-1};\,i=1,2,\dots,N\}$, approximate Laplace operator (12) coincides with (150).

Corollary 15**.**

Suppose that

[TABLE]

Moreover, suppose that $\{(\mathcal{X}_{N},\mathcal{V}_{N},h_{N})\}_{N\rightarrow\infty}$ is regular with order $m$, where $\mathcal{V}_{N}=\{V_{i}=\widehat{n}^{-1};\,i=1,2,\dots,N\}$. Then, there exists a positive constant $c$ independent of $N$ such that

[TABLE]

Furthermore, when $w^{\mbox{\rm\tiny MPS}}$ satisfies Hypothesis 6 with $k=0$,

[TABLE]

Remark 16**.**

Note that the reference weight function, which is commonly used in the MPS and defined as

[TABLE]

does not satisfy $w^{\mbox{\rm\tiny MPS}}\in\mathcal{W}$. In contrast, the continuous reference weight function as introduced in [28] satisfies $w^{\mbox{\rm\tiny MPS}}\in\mathcal{W}$. However, as far as we know, no reference weight functions that also satisfy Hypothesis 6 with $k=0$ are proposed in MPS.

Appendix B Computational procedure of the indicators

This appendix introduces the procedures for computing the indicators introduced in this paper, namely, the covering radius (17) and Voronoi deviation (18).

The covering radius $r_{N}$ can be computed as follows. As per the methods used to construct Voronoi decompositions, such as the increment method [6], we first draw the boundaries of the Voronoi region in $\Omega_{H}$. Next, for each particle, we compute the maximum distance from particle $x_{i}$ to the boundary of its Voronoi region $\sigma_{i}$ (i.e., $\max_{y\in\overline{\sigma_{i}}}|x_{i}-y|$). Finally, we obtain the covering radius $r_{N}$ by computing the maximum of these distances.

Next, we consider the Voronoi deviation $d_{N}$. Let $\zeta\in\mathbb{R}^{3N}$ be

[TABLE]

Using parameters $q,s_{i},a_{ij}\in\mathbb{R}^{+}$ $(i,j=1,2,\dots,N)$, we set $z\in\mathbb{R}^{N^{2}+N+1}$ as

[TABLE]

Moreover, we set $M\in\mathbb{R}^{3N\times(N^{2}+N+1)}$ so that equation $Mz=\zeta$ represents

[TABLE]

and

[TABLE]

Then, by considering $a_{ij}$ to be $\left|\sigma_{i}\cap\xi_{j}\right|$, we find that the minimum value of $q$ with condition $Mz=\zeta$ coincides with the Voronoi deviation $d_{N}$. We therefore consider the linear problem:

[TABLE]

Here, $b\mathrel{\mathop{:}}=(0,0,\dots,0,1)^{T}\in\mathbb{R}^{N^{2}+N+1}$. The solution $b^{T}z$ of (159) is equivalent to the Voronoi deviation $d_{N}$. Because $Mz=\zeta$ is unique for $(\mathcal{X}_{N},\mathcal{V}_{N},h_{N})$, the linear problem is computable via numerical methods for linear programming problems, such as the simplex method [7].

Appendix C Construction of reference weight functions

For all $n\in\mathbb{N}$ $(n\geq 2)$, it is possible to construct a reference weight function satisfying Hypothesis 5 with $n$ as the condition of Hypothesis 5 can be rewritten to include a finite number of conditions

[TABLE]

Here, the Gauss symbol $\lfloor a\rfloor$ denotes the largest integer that is less than or equal to $a$. For example, function $w$ is set as the $p$th polynomial function:

[TABLE]

Then, if coefficients $a_{\ell}$ satisfy the linear equations

[TABLE]

then $w$ satisfies $w\in\mathcal{W}$ and Hypothesis 5 with $n$. Therefore, to construct reference functions with Hypothesis 5 with $n$ represented by polynomial functions, the degree of the polynomial functions must be at least $\lfloor n/2\rfloor+2$.

Moreover, for all $k\in\mathbb{N}_{0}$, reference weight functions satisfying Hypothesis 6 with $k$ can be constructed based on the following proposition.

Proposition 17**.**

Assume that reference the weight function $w$ defined in $\mathbb{R}^{+}_{0}$ satisfies $w\in C^{1}(\mathbb{R}^{+}_{0})$ and is represented by a polynomial function in $[0,s]$ for $s\in(0,1]$. Let $p_{0}$ be the minimum degree of $w$ in $[0,s]$. Then, if $p_{0}-k\geq 1$, $w$ satisfies Hypothesis 6 with $k$.

Proof.

From the assumption, $w$ can be represented by

[TABLE]

where $p\in\mathbb{N}$ and $a_{\ell}\in\mathbb{R}$ $(\ell=p_{0},p_{0}+1,\dots,p)$. Set $w^{(k)}$ as (23). Since

[TABLE]

and

[TABLE]

if $p_{0}-k\geq 1$, we have $w$ satisfies Hypothesis 6 with $k$. ∎

This proposition means that the regularity of the reference functions around zero is important.