# A Numerical Convolution Representation of Potential for a Disk Surface Density in 3D

Authors: Chien-Chang Yen
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

This study presents a numerical convolution representation of a potential induced from a disk of surface density, which has often been investigated in the center region of galaxies. The advantage of this representation is to release the softening length for theN-body method and artificial boundary conditions for the spectral methods. With the help of fast Fourier transform, the computational complexity is only O
(
M
2
(
log
2
⁡
M
)
), where M is the number of zones in one dimension. Numerical results show an almost second order of accuracy on a Cartesian coordinate system. A comparison study also demonstrates that  this method can calculate the potential for disk surface density based on the uniform grids.


---

## Body

## 1. Introduction

Let us consider the Poisson equation with singular disk surface density [1–3]
(1)
Δ
Φ
(
x
)
=
4
π
G
ρ
(
x
)
=
4
π
G
σ
(
x
,
y
)
δ
(
z
)
,
in
R
3
,
wherex
=
(
x
,
y
,
z
),G,δ, and
(2)
ρ
(
x
)
=
σ
(
x
,
y
)
δ
(
z
)
,
are position, gravitational constant, Dirac symbol, and density, respectively. Hereafter, we setG
=
1without loss of generality.One possible method of solving (1) is to discretize the partial differential equation (1) using the finite difference method [4–6]. This discretization is
(3)
-
-
Φ
i
+
1
,
j
,
k
+
2
Φ
i
,
j
,
k
-
Φ
i
-
1
,
j
,
k
(
Δ
x
)
2
,
-
-
Φ
i
,
j
+
1
,
k
+
2
Φ
i
,
j
,
k
-
Φ
i
,
j
-
1
,
k
(
Δ
y
)
2
,
-
-
Φ
i
,
j
,
k
+
1
+
2
Φ
i
,
j
,
k
-
Φ
i
,
j
,
k
-
1
(
Δ
z
)
2
=
f
i
,
j
,
k
,
whereΦ
i
,
j
,
k
=
Φ
(
x
i
,
y
j
,
z
k
),f
i
,
j
,
0
=
σ
i
,
j, andf
i
,
j
,
k
=
0fork
≠
0, based on the uniform mesh grids(
x
i
,
y
j
,
z
k
)with the mesh sizesΔ
x,Δ
y, andΔ
zinx,y, andzdirections, respectively. In this approach, the values of the potential on imposed boundaries should be calculated and a fully 3D calculation must be performed. Any numerical solutions of this partial differential problem involveO
(
M
3
)unknowns, whereMis the number of zones in one dimension. The linear complexity of this approach is at leastO
(
M
3
), for example, [7]. It seems unfeasible that a thin gaseous disk problem can directly solve the partial differential equation.The other approach for solving (1) uses the integral form. The solution to (1) can be represented in terms of the fundamental solution,(
1
/
4
π
)
K
(
x
), where
(4)
K
(
x
)
=
1
x
2
+
y
2
+
z
2
,
as
(5)
Φ
(
x
,
y
,
z
)
=
-
∬
K
(
x
¯
-
x
,
y
¯
-
y
,
-
z
)
σ
(
x
¯
,
y
¯
)
d
x
¯
d
y
¯
.
The difficulties encountered in the numerical approach for solving (1) or (5) are related to the extent of the domainR
3. The spectral method uses trigonometric bases functions and the artificial periodic boundary conditions [8]. However, the potential is not periodic in reality. A direct calculation using theN-body method [1] has the total amount of complexity O
(
M
4
)based on the numberO
(
M
2
)of mesh zones. For theN-body method based on a uniform grid and using the FFT technique, the numerical complexity can be reduced fromO
(
M
4
)toO
(
M
2
(
log
⁡
2
(
M
)
)
2
). In other words, this produces a fast algorithm of linear complexity. Unfortunately, this method has a drawback, which is to introduce the softening parameter to overcome the singular integration
(6)
∬
D
K
(
x
,
y
,
0
)
d
x
d
y
,
whereDcontains the origin.The proposed method does not require the artificial periodic boundary condition and softening parameter and further improves the accuracy of theN-body method. This study uses the uniform and logarithmic grid generations in Cartesian and polar coordinates, respectively. The other treatment is linear approximation for the surface density on the cell. This ensures that the proposed method is of second order in Cartesian coordinates. This study uses two examples, aD
2disk [9] and its variationD
2,2, to show numerically that the method achieves second-order accuracy. A comparison study also demonstrates that the proposed method is suitable for potentials on a thin disk.Before the end of the section, another approach for an improvement ofN-body method is the fast multipole methods (FMM) [1, 10] which decomposes the region into local interactions and far-field interactions for every position. The FMM is efficient for fully three-dimensional calculations under the assumption that the far field is smooth. The FFM also relies on the ability to manipulate local multipole expansions for every box in the tree hierarchy [11–13]. The proposed presentation based on the uniform grids is without the assumption of smoothness in the far-field region and does not require extra efforts on the tree hierarchy.This paper is organized as follows. Sections2.1 and 2.2 present the numerical methods for Cartesian and spherical coordinates, respectively. Section 4.1 presents the order of accuracy of the proposed methods as verified by a family of finite disksD
2[9]. Section 4.2 presents a comparison with theN-body calculation. Finally, we conclude this study in Section 5.


## 2. Methods

### 2.1. Cartesian Coordinates

Assume that the numerical computation domainD
=
[
-
L
,
L
]
×
[
-
L
,
L
]for some numberL
>
0, which contains the support of the surface density, and further discretize the region uniformly as follows. Given a positive integerM, defineΔ
x
=
2
L
/
M,x
i
+
1
/
2
=
-
L
+
i
Δ
x,Δ
y
=
Δ
x, andy
j
+
1
/
2
=
-
L
+
j
Δ
y, wherei
,
j
=
0
,
…
,
M. Next, define the center of the cellD
i
,
j
=
[
x
i
-
1
/
2
,
x
i
+
1
/
2
]
×
[
y
j
-
1
/
2
,
y
j
+
1
/
2
]asx
i
=
(
x
i
-
1
/
2
+
x
i
+
1
/
2
)
/
2andy
j
=
(
y
j
-
1
/
2
+
y
j
+
1
/
2
)
/
2, wherei
,
j
=
1
,
…
,
M. Hence, the domainDis discretized into theM
2cells.We recall theN-body method with fast Fourier transform (FFT). Represent the potential functionΦof (1) as
(7)
Φ
(
x
,
y
,
z
)
=
-
∭
K
(
x
¯
-
x
,
y
¯
-
y
,
z
¯
-
z
)
ρ
(
x
¯
,
y
¯
,
z
¯
)
d
x
¯
d
y
¯
d
z
¯
,
Where  K
(
x
,
y
,
z
)
=
1
/
x
2
+
y
2
+
z
2
.According to (2), the potential is
(8)
Φ
(
x
,
y
,
z
)
=
∬
K
(
x
¯
-
x
,
y
¯
-
y
,
-
z
)
σ
(
x
¯
,
y
¯
)
d
x
¯
d
y
¯
.
Calculate (8) using a numerical approach for the planez
=
0assuming that the support of the surface density is compact. Based on a restriction to the calculation of potentials on thez
=
0plane, derive the formulae for allz. Equation (8) is equal to
(9)
Φ
(
x
,
y
,
z
)
=
∑
i
′
=
1
M
‍
∑
j
′
=
1
M
‍
∬
D
i
′
,
j
′
K
(
x
¯
-
x
,
y
¯
-
y
,
-
z
)
σ
(
x
¯
,
y
¯
)
d
x
¯
d
y
¯
.
The potential at(
x
i
,
y
j
,
z
k
)for theN-body method is approximated by
(10)
Φ
(
x
i
,
y
j
,
z
k
)
≈
∑
i
′
=
1
M
‍
∑
j
′
=
1
M
‍
K
(
x
i
′
-
x
i
,
y
j
′
-
y
j
,
-
z
k
)
σ
(
x
i
′
,
y
j
′
)
Δ
x
Δ
y
.
That is,
(11)
∬
D
i
′
,
j
′
K
(
x
¯
-
x
,
y
¯
-
y
,
-
z
)
σ
(
x
¯
,
y
¯
)
≈
K
(
x
i
′
-
x
i
,
y
j
′
-
y
j
,
-
z
k
)
σ
(
x
i
′
,
y
j
′
)
Δ
x
Δ
y
.
The approximated integral is of the first order and the fundamental function  K  is singular whenever(
i
′
,
j
′
,
0
)
=
(
i
,
j
,
k
). TheN-body method applies a softening parameter technique forz
=
0withϵ
>
0,
(12)
K
≈
K
ϵ
=
1
x
2
+
y
2
+
ϵ
2
to avoid the singularity.To improve the calculation of (11), approximate the surface densityσonD
i
,
jin (8) linearly by
(13)
σ
(
x
¯
,
y
¯
)
≈
σ
i
,
j
+
δ
i
,
j
x
(
x
¯
-
x
i
)
+
δ
i
,
j
y
(
y
¯
-
y
j
)
,
whereσ
i
,
j
=
σ
(
x
i
,
y
j
)andδ
i
,
j
x
=
(
∂
/
∂
x
)
σ
(
x
i
,
y
j
)andδ
i
,
j
y
=
(
∂
/
∂
y
)
σ
(
x
i
,
y
j
)are constant in the cellD
i
,
j. The error of the discretization isO
(
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
). If a higher order accuracy is required, additional terms in the Taylor expansion in (13) can be considered.Let(14)
K
i
-
i
′
,
j
-
j
′
,
k
0
=
∬
D
i
′
,
j
′
1
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
d
x
¯
d
y
¯
,
K
i
-
i
′
,
j
-
j
′
,
k
x
=
∬
D
i
′
,
j
′
x
¯
-
x
i
′
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
d
x
¯
d
y
¯
,
K
i
-
i
′
,
j
-
j
′
,
k
y
=
∬
D
i
′
,
j
′
y
¯
-
y
j
′
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
d
x
¯
d
y
¯
.
By substituting (13) into (9), the potential can be approximated by
(15)
Φ
i
,
j
,
k
≈
∑
i
′
=
1
M
‍
∑
j
′
=
1
M
∬
D
i
′
,
j
′
‍
K
(
x
¯
-
x
i
,
y
¯
-
y
j
,
-
z
k
)
×
(
∬
D
i
′
,
j
′
σ
i
′
,
j
′
+
δ
i
′
,
j
′
x
(
x
¯
-
x
i
′
)
+
δ
i
′
,
j
′
y
(
y
¯
-
y
j
′
)
∬
D
i
′
,
j
′
)
d
x
¯
d
y
¯
≡
Φ
i
,
j
,
k
0
+
Φ
i
,
j
,
k
x
+
Φ
i
,
j
,
k
y
,
where
(16)
Φ
i
,
j
,
k
0
=
∑
i
′
=
1
M
‍
∑
j
′
=
1
M
‍
σ
i
′
,
j
′
∬
D
i
′
,
j
′
1
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
d
x
¯
d
y
¯
=
∑
i
′
=
1
M
‍
∑
j
′
=
1
M
‍
σ
i
′
,
j
′
K
i
-
i
′
,
j
-
j
′
,
k
0
,
Φ
i
,
j
,
k
x
=
∑
i
′
=
1
M
‍
∑
j
′
=
1
M
‍
δ
i
′
,
j
′
x
∬
D
i
′
,
j
′
x
¯
-
x
i
′
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
d
x
¯
d
y
¯
=
∑
i
′
=
1
M
‍
∑
j
′
=
1
M
‍
δ
i
′
,
j
′
x
K
i
-
i
′
,
j
-
j
′
,
k
x
,
Φ
i
,
j
,
k
y
=
∑
i
′
=
1
M
‍
∑
j
′
=
1
M
‍
δ
i
′
,
j
′
y
∬
D
i
′
,
j
′
y
¯
-
y
j
′
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
d
x
¯
d
y
¯
=
∑
i
′
=
1
M
‍
∑
j
′
=
1
M
‍
δ
i
′
,
j
′
y
K
i
-
i
′
,
j
-
j
′
,
k
y
.Fortunately, the integrals of (14) can be obtained analytically with the help of the following simple formulae:(17)
∬
1
x
2
+
y
2
+
z
2
d
x
d
y
=
y
ln
⁡
(
x
+
x
2
+
y
2
+
z
2
)
+
x
ln
⁡
(
y
+
x
2
+
y
2
+
z
2
)
-
z
arctan
(
x
y
z
x
2
+
y
2
+
z
2
)
+
C
,
∬
x
x
2
+
y
2
+
z
2
d
x
d
y
=
1
2
y
x
2
+
y
2
+
z
2
+
1
2
(
x
2
+
z
2
)
ln
⁡
(
y
+
x
2
+
y
2
+
z
2
)
+
C
.
The valueK
i
-
i
′
,
j
-
j
′
,
k
0is then equal to
(18)
K
i
-
i
′
,
j
-
j
′
,
k
0
=
(
y
¯
-
y
j
)
ln
⁡
⁡
(
(
x
¯
-
x
i
)
+
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
)
+
(
x
¯
-
x
i
)
ln
⁡
⁡
(
(
y
¯
-
y
j
)
+
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
)
-
z
k
arctan
(
(
x
¯
-
x
i
)
(
y
¯
-
y
j
)
z
k
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
)
|
x
i
′
-
1
/
2
x
i
′
+
1
/
2
(
(
x
¯
-
x
i
)
(
y
¯
-
y
j
)
z
k
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
)
|
y
j
′
-
1
/
2
y
j
′
+
1
/
2
,
where the notationg
(
x
)
|
a
b
=
g
(
b
)
-
g
(
a
). The calculations ofK
i
-
i
′
,
j
-
j
′
,
k
xandK
i
-
i
′
,
j
-
j
′
,
k
yare split into two parts by the identityx
¯
-
x
i
′
=
(
x
¯
-
x
i
)
+
(
x
i
-
x
i
′
) andy
¯
-
y
j
′
=
(
y
¯
-
y
j
)
+
(
y
j
-
y
j
′
), respectively. This leads to
(19)
K
i
-
i
′
,
j
-
j
′
,
k
x
=
(
x
i
-
x
i
′
)
K
i
-
i
′
,
j
-
j
′
,
k
0
+
1
2
(
y
¯
-
y
j
)
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
+
1
2
(
(
x
¯
-
x
i
)
2
+
z
k
2
)
×
ln
⁡
(
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
(
y
¯
-
y
j
)
+
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
)
|
x
i
′
-
1
/
2
x
i
′
+
1
/
2
+
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
)
|
y
j
′
-
1
/
2
y
j
′
+
1
/
2
,
K
i
-
i
′
,
j
-
j
′
,
k
y
=
(
y
j
-
y
j
′
)
K
i
-
i
′
,
j
-
j
′
,
k
0
+
1
2
(
x
¯
-
x
i
)
(
y
¯
-
y
j
)
2
+
(
x
¯
-
x
i
)
2
+
z
k
2
+
1
2
(
(
y
¯
-
y
j
)
2
+
z
k
2
)
×
ln
⁡
(
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
(
x
¯
-
x
i
)
+
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
)
|
x
i
′
-
1
/
2
x
i
′
+
1
/
2
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
|
y
j
′
-
1
/
2
y
j
′
+
1
/
2
.The forms ofΦ
i
,
j
,
k
0,Φ
i
,
j
,
k
x, andΦ
i
,
j
,
k
yin (16) are a type of convolution and the computational complexity can be reduced toO
(
M
log
⁡
2
M
)with the help of FFT in thexandydirections. In thezdirection, the computational complexity isO
(
M
). This implies that the total complexity of this method isO
(
M
2
(
log
⁡
2
M
)
2
)andO
(
M
3
(
log
⁡
2
M
)
2
)for the planez
=
0and the whole 3D space, respectively.


### 2.2. Spherical and Polar Coordinates

While the method appears to be second-order-accurate for discretization on a Cartesian grid, it loses second-order accuracy using spherical polar coordinates. An algorithm for solving the Poisson equation in spherical polar coordinates has been used extensively in protostellar collapse and fragmentation calculations, where often a thin disk formed at the central core of the collapsing cloud [14]. The algorithm implemented in [15] is subsequently modified, see [16], to deal with filamentary density distributions. This method in [16] employed finite differences for the radial dependence of the potential and spherical harmonics for the angular dependence, resulting in a second-order method through convergence testing.We are going to develop our approach for spherical coordinates(
r
,
θ
,
φ
), which can be reduced to polar coordinates (φ
=
π
/
2). The singular integral disappears, but the high order of accuracy is not attained because there is no explicit form for the integral of an elliptic function. Although the proposed calculation of potential has some drawbacks, it is still better thanN-body calculations and the method in [17].The support of the surface density is compact, contained in a regionR
=
[
0
,
L
]
×
[
0,2
π
]for some numberL
>
0. Confine the discretization on the radial region in logarithmic form and the azimuthal region uniformly as follows to achieve a fast calculation. Given a positive integerM, defineΔ
θ
=
2
π
/
M,0
<
β
0
<
1,β
=
β
0
(
1
-
Δ
θ
),r
i
+
1
/
2
=
β
M
-
i
L,θ
j
+
1
/
2
=
j
Δ
θ,i
,
j
=
0
,
…
,
M,r
i
=
(
1
/
2
)
(
r
i
-
1
/
2
+
r
i
+
1
/
2
), andθ
j
=
(
1
/
2
)
(
θ
j
-
1
/
2
+
θ
j
+
1
/
2
)wherei
,
j
=
1
,
…
,
M. The effect ofβ
0is to refine the mesh. Here, the proposed method for polar coordinates is of first order because a singular integration occurs (see below). Furthermore, the region  R  is discretized into theM
2cells,R
i
j
=
[
r
i
-
1
/
2
,
r
i
+
1
/
2
]
×
[
θ
j
-
1
/
2
,
θ
j
+
1
/
2
]fori
,
j
=
1
,
…
,
M. Forj
=
1
,
…
,
M, the cellsR
1
,
jdo not cover the origin and extra cellsR
^
j
=
[
0
,
r
1
/
2
]
×
[
θ
j
-
1
/
2
,
θ
j
+
1
/
2
]should be included. To simplify notation, denoteR
0
,
j
=
R
^
j,j
=
1
,
…
,
M.The potential functionΦof (1) under the assumptionG
=
1in spherical coordinates can be expressed as
(20)
Φ
(
r
,
θ
,
φ
)
=
-
∭
K
(
r
¯
,
r
,
θ
¯
-
θ
,
φ
¯
,
φ
)
ρ
(
r
¯
,
θ
¯
,
φ
¯
)
r
¯
2
sin
φ
¯
d
r
¯
d
θ
¯
d
φ
¯
,
where
(21)
K
(
r
¯
,
r
,
θ
¯
-
θ
,
φ
¯
,
φ
)
=
1
×
(
r
¯
2
-
2
r
¯
r
cos
⁡
(
θ
¯
-
θ
)
sin
⁡
(
φ
¯
)
sin
⁡
(
φ
)
+
r
2
-
2
r
¯
r
cos
⁡
(
φ
¯
)
cos
⁡
⁡
(
φ
)
(
θ
¯
-
θ
)
)
-
1
/
2
.The surface densityσonR
i
,
jin (20) can be linearly approximated by
(22)
σ
(
r
¯
,
θ
¯
)
≈
σ
i
,
j
+
δ
i
,
j
r
(
r
¯
-
r
i
)
+
δ
i
,
j
θ
(
θ
¯
-
θ
j
)
,
whereσ
i
,
j
=
σ
(
r
i
,
θ
j
)andδ
i
,
j
r
=
(
∂
/
∂
r
)
σ
(
r
i
,
θ
j
)andδ
i
,
j
θ
=
(
∂
/
∂
θ
)
σ
(
r
i
,
θ
j
)are constants in the cellR
i
,
j. Equation (22) is the Taylor expansion in two dimensions and it follows that the error of the discretization isO
(
(
r
¯
-
r
i
)
2
+
(
θ
¯
-
θ
j
)
2
).In this section, let(23)
K
i
-
i
′
,
j
-
j
′
,
k
0
=
∬
R
i
′
,
j
′
r
¯
r
¯
2
+
r
i
2
-
2
r
¯
r
i
cos
⁡
(
θ
¯
-
θ
j
)
sin
(
φ
k
)
d
r
¯
d
θ
¯
,
(24)
r
i
K
i
-
i
′
,
j
-
j
′
,
k
r
=
∬
R
i
′
,
j
′
r
¯
(
r
¯
-
r
i
′
)
r
¯
2
+
r
i
2
-
2
r
¯
r
i
cos
⁡
(
θ
¯
-
θ
j
)
sin
(
φ
k
)
d
r
¯
d
θ
¯
,
(25)
K
i
-
i
′
,
j
-
j
′
,
k
θ
=
∬
R
i
′
,
j
′
r
¯
(
θ
¯
-
θ
j
′
)
r
¯
2
+
r
i
2
-
2
r
¯
r
i
cos
⁡
(
θ
¯
-
θ
j
)
sin
(
φ
k
)
d
r
¯
d
θ
¯
.
Similarly forK
¯
i
,
j
-
j
′
,
k
0,K
¯
i
,
j
-
j
′
,
k
r, andK
¯
i
,
j
-
j
′
,
k
θon the central region R
0
,
j. The termr
iin (24) is for the formulation of a convolution type. According to (20) and (22),
(26)
Φ
i
,
j
,
k
≈
∑
i
′
=
0
M
‍
∑
j
′
=
1
M
‍
∬
R
i
′
,
j
′
K
(
r
¯
,
r
i
,
θ
¯
-
θ
j
,
φ
k
,
π
2
)
×
(
σ
i
′
,
j
′
+
δ
i
′
,
j
′
r
(
r
¯
-
r
i
′
)
+
δ
i
′
,
j
′
θ
(
θ
¯
-
θ
j
′
)
)
r
¯
d
r
¯
d
θ
¯
∶
=
Φ
i
,
j
,
k
0
,
p
+
Φ
i
,
j
,
k
r
+
Φ
i
,
j
,
k
θ
,
where
(27)
Φ
i
,
j
,
k
0
,
p
=
∑
i
′
=
0
M
‍
∑
j
′
=
1
M
‍
σ
i
′
,
j
′
×
∬
R
i
′
,
j
′
r
¯
r
¯
2
+
r
i
2
-
2
r
¯
r
i
cos
⁡
⁡
(
θ
¯
-
θ
j
)
sin
⁡
(
φ
k
)
d
r
¯
d
θ
¯
,
Φ
i
,
j
,
k
r
=
∑
i
′
=
0
M
‍
∑
j
′
=
1
M
‍
δ
i
′
,
j
′
r
×
∬
R
i
′
,
j
′
r
¯
(
r
¯
-
r
i
′
)
r
¯
2
+
r
i
2
-
2
r
¯
r
i
cos
⁡
⁡
(
θ
¯
-
θ
j
)
sin
⁡
(
φ
k
)
d
r
¯
d
θ
¯
,
Φ
i
,
j
,
k
θ
=
∑
i
′
=
0
M
‍
∑
j
′
=
1
M
‍
δ
i
′
,
j
′
θ
×
∬
R
i
′
,
j
′
r
¯
(
θ
¯
-
θ
j
′
)
r
¯
2
+
r
i
2
-
2
r
¯
r
i
cos
⁡
⁡
(
θ
¯
-
θ
j
)
sin
⁡
(
φ
k
)
d
r
¯
d
θ
¯
.
Equation (27) can then be rewritten as
(28)
Φ
i
,
j
,
k
0
,
p
=
∑
i
′
=
1
M
‍
∑
j
′
=
1
M
‍
σ
i
′
,
j
′
K
i
-
i
′
,
j
-
j
′
,
k
0
+
∑
j
′
=
1
M
‍
σ
0
,
j
′
K
¯
i
,
j
-
j
′
,
k
0
,
Φ
i
,
j
,
k
r
=
r
i
∑
i
′
=
1
M
‍
∑
j
′
=
1
M
‍
δ
i
′
,
j
′
r
K
i
-
i
,
j
-
j
′
,
k
r
+
r
i
∑
j
′
=
1
M
‍
δ
0
,
j
′
r
K
¯
i
,
j
-
j
′
,
k
r
,
Φ
i
,
j
,
k
θ
=
∑
i
′
=
1
M
‍
∑
j
′
=
1
M
‍
δ
i
′
,
j
′
θ
K
i
-
i
′
,
j
-
j
′
,
k
θ
+
∑
j
′
=
1
M
‍
δ
0
,
j
′
θ
K
¯
i
,
j
-
j
′
,
k
θ
.DefineF
(
r
~
,
θ
,
φ
)
=
1
+
r
~
2
-
2
r
~
cos
⁡
(
θ
)
sin
(
φ
), wherer
~is the dimensionless radius. The evaluation of (23), (24), and (25) can be obtained with the help of the following simple integrals:
(29)
∫
‍
r
¯
r
¯
2
+
r
2
-
2
r
r
¯
cos
⁡
(
θ
)
sin
(
φ
)
d
r
¯
=
r
F
(
r
¯
r
,
θ
,
φ
)
+
r
cos
⁡
(
θ
)
sin
⁡
(
φ
)
ln
⁡
(
-
cos
⁡
⁡
(
θ
)
sin
⁡
(
φ
)
+
r
¯
r
+
F
(
r
¯
r
,
θ
,
φ
)
)
+
C
∶
=
r
H
1
(
r
¯
r
,
θ
,
φ
)
+
C
,
∫
‍
r
¯
2
r
¯
2
+
r
2
-
2
r
r
¯
cos
⁡
(
θ
)
sin
(
φ
)
d
r
¯
=
r
2
(
1
2
r
¯
r
+
3
2
cos
⁡
(
θ
)
sin
(
φ
)
)
F
(
r
¯
r
,
θ
,
φ
)
+
r
2
(
3
2
cos
⁡
2
(
θ
)
sin
2
(
φ
)
-
1
2
)
×
ln
⁡
(
-
cos
⁡
⁡
(
θ
)
sin
⁡
(
φ
)
+
r
¯
r
+
F
(
r
¯
r
,
θ
,
φ
)
)
+
C
∶
=
r
2
H
2
(
r
¯
r
,
θ
,
φ
)
+
C
.The definitions ofr
i
′
+
1
/
2andr
ilead to
(30)
r
i
′
+
1
/
2
r
i
=
2
β
i
-
i
′
1
+
β
,
r
i
′
r
i
=
β
i
-
i
′
.
Calculate the value of the integral
(31)
K
i
-
i
′
,
j
-
j
′
,
k
0
=
∫
θ
j
′
-
1
/
2
θ
j
′
+
1
/
2
‍
∫
r
i
′
-
1
/
2
r
i
′
+
1
/
2
‍
r
¯
r
¯
2
+
r
i
2
-
2
r
¯
r
i
cos
⁡
(
θ
¯
-
θ
j
)
sin
(
φ
k
)
d
r
¯
d
θ
¯
=
∫
θ
j
′
-
1
/
2
θ
j
′
+
1
/
2
‍
r
i
F
(
r
¯
r
i
,
θ
¯
-
θ
j
,
φ
k
)
+
r
i
cos
⁡
⁡
(
θ
¯
-
θ
j
)
sin
⁡
(
φ
k
)
×
ln
⁡
(
(
r
¯
r
i
,
θ
¯
-
θ
j
,
φ
k
)
-
cos
⁡
(
θ
)
sin
(
φ
k
)
+
r
¯
r
i
+
F
(
r
¯
r
i
,
θ
¯
-
θ
j
,
φ
k
)
)
|
r
i
′
-
1
/
2
r
i
′
+
1
/
2
d
θ
¯
.
The last integral in this equation is an elliptic integral, and this study uses a trapezoidal rule for its evaluation. This numerical approach approximates the valueK
i
-
i
′
,
j
-
j
′
,
k
0as follows:
(32)
K
i
-
i
′
,
j
-
j
′
,
k
0
≈
1
2
(
r
i
H
1
(
r
i
′
+
1
/
2
r
i
,
θ
j
′
+
1
/
2
-
θ
j
,
φ
k
)
-
r
i
H
1
(
r
i
′
-
1
/
2
r
i
,
θ
j
′
-
1
/
2
-
θ
j
,
φ
k
)
+
r
i
H
1
(
r
i
′
+
1
/
2
r
i
,
θ
j
′
-
1
/
2
-
θ
j
,
φ
k
)
-
r
i
H
1
(
r
i
′
-
1
/
2
r
i
,
θ
j
′
-
1
/
2
-
θ
j
,
φ
k
)
)
×
(
θ
j
′
+
1
/
2
-
θ
j
′
-
1
/
2
)
∶
=
r
i
H
1
(
r
¯
r
i
,
θ
¯
-
θ
j
,
φ
k
)
|
r
i
′
-
1
/
2
r
i
′
+
1
/
2
]
θ
j
′
-
1
/
2
θ
j
′
+
1
/
2
,
where the notationf
(
·
)
]
a
b
=
(
1
/
2
)
(
f
(
b
)
+
f
(
a
)
)
(
b
-
a
). Similarly,
(33)
K
i
-
i
′
,
j
-
j
′
,
k
r
≈
r
i
2
(
H
2
(
r
¯
r
i
,
θ
¯
-
θ
j
,
φ
k
)
|
r
i
′
-
1
/
2
r
i
′
+
1
/
2
]
θ
j
′
-
1
/
2
θ
j
′
+
1
/
2
-
r
i
′
r
i
K
i
-
i
′
,
j
-
j
′
,
k
0
)
,
K
i
-
i
′
,
j
-
j
′
,
k
θ
≈
(
θ
¯
-
θ
j
)
r
i
H
1
(
r
¯
r
i
,
θ
¯
-
θ
j
,
φ
k
)
|
r
i
′
-
1
/
2
r
i
′
+
1
/
2
]
θ
j
′
-
1
/
2
θ
j
′
+
1
/
2
.


## 2.1. Cartesian Coordinates

Assume that the numerical computation domainD
=
[
-
L
,
L
]
×
[
-
L
,
L
]for some numberL
>
0, which contains the support of the surface density, and further discretize the region uniformly as follows. Given a positive integerM, defineΔ
x
=
2
L
/
M,x
i
+
1
/
2
=
-
L
+
i
Δ
x,Δ
y
=
Δ
x, andy
j
+
1
/
2
=
-
L
+
j
Δ
y, wherei
,
j
=
0
,
…
,
M. Next, define the center of the cellD
i
,
j
=
[
x
i
-
1
/
2
,
x
i
+
1
/
2
]
×
[
y
j
-
1
/
2
,
y
j
+
1
/
2
]asx
i
=
(
x
i
-
1
/
2
+
x
i
+
1
/
2
)
/
2andy
j
=
(
y
j
-
1
/
2
+
y
j
+
1
/
2
)
/
2, wherei
,
j
=
1
,
…
,
M. Hence, the domainDis discretized into theM
2cells.We recall theN-body method with fast Fourier transform (FFT). Represent the potential functionΦof (1) as
(7)
Φ
(
x
,
y
,
z
)
=
-
∭
K
(
x
¯
-
x
,
y
¯
-
y
,
z
¯
-
z
)
ρ
(
x
¯
,
y
¯
,
z
¯
)
d
x
¯
d
y
¯
d
z
¯
,
Where  K
(
x
,
y
,
z
)
=
1
/
x
2
+
y
2
+
z
2
.According to (2), the potential is
(8)
Φ
(
x
,
y
,
z
)
=
∬
K
(
x
¯
-
x
,
y
¯
-
y
,
-
z
)
σ
(
x
¯
,
y
¯
)
d
x
¯
d
y
¯
.
Calculate (8) using a numerical approach for the planez
=
0assuming that the support of the surface density is compact. Based on a restriction to the calculation of potentials on thez
=
0plane, derive the formulae for allz. Equation (8) is equal to
(9)
Φ
(
x
,
y
,
z
)
=
∑
i
′
=
1
M
‍
∑
j
′
=
1
M
‍
∬
D
i
′
,
j
′
K
(
x
¯
-
x
,
y
¯
-
y
,
-
z
)
σ
(
x
¯
,
y
¯
)
d
x
¯
d
y
¯
.
The potential at(
x
i
,
y
j
,
z
k
)for theN-body method is approximated by
(10)
Φ
(
x
i
,
y
j
,
z
k
)
≈
∑
i
′
=
1
M
‍
∑
j
′
=
1
M
‍
K
(
x
i
′
-
x
i
,
y
j
′
-
y
j
,
-
z
k
)
σ
(
x
i
′
,
y
j
′
)
Δ
x
Δ
y
.
That is,
(11)
∬
D
i
′
,
j
′
K
(
x
¯
-
x
,
y
¯
-
y
,
-
z
)
σ
(
x
¯
,
y
¯
)
≈
K
(
x
i
′
-
x
i
,
y
j
′
-
y
j
,
-
z
k
)
σ
(
x
i
′
,
y
j
′
)
Δ
x
Δ
y
.
The approximated integral is of the first order and the fundamental function  K  is singular whenever(
i
′
,
j
′
,
0
)
=
(
i
,
j
,
k
). TheN-body method applies a softening parameter technique forz
=
0withϵ
>
0,
(12)
K
≈
K
ϵ
=
1
x
2
+
y
2
+
ϵ
2
to avoid the singularity.To improve the calculation of (11), approximate the surface densityσonD
i
,
jin (8) linearly by
(13)
σ
(
x
¯
,
y
¯
)
≈
σ
i
,
j
+
δ
i
,
j
x
(
x
¯
-
x
i
)
+
δ
i
,
j
y
(
y
¯
-
y
j
)
,
whereσ
i
,
j
=
σ
(
x
i
,
y
j
)andδ
i
,
j
x
=
(
∂
/
∂
x
)
σ
(
x
i
,
y
j
)andδ
i
,
j
y
=
(
∂
/
∂
y
)
σ
(
x
i
,
y
j
)are constant in the cellD
i
,
j. The error of the discretization isO
(
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
). If a higher order accuracy is required, additional terms in the Taylor expansion in (13) can be considered.Let(14)
K
i
-
i
′
,
j
-
j
′
,
k
0
=
∬
D
i
′
,
j
′
1
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
d
x
¯
d
y
¯
,
K
i
-
i
′
,
j
-
j
′
,
k
x
=
∬
D
i
′
,
j
′
x
¯
-
x
i
′
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
d
x
¯
d
y
¯
,
K
i
-
i
′
,
j
-
j
′
,
k
y
=
∬
D
i
′
,
j
′
y
¯
-
y
j
′
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
d
x
¯
d
y
¯
.
By substituting (13) into (9), the potential can be approximated by
(15)
Φ
i
,
j
,
k
≈
∑
i
′
=
1
M
‍
∑
j
′
=
1
M
∬
D
i
′
,
j
′
‍
K
(
x
¯
-
x
i
,
y
¯
-
y
j
,
-
z
k
)
×
(
∬
D
i
′
,
j
′
σ
i
′
,
j
′
+
δ
i
′
,
j
′
x
(
x
¯
-
x
i
′
)
+
δ
i
′
,
j
′
y
(
y
¯
-
y
j
′
)
∬
D
i
′
,
j
′
)
d
x
¯
d
y
¯
≡
Φ
i
,
j
,
k
0
+
Φ
i
,
j
,
k
x
+
Φ
i
,
j
,
k
y
,
where
(16)
Φ
i
,
j
,
k
0
=
∑
i
′
=
1
M
‍
∑
j
′
=
1
M
‍
σ
i
′
,
j
′
∬
D
i
′
,
j
′
1
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
d
x
¯
d
y
¯
=
∑
i
′
=
1
M
‍
∑
j
′
=
1
M
‍
σ
i
′
,
j
′
K
i
-
i
′
,
j
-
j
′
,
k
0
,
Φ
i
,
j
,
k
x
=
∑
i
′
=
1
M
‍
∑
j
′
=
1
M
‍
δ
i
′
,
j
′
x
∬
D
i
′
,
j
′
x
¯
-
x
i
′
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
d
x
¯
d
y
¯
=
∑
i
′
=
1
M
‍
∑
j
′
=
1
M
‍
δ
i
′
,
j
′
x
K
i
-
i
′
,
j
-
j
′
,
k
x
,
Φ
i
,
j
,
k
y
=
∑
i
′
=
1
M
‍
∑
j
′
=
1
M
‍
δ
i
′
,
j
′
y
∬
D
i
′
,
j
′
y
¯
-
y
j
′
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
d
x
¯
d
y
¯
=
∑
i
′
=
1
M
‍
∑
j
′
=
1
M
‍
δ
i
′
,
j
′
y
K
i
-
i
′
,
j
-
j
′
,
k
y
.Fortunately, the integrals of (14) can be obtained analytically with the help of the following simple formulae:(17)
∬
1
x
2
+
y
2
+
z
2
d
x
d
y
=
y
ln
⁡
(
x
+
x
2
+
y
2
+
z
2
)
+
x
ln
⁡
(
y
+
x
2
+
y
2
+
z
2
)
-
z
arctan
(
x
y
z
x
2
+
y
2
+
z
2
)
+
C
,
∬
x
x
2
+
y
2
+
z
2
d
x
d
y
=
1
2
y
x
2
+
y
2
+
z
2
+
1
2
(
x
2
+
z
2
)
ln
⁡
(
y
+
x
2
+
y
2
+
z
2
)
+
C
.
The valueK
i
-
i
′
,
j
-
j
′
,
k
0is then equal to
(18)
K
i
-
i
′
,
j
-
j
′
,
k
0
=
(
y
¯
-
y
j
)
ln
⁡
⁡
(
(
x
¯
-
x
i
)
+
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
)
+
(
x
¯
-
x
i
)
ln
⁡
⁡
(
(
y
¯
-
y
j
)
+
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
)
-
z
k
arctan
(
(
x
¯
-
x
i
)
(
y
¯
-
y
j
)
z
k
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
)
|
x
i
′
-
1
/
2
x
i
′
+
1
/
2
(
(
x
¯
-
x
i
)
(
y
¯
-
y
j
)
z
k
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
)
|
y
j
′
-
1
/
2
y
j
′
+
1
/
2
,
where the notationg
(
x
)
|
a
b
=
g
(
b
)
-
g
(
a
). The calculations ofK
i
-
i
′
,
j
-
j
′
,
k
xandK
i
-
i
′
,
j
-
j
′
,
k
yare split into two parts by the identityx
¯
-
x
i
′
=
(
x
¯
-
x
i
)
+
(
x
i
-
x
i
′
) andy
¯
-
y
j
′
=
(
y
¯
-
y
j
)
+
(
y
j
-
y
j
′
), respectively. This leads to
(19)
K
i
-
i
′
,
j
-
j
′
,
k
x
=
(
x
i
-
x
i
′
)
K
i
-
i
′
,
j
-
j
′
,
k
0
+
1
2
(
y
¯
-
y
j
)
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
+
1
2
(
(
x
¯
-
x
i
)
2
+
z
k
2
)
×
ln
⁡
(
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
(
y
¯
-
y
j
)
+
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
)
|
x
i
′
-
1
/
2
x
i
′
+
1
/
2
+
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
)
|
y
j
′
-
1
/
2
y
j
′
+
1
/
2
,
K
i
-
i
′
,
j
-
j
′
,
k
y
=
(
y
j
-
y
j
′
)
K
i
-
i
′
,
j
-
j
′
,
k
0
+
1
2
(
x
¯
-
x
i
)
(
y
¯
-
y
j
)
2
+
(
x
¯
-
x
i
)
2
+
z
k
2
+
1
2
(
(
y
¯
-
y
j
)
2
+
z
k
2
)
×
ln
⁡
(
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
(
x
¯
-
x
i
)
+
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
)
|
x
i
′
-
1
/
2
x
i
′
+
1
/
2
(
x
¯
-
x
i
)
2
+
(
y
¯
-
y
j
)
2
+
z
k
2
|
y
j
′
-
1
/
2
y
j
′
+
1
/
2
.The forms ofΦ
i
,
j
,
k
0,Φ
i
,
j
,
k
x, andΦ
i
,
j
,
k
yin (16) are a type of convolution and the computational complexity can be reduced toO
(
M
log
⁡
2
M
)with the help of FFT in thexandydirections. In thezdirection, the computational complexity isO
(
M
). This implies that the total complexity of this method isO
(
M
2
(
log
⁡
2
M
)
2
)andO
(
M
3
(
log
⁡
2
M
)
2
)for the planez
=
0and the whole 3D space, respectively.


## 2.2. Spherical and Polar Coordinates

While the method appears to be second-order-accurate for discretization on a Cartesian grid, it loses second-order accuracy using spherical polar coordinates. An algorithm for solving the Poisson equation in spherical polar coordinates has been used extensively in protostellar collapse and fragmentation calculations, where often a thin disk formed at the central core of the collapsing cloud [14]. The algorithm implemented in [15] is subsequently modified, see [16], to deal with filamentary density distributions. This method in [16] employed finite differences for the radial dependence of the potential and spherical harmonics for the angular dependence, resulting in a second-order method through convergence testing.We are going to develop our approach for spherical coordinates(
r
,
θ
,
φ
), which can be reduced to polar coordinates (φ
=
π
/
2). The singular integral disappears, but the high order of accuracy is not attained because there is no explicit form for the integral of an elliptic function. Although the proposed calculation of potential has some drawbacks, it is still better thanN-body calculations and the method in [17].The support of the surface density is compact, contained in a regionR
=
[
0
,
L
]
×
[
0,2
π
]for some numberL
>
0. Confine the discretization on the radial region in logarithmic form and the azimuthal region uniformly as follows to achieve a fast calculation. Given a positive integerM, defineΔ
θ
=
2
π
/
M,0
<
β
0
<
1,β
=
β
0
(
1
-
Δ
θ
),r
i
+
1
/
2
=
β
M
-
i
L,θ
j
+
1
/
2
=
j
Δ
θ,i
,
j
=
0
,
…
,
M,r
i
=
(
1
/
2
)
(
r
i
-
1
/
2
+
r
i
+
1
/
2
), andθ
j
=
(
1
/
2
)
(
θ
j
-
1
/
2
+
θ
j
+
1
/
2
)wherei
,
j
=
1
,
…
,
M. The effect ofβ
0is to refine the mesh. Here, the proposed method for polar coordinates is of first order because a singular integration occurs (see below). Furthermore, the region  R  is discretized into theM
2cells,R
i
j
=
[
r
i
-
1
/
2
,
r
i
+
1
/
2
]
×
[
θ
j
-
1
/
2
,
θ
j
+
1
/
2
]fori
,
j
=
1
,
…
,
M. Forj
=
1
,
…
,
M, the cellsR
1
,
jdo not cover the origin and extra cellsR
^
j
=
[
0
,
r
1
/
2
]
×
[
θ
j
-
1
/
2
,
θ
j
+
1
/
2
]should be included. To simplify notation, denoteR
0
,
j
=
R
^
j,j
=
1
,
…
,
M.The potential functionΦof (1) under the assumptionG
=
1in spherical coordinates can be expressed as
(20)
Φ
(
r
,
θ
,
φ
)
=
-
∭
K
(
r
¯
,
r
,
θ
¯
-
θ
,
φ
¯
,
φ
)
ρ
(
r
¯
,
θ
¯
,
φ
¯
)
r
¯
2
sin
φ
¯
d
r
¯
d
θ
¯
d
φ
¯
,
where
(21)
K
(
r
¯
,
r
,
θ
¯
-
θ
,
φ
¯
,
φ
)
=
1
×
(
r
¯
2
-
2
r
¯
r
cos
⁡
(
θ
¯
-
θ
)
sin
⁡
(
φ
¯
)
sin
⁡
(
φ
)
+
r
2
-
2
r
¯
r
cos
⁡
(
φ
¯
)
cos
⁡
⁡
(
φ
)
(
θ
¯
-
θ
)
)
-
1
/
2
.The surface densityσonR
i
,
jin (20) can be linearly approximated by
(22)
σ
(
r
¯
,
θ
¯
)
≈
σ
i
,
j
+
δ
i
,
j
r
(
r
¯
-
r
i
)
+
δ
i
,
j
θ
(
θ
¯
-
θ
j
)
,
whereσ
i
,
j
=
σ
(
r
i
,
θ
j
)andδ
i
,
j
r
=
(
∂
/
∂
r
)
σ
(
r
i
,
θ
j
)andδ
i
,
j
θ
=
(
∂
/
∂
θ
)
σ
(
r
i
,
θ
j
)are constants in the cellR
i
,
j. Equation (22) is the Taylor expansion in two dimensions and it follows that the error of the discretization isO
(
(
r
¯
-
r
i
)
2
+
(
θ
¯
-
θ
j
)
2
).In this section, let(23)
K
i
-
i
′
,
j
-
j
′
,
k
0
=
∬
R
i
′
,
j
′
r
¯
r
¯
2
+
r
i
2
-
2
r
¯
r
i
cos
⁡
(
θ
¯
-
θ
j
)
sin
(
φ
k
)
d
r
¯
d
θ
¯
,
(24)
r
i
K
i
-
i
′
,
j
-
j
′
,
k
r
=
∬
R
i
′
,
j
′
r
¯
(
r
¯
-
r
i
′
)
r
¯
2
+
r
i
2
-
2
r
¯
r
i
cos
⁡
(
θ
¯
-
θ
j
)
sin
(
φ
k
)
d
r
¯
d
θ
¯
,
(25)
K
i
-
i
′
,
j
-
j
′
,
k
θ
=
∬
R
i
′
,
j
′
r
¯
(
θ
¯
-
θ
j
′
)
r
¯
2
+
r
i
2
-
2
r
¯
r
i
cos
⁡
(
θ
¯
-
θ
j
)
sin
(
φ
k
)
d
r
¯
d
θ
¯
.
Similarly forK
¯
i
,
j
-
j
′
,
k
0,K
¯
i
,
j
-
j
′
,
k
r, andK
¯
i
,
j
-
j
′
,
k
θon the central region R
0
,
j. The termr
iin (24) is for the formulation of a convolution type. According to (20) and (22),
(26)
Φ
i
,
j
,
k
≈
∑
i
′
=
0
M
‍
∑
j
′
=
1
M
‍
∬
R
i
′
,
j
′
K
(
r
¯
,
r
i
,
θ
¯
-
θ
j
,
φ
k
,
π
2
)
×
(
σ
i
′
,
j
′
+
δ
i
′
,
j
′
r
(
r
¯
-
r
i
′
)
+
δ
i
′
,
j
′
θ
(
θ
¯
-
θ
j
′
)
)
r
¯
d
r
¯
d
θ
¯
∶
=
Φ
i
,
j
,
k
0
,
p
+
Φ
i
,
j
,
k
r
+
Φ
i
,
j
,
k
θ
,
where
(27)
Φ
i
,
j
,
k
0
,
p
=
∑
i
′
=
0
M
‍
∑
j
′
=
1
M
‍
σ
i
′
,
j
′
×
∬
R
i
′
,
j
′
r
¯
r
¯
2
+
r
i
2
-
2
r
¯
r
i
cos
⁡
⁡
(
θ
¯
-
θ
j
)
sin
⁡
(
φ
k
)
d
r
¯
d
θ
¯
,
Φ
i
,
j
,
k
r
=
∑
i
′
=
0
M
‍
∑
j
′
=
1
M
‍
δ
i
′
,
j
′
r
×
∬
R
i
′
,
j
′
r
¯
(
r
¯
-
r
i
′
)
r
¯
2
+
r
i
2
-
2
r
¯
r
i
cos
⁡
⁡
(
θ
¯
-
θ
j
)
sin
⁡
(
φ
k
)
d
r
¯
d
θ
¯
,
Φ
i
,
j
,
k
θ
=
∑
i
′
=
0
M
‍
∑
j
′
=
1
M
‍
δ
i
′
,
j
′
θ
×
∬
R
i
′
,
j
′
r
¯
(
θ
¯
-
θ
j
′
)
r
¯
2
+
r
i
2
-
2
r
¯
r
i
cos
⁡
⁡
(
θ
¯
-
θ
j
)
sin
⁡
(
φ
k
)
d
r
¯
d
θ
¯
.
Equation (27) can then be rewritten as
(28)
Φ
i
,
j
,
k
0
,
p
=
∑
i
′
=
1
M
‍
∑
j
′
=
1
M
‍
σ
i
′
,
j
′
K
i
-
i
′
,
j
-
j
′
,
k
0
+
∑
j
′
=
1
M
‍
σ
0
,
j
′
K
¯
i
,
j
-
j
′
,
k
0
,
Φ
i
,
j
,
k
r
=
r
i
∑
i
′
=
1
M
‍
∑
j
′
=
1
M
‍
δ
i
′
,
j
′
r
K
i
-
i
,
j
-
j
′
,
k
r
+
r
i
∑
j
′
=
1
M
‍
δ
0
,
j
′
r
K
¯
i
,
j
-
j
′
,
k
r
,
Φ
i
,
j
,
k
θ
=
∑
i
′
=
1
M
‍
∑
j
′
=
1
M
‍
δ
i
′
,
j
′
θ
K
i
-
i
′
,
j
-
j
′
,
k
θ
+
∑
j
′
=
1
M
‍
δ
0
,
j
′
θ
K
¯
i
,
j
-
j
′
,
k
θ
.DefineF
(
r
~
,
θ
,
φ
)
=
1
+
r
~
2
-
2
r
~
cos
⁡
(
θ
)
sin
(
φ
), wherer
~is the dimensionless radius. The evaluation of (23), (24), and (25) can be obtained with the help of the following simple integrals:
(29)
∫
‍
r
¯
r
¯
2
+
r
2
-
2
r
r
¯
cos
⁡
(
θ
)
sin
(
φ
)
d
r
¯
=
r
F
(
r
¯
r
,
θ
,
φ
)
+
r
cos
⁡
(
θ
)
sin
⁡
(
φ
)
ln
⁡
(
-
cos
⁡
⁡
(
θ
)
sin
⁡
(
φ
)
+
r
¯
r
+
F
(
r
¯
r
,
θ
,
φ
)
)
+
C
∶
=
r
H
1
(
r
¯
r
,
θ
,
φ
)
+
C
,
∫
‍
r
¯
2
r
¯
2
+
r
2
-
2
r
r
¯
cos
⁡
(
θ
)
sin
(
φ
)
d
r
¯
=
r
2
(
1
2
r
¯
r
+
3
2
cos
⁡
(
θ
)
sin
(
φ
)
)
F
(
r
¯
r
,
θ
,
φ
)
+
r
2
(
3
2
cos
⁡
2
(
θ
)
sin
2
(
φ
)
-
1
2
)
×
ln
⁡
(
-
cos
⁡
⁡
(
θ
)
sin
⁡
(
φ
)
+
r
¯
r
+
F
(
r
¯
r
,
θ
,
φ
)
)
+
C
∶
=
r
2
H
2
(
r
¯
r
,
θ
,
φ
)
+
C
.The definitions ofr
i
′
+
1
/
2andr
ilead to
(30)
r
i
′
+
1
/
2
r
i
=
2
β
i
-
i
′
1
+
β
,
r
i
′
r
i
=
β
i
-
i
′
.
Calculate the value of the integral
(31)
K
i
-
i
′
,
j
-
j
′
,
k
0
=
∫
θ
j
′
-
1
/
2
θ
j
′
+
1
/
2
‍
∫
r
i
′
-
1
/
2
r
i
′
+
1
/
2
‍
r
¯
r
¯
2
+
r
i
2
-
2
r
¯
r
i
cos
⁡
(
θ
¯
-
θ
j
)
sin
(
φ
k
)
d
r
¯
d
θ
¯
=
∫
θ
j
′
-
1
/
2
θ
j
′
+
1
/
2
‍
r
i
F
(
r
¯
r
i
,
θ
¯
-
θ
j
,
φ
k
)
+
r
i
cos
⁡
⁡
(
θ
¯
-
θ
j
)
sin
⁡
(
φ
k
)
×
ln
⁡
(
(
r
¯
r
i
,
θ
¯
-
θ
j
,
φ
k
)
-
cos
⁡
(
θ
)
sin
(
φ
k
)
+
r
¯
r
i
+
F
(
r
¯
r
i
,
θ
¯
-
θ
j
,
φ
k
)
)
|
r
i
′
-
1
/
2
r
i
′
+
1
/
2
d
θ
¯
.
The last integral in this equation is an elliptic integral, and this study uses a trapezoidal rule for its evaluation. This numerical approach approximates the valueK
i
-
i
′
,
j
-
j
′
,
k
0as follows:
(32)
K
i
-
i
′
,
j
-
j
′
,
k
0
≈
1
2
(
r
i
H
1
(
r
i
′
+
1
/
2
r
i
,
θ
j
′
+
1
/
2
-
θ
j
,
φ
k
)
-
r
i
H
1
(
r
i
′
-
1
/
2
r
i
,
θ
j
′
-
1
/
2
-
θ
j
,
φ
k
)
+
r
i
H
1
(
r
i
′
+
1
/
2
r
i
,
θ
j
′
-
1
/
2
-
θ
j
,
φ
k
)
-
r
i
H
1
(
r
i
′
-
1
/
2
r
i
,
θ
j
′
-
1
/
2
-
θ
j
,
φ
k
)
)
×
(
θ
j
′
+
1
/
2
-
θ
j
′
-
1
/
2
)
∶
=
r
i
H
1
(
r
¯
r
i
,
θ
¯
-
θ
j
,
φ
k
)
|
r
i
′
-
1
/
2
r
i
′
+
1
/
2
]
θ
j
′
-
1
/
2
θ
j
′
+
1
/
2
,
where the notationf
(
·
)
]
a
b
=
(
1
/
2
)
(
f
(
b
)
+
f
(
a
)
)
(
b
-
a
). Similarly,
(33)
K
i
-
i
′
,
j
-
j
′
,
k
r
≈
r
i
2
(
H
2
(
r
¯
r
i
,
θ
¯
-
θ
j
,
φ
k
)
|
r
i
′
-
1
/
2
r
i
′
+
1
/
2
]
θ
j
′
-
1
/
2
θ
j
′
+
1
/
2
-
r
i
′
r
i
K
i
-
i
′
,
j
-
j
′
,
k
0
)
,
K
i
-
i
′
,
j
-
j
′
,
k
θ
≈
(
θ
¯
-
θ
j
)
r
i
H
1
(
r
¯
r
i
,
θ
¯
-
θ
j
,
φ
k
)
|
r
i
′
-
1
/
2
r
i
′
+
1
/
2
]
θ
j
′
-
1
/
2
θ
j
′
+
1
/
2
.


## 3. Convergence Results

Let us denote|
A
|as the area of the regionAinR
2and discuss the convergence of the proposed method. The integrals of  K
(
x
,
y
,
0
)over a region  R
=
{
(
x
,
y
)
:
x
2
+
y
2
≤
a
}andD
=
{
(
x
,
y
)
:
-
a
≤
x
≤
a
,
-
a
≤
y
≤
a
}, for some constanta, are
(34)
∬
R
K
(
x
,
y
,
0
)
d
x
d
y
=
1
2
a
=
1
2
|
R
|
π
=
C
|
R
|
,
(35)
∬
D
K
(
x
,
y
,
0
)
d
x
d
y
=
y
ln
⁡
(
x
+
x
2
+
y
2
)
+
x
ln
⁡
(
y
+
x
2
+
y
2
)
|
-
a
a
|
-
a
a
=
4
a
ln
⁡
(
3
+
2
2
)
=
C
|
D
|
,
whereCis a constant. Furthermore, ifD
=
R
=
R
2, then this integral is∞.Lemma 1.
If two compact disk densitiesσ
1 and σ
2 satisfy |
σ
1
(
x
,
y
)
-
σ
2
(
x
,
y
)
|
<
h then their corresponding potentials Φ
1 and Φ
2, respectively, satisfy(36)
|
Φ
1
(
x
,
y
)
-
Φ
2
(
x
,
y
)
|
≤
C
h
,
whereCis a positive constant and depends on the supports ofσ
1andσ
2.Proof.
Let the supports ofσ
1andσ
2be contained in a finite regionD. Following (5) and (35), it yields that
(37)
|
Φ
1
(
x
,
y
)
-
Φ
2
(
x
,
y
)
|
=
|
∬
D
K
(
x
¯
-
x
,
y
¯
-
y
,
0
)
(
σ
1
(
x
¯
,
y
¯
)
-
σ
2
(
x
¯
,
y
¯
)
)
d
x
¯
d
y
¯
|
≤
∬
D
K
(
x
¯
-
x
,
y
¯
-
y
,
0
)
|
(
σ
1
(
x
¯
,
y
¯
)
-
σ
2
(
x
¯
,
y
¯
)
)
|
d
x
¯
d
y
¯
≤
C
h
.
Since|
D
|is finite and the lemma is obtained.The order of accuracy between analytic and numerical solutions can be observed from the truncated errors of the disk surfaceσ
(
x
,
y
)and truncated densityσ
n
(
x
,
y
),
(38)
|
σ
n
(
x
,
y
)
-
σ
(
x
,
y
)
|
≤
C
(
Δ
x
)
s
.The potential calculated numerically is denoted byΦ
nand the analytic solution isΦ. We estimate
(39)
|
Φ
n
(
x
,
y
)
-
Φ
(
x
,
y
)
|
≤
∑
i
=
1
M
‍
∑
j
=
1
M
∬
D
i
,
j
‍
K
(
x
¯
-
x
,
y
¯
-
y
,
0
)
×
|
(
σ
1
(
x
¯
,
y
¯
)
-
σ
2
(
x
¯
,
y
¯
)
)
|
d
x
¯
d
y
¯
≤
C
(
Δ
x
)
2
∑
i
=
1
M
‍
∑
j
=
1
M
∬
D
i
,
j
‍
K
(
x
¯
-
x
,
y
¯
-
y
,
0
)
d
x
¯
d
y
¯
≤
C
(
Δ
x
)
s
∬
D
K
(
x
¯
-
x
,
y
¯
-
y
,
0
)
d
x
¯
d
y
¯
≤
C
~
(
Δ
x
)
s
.We deduce the following theorem by aforementioned argument and Lemma1.Theorem 2.
If the truncated errors of the Taylor series of the disk surface densityσ
(
x
,
y
)areO
(
h
s
), wherehis the mesh size in one dimension ands
>
0, then the errors between the numerical and analytic potential also are O
(
h
s
).


## 4. Numerical Results

This study verifies that the proposed method achieved is of second-order accuracy through the following two examples: aD
2disk [9] and a nonaxisymmetric diskD
2,2consisting of two superposedD
2disks. TheD
2disk has the surface density
(40)
Σ
D
2
(
R
;
α
)
=
{
σ
0
(
1
-
R
α
2
)
3
/
2
,
for
R
<
α
,
0
,
for
R
>
α
,
whereR
=
x
2
+
y
2andαis a given constant. The corresponding potential is
(41)
Φ
D
2
(
R
,
z
;
α
)
=
-
π
σ
0
G
64
α
3
[
(
f
1
-
f
2
2
R
)
6
(
8
α
4
-
8
α
2
R
2
+
16
α
2
z
2
+
3
R
4
-
24
R
2
z
2
+
8
z
4
)
sin
-
1
(
f
1
-
f
2
2
R
)
+
2
α
(
18
α
2
-
9
R
2
+
26
z
2
)
f
1
f
2
-
f
3
-
2
|
z
|
(
58
α
2
-
55
R
2
+
50
z
2
)
f
1
f
2
+
f
3
(
f
1
-
f
2
2
R
)
]
,
where
(42)
f
1
=
z
2
+
(
R
+
α
)
2
,
f
2
=
z
2
+
(
R
-
α
)
2
,
f
3
=
α
2
-
R
2
-
z
2
.


### 4.1. Order of Accuracy

This study investigates the order of accuracy of the proposed method in normsL
1,L
2, andL
∞to demonstrate the convergence in total variation, energy, and pointwise senses, respectively. These norms for a functionfon a domainΩare defined as
(43)
∥
f
∥
p
=
(
∫
Ω
‍
|
f
(
x
)
|
p
d
x
)
1
/
p
,
for
p
=
1,2
,
∥
f
∥
∞
=
ess
Ω
sup
⁡
|
f
(
x
)
|
.
Without loss of generality for studying the order of accuracy, set the computational domainΩ
=
[
-
1,1
]
×
[
-
1,1
], and σ
0
=
G
=
1. Setα
=
0.25for aD
2disk simulation andα
=
0.4for aD
2,2disk simulation, whereΣ
D
2,2
=
Σ
D
2
(
R
1
;
α
)
+
Σ
D
2
(
R
2
;
α
)withR
1
=
(
x
-
0.2
)
2
+
y
2andR
2
=
(
x
+
0.2
)
2
+
y
2. Figure 1 presents the profiles of the surface density, the potential onz
=
0plane and the residual between analytic and numerical solutions forN
=
512. In Table 1, columnE
prepresents the error between the analytic and numerical solutions inp-norm,p
=
1,2,   and∞. ColumnO
pin Table 1 represents the order of accuracy as measured bylog
⁡
2
(
E
p
(
2
k
-
1
)
/
E
p
(
2
k
)
)fork
=
6to9and similarly forO
p. These errors show that the proposed method achieves nearly second-order accuracy. More precisely, the order of convergence is approximately1.9or2.0.Table 1
The errors and order accuracy of the proposed method on Cartesian coordinates for various number of zones,N
=
2
k from k
=
5 to 9 for the D
2 (up table) and D
2,2 (down table) disks. This table shows that the accuracy is nearly second order.
N
k
E
1
E
2
E
∞
N
k
-
1
/
N
k
O
1
O
2
O
∞
32
5.617
E
-
3
5.836
E
-
3
1.309
E
-
2
64
1.573
E
-
3
1.555
E
-
3
3.435
E
-
3
32/64
1.836
1.908
1.930
128
3.990
E
-
4
4.022
E
-
4
8.841
E
-
4
64/128
1.979
1.951
1.958
256
1.025
E
-
4
1.023
E
-
4
2.322
E
-
4
128/256
1.961
1.975
1.928
512
2.584
E
-
5
2.578
E
-
5
5.614
E
-
5
256/512
1.988
1.988
2.048
N
k
E
1
E
2
E
∞
N
k
-
1
/
N
k
O
1
O
2
O
∞
32
8.008
E
-
3
7.331
E
-
3
1.822
E
-
2
64
2.129
E
-
3
2.046
E
-
3
5.179
E
-
3
32/64
1.911
1.841
1.815
128
5.344
E
-
4
5.415
E
-
4
1.374
E
-
3
64/128
1.994
1.918
1.914
256
1.396
E
-
4
1.374
E
-
4
3.468
E
-
4
128/256
1.937
1.979
1.986
512
3.516
E
-
5
3.462
E
-
5
8.734
E
-
5
256/512
1.989
1.989
1.989A simulation of theD
2(up row) andD
2,2(down row) disks forN
=
512. The profiles are surface density (left), the potential on the planez
=
0(middle), and the residual between analytic and numerical solutions in natural common logarithm, respectively.
(a)
(b)
(c)
(d)
(e)
(f)Continue to use theD
2disk as an example and a unit diskD
(
0,1
)
=
Ω
=
[
0,1
]
×
[
0,2
π
]as the computational domain to investigate the potential in polar coordinates. Set the valueβ
0
=
0.99. Figure 2 presents the profiles of the surface density, potential, and the difference between analytic and numerical solutions forN
=
512. According to Table 2, the order of accuracy is only approximately 1. Although the surface density at the origin is smooth, the singular elliptic integral introduces significant error there.Table 2
The errors and order accuracy of the proposed method on polar coordinates for various number of zones,N
=
2
k from k
=
5 to 9 for the D
2 disk. This table shows that the accuracy is only first order.
N
k
E
1
E
2
E
∞
N
k
-
1
/
N
k
O
1
O
2
O
∞
32
4.43
E
-
2
3.62
E
-
2
5.25
E
-
2
64
1.38
E
-
2
1.47
E
-
2
2.33
E
-
2
32/64
1.68
1.30
1.17
128
5.89
E
-
3
7.00
E
-
3
1.13
E
-
2
64/128
1.23
1.07
1.04
256
2.77
E
-
3
3.42
E
-
3
5.58
E
-
3
128/256
1.08
1.03
1.02
512
1.48
E
-
3
1.68
E
-
4
2.76
E
-
3
256/512
0.90
1.03
1.02The numerical solutions of theD
2disk forN
=
512to investigate the calculation of potentials in polar coordinates. From left to right, the profiles are surface density, the potential on the planez
=
0, and the difference between analytic and numerical solutions in natural common logarithm, respectively.
(a)
(b)
(c)


### 4.2. A Comparison Study

Start with(44)
Φ
(
x
,
y
,
0
)
=
-
G
∬
K
(
x
¯
-
x
,
y
¯
-
y
,
0
)
σ
(
x
¯
,
y
¯
)
d
x
¯
d
y
¯
=
-
G
∑
i
=
1
N
‍
∑
j
=
1
N
‍
∬
D
i
,
j
K
(
x
¯
-
x
,
y
¯
-
y
,
0
)
σ
(
x
¯
,
y
¯
)
d
x
¯
d
y
¯
.
And introduce a softening parameterϵto approximate
(45)
∬
D
i
,
j
K
(
x
¯
-
x
,
y
¯
-
y
)
σ
(
x
¯
,
y
¯
)
≈
-
G
ϵ
2
+
(
x
i
′
-
x
i
)
2
+
(
y
j
′
-
y
j
)
2
∬
D
i
,
j
σ
(
x
¯
,
y
¯
)
d
x
¯
d
y
¯
.
For polar coordinates [1], the value ofKcan be approximated by
(46)
K
i
′
-
i
,
j
′
-
j
∶
=
-
G
2
(
cosh
⁡
(
u
i
′
-
u
i
)
-
cos
⁡
(
θ
j
′
-
θ
j
)
)
,
whereu
i
′
=
ln
⁡
(
x
i
′
)andu
i
=
ln
⁡
(
x
i
). Note that when(
i
′
,
j
′
)
=
(
i
,
j
),   K  is undefined. Previous research has presented a way to avoid the singularity problem [1]. However, the proposed method avoids the singularity problem by directly evaluating the forces, thereby increasing the order of accuracy.This study compares the proposed method with theN-body method using the simulations of theD
2andD
2,2disks in the previous section. For Cartesian coordinates, choose the softening parameters as the mesh sizeϵ
=
Δ
x. Table 3 shows the errors for the disksD
2 and D
2,2 that the accuracy when using the softening parameter approach for theD
2andD
2,2disks is of first order. This comparison confirms that the proposed method is more accurate and has a higher order of accuracy.Table 3
The top (bottom) table demonstrates the errors and order accuracy of theN-body simulations in Cartesian coordinates for various number of zones, N
=
2
k from k
=
5 to 9 for the D
2 (D
2,2) disk. This table shows that the accuracy is only first order.
N
k
E
1
E
2
L
∞
N
k
-
1
/
N
k
O
1
O
2
O
∞
32
1.199
E
-
1
1.532
E
-
1
3.362
E
-
1
64
6.069
E
-
2
8.141
E
-
2
1.818
E
-
1
32/64
0.9823
0.9121
0.8870
128
3.058
E
-
2
4.203
E
-
2
9.446
E
-
1
64/128
0.9889
0.9538
0.9446
256
1.535
E
-
2
2.136
E
-
2
4.813
E
-
1
128/256
0.9943
0.9765
0.9728
512
7.687
E
-
2
1.077
E
-
2
2.429
E
-
2
256/512
0.9977
0.9879
0.9866
N
k
E
1
E
2
L
∞
N
k
-
1
/
N
k
O
1
O
2
O
∞
32
1.522
E
-
1
1.904
E
-
1
4.227
E
-
1
64
7.762
E
-
2
1.016
E
-
1
2.326
E
-
1
32/64
0.9714
0.9065
0.8617
128
3.916
E
-
2
5.252
E
-
2
1.218
E
-
1
64/128
0.9871
0.9514
0.9333
256
1.964
E
-
2
2.671
E
-
2
6.228
E
-
2
128/256
0.9949
0.9756
0.9676
512
9.839
E
-
3
1.347
E
-
2
3.149
E
-
2
256/512
0.9973
0.9875
0.9839This study calculates potentials on a disk of surface density with as few restrictions as possible. Alternatively, it is possible to solve the reduced equation given by(47)
Φ
(
x
,
y
,
0
)
=
-
G
∬
σ
(
x
¯
,
y
¯
)
(
x
¯
-
x
)
2
+
(
y
¯
-
y
)
2
d
x
¯
d
y
¯
or
(48)
Φ
(
r
,
θ
,
0
)
=
-
G
∬
σ
(
r
¯
,
θ
¯
)
r
¯
2
+
r
2
-
2
r
¯
r
cos
⁡
(
θ
¯
-
θ
)
r
¯
d
r
¯
d
θ
¯
.
The approach in [17] transforms the polar coordinate(
r
,
θ
)into the coordinate(
u
,
θ
)by settingr
=
e
uoru
=
ln
⁡
(
r
). The potential-density pair in terms of the reduced surface density and reduced potential is given in [17], and it is
(49)
e
3
u
/
2
σ
(
e
u
,
θ
)
=
1
4
π
2
∑
m
‍
∫
-
∞
∞
‍
A
m
(
α
)
e
i
(
m
θ
+
α
u
)
d
α
,
(50)
e
u
/
2
Φ
(
e
u
,
θ
)
=
-
1
2
π
G
∑
m
‍
∫
-
∞
∞
‍
K
(
α
,
m
)
A
m
(
α
)
exp
⁡
[
i
(
m
θ
+
α
u
)
]
d
α
,
whereKis real and positive and is defined as
(51)
K
(
α
,
m
)
≡
1
2
Γ
[
(
m
+
1
/
2
+
i
α
)
/
2
]
Γ
[
(
m
+
1
/
2
-
i
α
)
/
2
]
Γ
[
(
m
+
3
/
2
+
i
α
)
/
2
]
Γ
[
(
m
+
3
/
2
-
i
α
)
/
2
]
.
The bounded unit diskD
(
0,1
)
=
[
0,1
]
×
[
0,2
π
]can be transformed to the unbounded domainU
=
(
-
∞
,
0
]
×
[
0,2
π
]and apply this method to theD
2disk using the polar coordinates. In this special case, computem
=
0and truncate
(52)
A
0
(
α
)
=
∫
-
∞
0
‍
e
3
u
/
2
σ
(
e
u
)
e
-
i
α
u
d
u
≈
∫
u
min
⁡
0
‍
e
3
u
/
2
σ
(
e
u
)
e
-
i
α
u
d
u
,
where the valueu
min
⁡is to approximately-
∞. The truncation process produces a hole in the unit disk and can introduce significant errors at the origin. Given a positive integerN and based on the discretization for the radial region in the previous subsection, calculate (52) and (50) using the trapzoidal rule. Figure 3 shows the variation of the potential with respect to radius. The profile on the left panel shows that the numerical and analytic solutions for the Kalnajs method agree well except for those close to the origin ofN
=
512. The small window embedded within the panel zooms in on the residuals between the numerical and analytic solutions to the interval[
0,0.3
]. The truncated portion contributes to significant errors near the origin. In contrast, the application of the proposed method to the calculation of potentials leads to the results shown in the right panel of Figure 3. Although the singular integration remains because of the unbounded domain, the proposed method is preferable for both Cartesian and polar coordinates because a hole near the origin is not introduced.The variation of the potential with respect to radius using the Kalnajs method (a) and the proposed method (b). The residuals appear in the small window in each panel and show that the Kalnajs method produces significant errors near the origin, which are eliminated by the proposed method.
(a)
(b)


### 4.3. Potentials onz-Axis

Although this study concentrates on the calculation of potentials for a disk, the proposed method can also calculate the forces on the whole space compared with [18]. The simulation of Example1in Section 4.1 shows that the order of accuracy for thez-axis is of nearly second order (see Table 4), and Figure 4 shows the profiles, the potential onz-axis and the residual between analytic and numerical solutions.Table 4
The errors of the potential on thez-axis and order accuracy of the proposed method on Cartesian coordinates for various zones, N
=
2
k from k
=
5 to 9 for the D
2 disk. This table shows that the accuracy is nearly second order.
N
k
E
1
(
z
)
E
2
(
z
)
E
∞
(
z
)
N
k
-
1
/
N
k
O
1
O
2
O
∞
32
1.175
E
-
2
1.129
E
-
2
1.518
E
-
2
64
2.998
E
-
3
3.093
E
-
3
5.230
E
-
3
32/64
1.971
1.868
1.537
128
7.672
E
-
4
7.963
E
-
4
1.541
E
-
4
64/128
1.966
1.958
1.958
256
1.910
E
-
4
2.000
E
-
4
4.171
E
-
4
128/256
2.006
1.993
1.993
512
4.775
E
-
5
5.006
E
-
5
1.085
E
-
4
256/512
2.000
1.988
1.998The simulation of aD
2disk forN
=
512. The profiles are the potential onz-axis (a) and the residuals between analytic and numerical solutions (b) in natural common logarithm.
(a)
(b)


## 4.1. Order of Accuracy

This study investigates the order of accuracy of the proposed method in normsL
1,L
2, andL
∞to demonstrate the convergence in total variation, energy, and pointwise senses, respectively. These norms for a functionfon a domainΩare defined as
(43)
∥
f
∥
p
=
(
∫
Ω
‍
|
f
(
x
)
|
p
d
x
)
1
/
p
,
for
p
=
1,2
,
∥
f
∥
∞
=
ess
Ω
sup
⁡
|
f
(
x
)
|
.
Without loss of generality for studying the order of accuracy, set the computational domainΩ
=
[
-
1,1
]
×
[
-
1,1
], and σ
0
=
G
=
1. Setα
=
0.25for aD
2disk simulation andα
=
0.4for aD
2,2disk simulation, whereΣ
D
2,2
=
Σ
D
2
(
R
1
;
α
)
+
Σ
D
2
(
R
2
;
α
)withR
1
=
(
x
-
0.2
)
2
+
y
2andR
2
=
(
x
+
0.2
)
2
+
y
2. Figure 1 presents the profiles of the surface density, the potential onz
=
0plane and the residual between analytic and numerical solutions forN
=
512. In Table 1, columnE
prepresents the error between the analytic and numerical solutions inp-norm,p
=
1,2,   and∞. ColumnO
pin Table 1 represents the order of accuracy as measured bylog
⁡
2
(
E
p
(
2
k
-
1
)
/
E
p
(
2
k
)
)fork
=
6to9and similarly forO
p. These errors show that the proposed method achieves nearly second-order accuracy. More precisely, the order of convergence is approximately1.9or2.0.Table 1
The errors and order accuracy of the proposed method on Cartesian coordinates for various number of zones,N
=
2
k from k
=
5 to 9 for the D
2 (up table) and D
2,2 (down table) disks. This table shows that the accuracy is nearly second order.
N
k
E
1
E
2
E
∞
N
k
-
1
/
N
k
O
1
O
2
O
∞
32
5.617
E
-
3
5.836
E
-
3
1.309
E
-
2
64
1.573
E
-
3
1.555
E
-
3
3.435
E
-
3
32/64
1.836
1.908
1.930
128
3.990
E
-
4
4.022
E
-
4
8.841
E
-
4
64/128
1.979
1.951
1.958
256
1.025
E
-
4
1.023
E
-
4
2.322
E
-
4
128/256
1.961
1.975
1.928
512
2.584
E
-
5
2.578
E
-
5
5.614
E
-
5
256/512
1.988
1.988
2.048
N
k
E
1
E
2
E
∞
N
k
-
1
/
N
k
O
1
O
2
O
∞
32
8.008
E
-
3
7.331
E
-
3
1.822
E
-
2
64
2.129
E
-
3
2.046
E
-
3
5.179
E
-
3
32/64
1.911
1.841
1.815
128
5.344
E
-
4
5.415
E
-
4
1.374
E
-
3
64/128
1.994
1.918
1.914
256
1.396
E
-
4
1.374
E
-
4
3.468
E
-
4
128/256
1.937
1.979
1.986
512
3.516
E
-
5
3.462
E
-
5
8.734
E
-
5
256/512
1.989
1.989
1.989A simulation of theD
2(up row) andD
2,2(down row) disks forN
=
512. The profiles are surface density (left), the potential on the planez
=
0(middle), and the residual between analytic and numerical solutions in natural common logarithm, respectively.
(a)
(b)
(c)
(d)
(e)
(f)Continue to use theD
2disk as an example and a unit diskD
(
0,1
)
=
Ω
=
[
0,1
]
×
[
0,2
π
]as the computational domain to investigate the potential in polar coordinates. Set the valueβ
0
=
0.99. Figure 2 presents the profiles of the surface density, potential, and the difference between analytic and numerical solutions forN
=
512. According to Table 2, the order of accuracy is only approximately 1. Although the surface density at the origin is smooth, the singular elliptic integral introduces significant error there.Table 2
The errors and order accuracy of the proposed method on polar coordinates for various number of zones,N
=
2
k from k
=
5 to 9 for the D
2 disk. This table shows that the accuracy is only first order.
N
k
E
1
E
2
E
∞
N
k
-
1
/
N
k
O
1
O
2
O
∞
32
4.43
E
-
2
3.62
E
-
2
5.25
E
-
2
64
1.38
E
-
2
1.47
E
-
2
2.33
E
-
2
32/64
1.68
1.30
1.17
128
5.89
E
-
3
7.00
E
-
3
1.13
E
-
2
64/128
1.23
1.07
1.04
256
2.77
E
-
3
3.42
E
-
3
5.58
E
-
3
128/256
1.08
1.03
1.02
512
1.48
E
-
3
1.68
E
-
4
2.76
E
-
3
256/512
0.90
1.03
1.02The numerical solutions of theD
2disk forN
=
512to investigate the calculation of potentials in polar coordinates. From left to right, the profiles are surface density, the potential on the planez
=
0, and the difference between analytic and numerical solutions in natural common logarithm, respectively.
(a)
(b)
(c)


## 4.2. A Comparison Study

Start with(44)
Φ
(
x
,
y
,
0
)
=
-
G
∬
K
(
x
¯
-
x
,
y
¯
-
y
,
0
)
σ
(
x
¯
,
y
¯
)
d
x
¯
d
y
¯
=
-
G
∑
i
=
1
N
‍
∑
j
=
1
N
‍
∬
D
i
,
j
K
(
x
¯
-
x
,
y
¯
-
y
,
0
)
σ
(
x
¯
,
y
¯
)
d
x
¯
d
y
¯
.
And introduce a softening parameterϵto approximate
(45)
∬
D
i
,
j
K
(
x
¯
-
x
,
y
¯
-
y
)
σ
(
x
¯
,
y
¯
)
≈
-
G
ϵ
2
+
(
x
i
′
-
x
i
)
2
+
(
y
j
′
-
y
j
)
2
∬
D
i
,
j
σ
(
x
¯
,
y
¯
)
d
x
¯
d
y
¯
.
For polar coordinates [1], the value ofKcan be approximated by
(46)
K
i
′
-
i
,
j
′
-
j
∶
=
-
G
2
(
cosh
⁡
(
u
i
′
-
u
i
)
-
cos
⁡
(
θ
j
′
-
θ
j
)
)
,
whereu
i
′
=
ln
⁡
(
x
i
′
)andu
i
=
ln
⁡
(
x
i
). Note that when(
i
′
,
j
′
)
=
(
i
,
j
),   K  is undefined. Previous research has presented a way to avoid the singularity problem [1]. However, the proposed method avoids the singularity problem by directly evaluating the forces, thereby increasing the order of accuracy.This study compares the proposed method with theN-body method using the simulations of theD
2andD
2,2disks in the previous section. For Cartesian coordinates, choose the softening parameters as the mesh sizeϵ
=
Δ
x. Table 3 shows the errors for the disksD
2 and D
2,2 that the accuracy when using the softening parameter approach for theD
2andD
2,2disks is of first order. This comparison confirms that the proposed method is more accurate and has a higher order of accuracy.Table 3
The top (bottom) table demonstrates the errors and order accuracy of theN-body simulations in Cartesian coordinates for various number of zones, N
=
2
k from k
=
5 to 9 for the D
2 (D
2,2) disk. This table shows that the accuracy is only first order.
N
k
E
1
E
2
L
∞
N
k
-
1
/
N
k
O
1
O
2
O
∞
32
1.199
E
-
1
1.532
E
-
1
3.362
E
-
1
64
6.069
E
-
2
8.141
E
-
2
1.818
E
-
1
32/64
0.9823
0.9121
0.8870
128
3.058
E
-
2
4.203
E
-
2
9.446
E
-
1
64/128
0.9889
0.9538
0.9446
256
1.535
E
-
2
2.136
E
-
2
4.813
E
-
1
128/256
0.9943
0.9765
0.9728
512
7.687
E
-
2
1.077
E
-
2
2.429
E
-
2
256/512
0.9977
0.9879
0.9866
N
k
E
1
E
2
L
∞
N
k
-
1
/
N
k
O
1
O
2
O
∞
32
1.522
E
-
1
1.904
E
-
1
4.227
E
-
1
64
7.762
E
-
2
1.016
E
-
1
2.326
E
-
1
32/64
0.9714
0.9065
0.8617
128
3.916
E
-
2
5.252
E
-
2
1.218
E
-
1
64/128
0.9871
0.9514
0.9333
256
1.964
E
-
2
2.671
E
-
2
6.228
E
-
2
128/256
0.9949
0.9756
0.9676
512
9.839
E
-
3
1.347
E
-
2
3.149
E
-
2
256/512
0.9973
0.9875
0.9839This study calculates potentials on a disk of surface density with as few restrictions as possible. Alternatively, it is possible to solve the reduced equation given by(47)
Φ
(
x
,
y
,
0
)
=
-
G
∬
σ
(
x
¯
,
y
¯
)
(
x
¯
-
x
)
2
+
(
y
¯
-
y
)
2
d
x
¯
d
y
¯
or
(48)
Φ
(
r
,
θ
,
0
)
=
-
G
∬
σ
(
r
¯
,
θ
¯
)
r
¯
2
+
r
2
-
2
r
¯
r
cos
⁡
(
θ
¯
-
θ
)
r
¯
d
r
¯
d
θ
¯
.
The approach in [17] transforms the polar coordinate(
r
,
θ
)into the coordinate(
u
,
θ
)by settingr
=
e
uoru
=
ln
⁡
(
r
). The potential-density pair in terms of the reduced surface density and reduced potential is given in [17], and it is
(49)
e
3
u
/
2
σ
(
e
u
,
θ
)
=
1
4
π
2
∑
m
‍
∫
-
∞
∞
‍
A
m
(
α
)
e
i
(
m
θ
+
α
u
)
d
α
,
(50)
e
u
/
2
Φ
(
e
u
,
θ
)
=
-
1
2
π
G
∑
m
‍
∫
-
∞
∞
‍
K
(
α
,
m
)
A
m
(
α
)
exp
⁡
[
i
(
m
θ
+
α
u
)
]
d
α
,
whereKis real and positive and is defined as
(51)
K
(
α
,
m
)
≡
1
2
Γ
[
(
m
+
1
/
2
+
i
α
)
/
2
]
Γ
[
(
m
+
1
/
2
-
i
α
)
/
2
]
Γ
[
(
m
+
3
/
2
+
i
α
)
/
2
]
Γ
[
(
m
+
3
/
2
-
i
α
)
/
2
]
.
The bounded unit diskD
(
0,1
)
=
[
0,1
]
×
[
0,2
π
]can be transformed to the unbounded domainU
=
(
-
∞
,
0
]
×
[
0,2
π
]and apply this method to theD
2disk using the polar coordinates. In this special case, computem
=
0and truncate
(52)
A
0
(
α
)
=
∫
-
∞
0
‍
e
3
u
/
2
σ
(
e
u
)
e
-
i
α
u
d
u
≈
∫
u
min
⁡
0
‍
e
3
u
/
2
σ
(
e
u
)
e
-
i
α
u
d
u
,
where the valueu
min
⁡is to approximately-
∞. The truncation process produces a hole in the unit disk and can introduce significant errors at the origin. Given a positive integerN and based on the discretization for the radial region in the previous subsection, calculate (52) and (50) using the trapzoidal rule. Figure 3 shows the variation of the potential with respect to radius. The profile on the left panel shows that the numerical and analytic solutions for the Kalnajs method agree well except for those close to the origin ofN
=
512. The small window embedded within the panel zooms in on the residuals between the numerical and analytic solutions to the interval[
0,0.3
]. The truncated portion contributes to significant errors near the origin. In contrast, the application of the proposed method to the calculation of potentials leads to the results shown in the right panel of Figure 3. Although the singular integration remains because of the unbounded domain, the proposed method is preferable for both Cartesian and polar coordinates because a hole near the origin is not introduced.The variation of the potential with respect to radius using the Kalnajs method (a) and the proposed method (b). The residuals appear in the small window in each panel and show that the Kalnajs method produces significant errors near the origin, which are eliminated by the proposed method.
(a)
(b)


## 4.3. Potentials onz-Axis

Although this study concentrates on the calculation of potentials for a disk, the proposed method can also calculate the forces on the whole space compared with [18]. The simulation of Example1in Section 4.1 shows that the order of accuracy for thez-axis is of nearly second order (see Table 4), and Figure 4 shows the profiles, the potential onz-axis and the residual between analytic and numerical solutions.Table 4
The errors of the potential on thez-axis and order accuracy of the proposed method on Cartesian coordinates for various zones, N
=
2
k from k
=
5 to 9 for the D
2 disk. This table shows that the accuracy is nearly second order.
N
k
E
1
(
z
)
E
2
(
z
)
E
∞
(
z
)
N
k
-
1
/
N
k
O
1
O
2
O
∞
32
1.175
E
-
2
1.129
E
-
2
1.518
E
-
2
64
2.998
E
-
3
3.093
E
-
3
5.230
E
-
3
32/64
1.971
1.868
1.537
128
7.672
E
-
4
7.963
E
-
4
1.541
E
-
4
64/128
1.966
1.958
1.958
256
1.910
E
-
4
2.000
E
-
4
4.171
E
-
4
128/256
2.006
1.993
1.993
512
4.775
E
-
5
5.006
E
-
5
1.085
E
-
4
256/512
2.000
1.988
1.998The simulation of aD
2disk forN
=
512. The profiles are the potential onz-axis (a) and the residuals between analytic and numerical solutions (b) in natural common logarithm.
(a)
(b)


## 5. Conclusion

This study presents an improved method of theN-body calculation for solving Poisson equation induced by a disk of surface density. The proposed method does not require artificial boundary conditions and it also does not require a softening length parameter and the computational complexity is O
(
N
2
(
log
⁡
2
(
N
)
2
)
), which is the same order as theN-body method. The proposed method also achieves second-order accuracy.



---

*Source: 693537-2014-07-24.xml*