Journal of Applied Mathematics

Volume 2014 (2014), Article ID 693537, 12 pages

http://dx.doi.org/10.1155/2014/693537

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

Department of Mathematics, Fu Jen Catholic University, New Taipei City 24205, Taiwan

Received 22 April 2014; Accepted 14 July 2014; Published 24 July 2014

Academic Editor: Chein-Shan Liu

Copyright © 2014 Chien-Chang Yen. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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 the -body method and artificial boundary conditions for the spectral methods. With the help of fast Fourier transform, the computational complexity is only , where 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.

#### 1. Introduction

Let us consider the Poisson equation with singular disk surface density [1–3] where,,, and are position, gravitational constant, Dirac symbol, and density, respectively. Hereafter, we setwithout 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 where,, andfor, based on the uniform mesh gridswith the mesh sizes,, andin,, anddirections, 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 involveunknowns, whereis the number of zones in one dimension. The linear complexity of this approach is at least, 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,, where as The difficulties encountered in the numerical approach for solving (1) or (5) are related to the extent of the domain. 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 the-body method [1] has the total amount of complexity based on the numberof mesh zones. For the-body method based on a uniform grid and using the FFT technique, the numerical complexity can be reduced fromto. 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 wherecontains the origin.

The proposed method does not require the artificial periodic boundary condition and softening parameter and further improves the accuracy of the-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, adisk [9] and its variation, 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 of-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. Sections 2.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 disks[9]. Section 4.2 presents a comparison with the-body calculation. Finally, we conclude this study in Section 5.

#### 2. Methods

##### 2.1. Cartesian Coordinates

Assume that the numerical computation domainfor some number, which contains the support of the surface density, and further discretize the region uniformly as follows. Given a positive integer, define,,, and, where. Next, define the center of the cellasand, where. Hence, the domainis discretized into thecells.

We recall the-body method with fast Fourier transform (FFT). Represent the potential functionof (1) as Where According to (2), the potential is Calculate (8) using a numerical approach for the planeassuming that the support of the surface density is compact. Based on a restriction to the calculation of potentials on theplane, derive the formulae for all. Equation (8) is equal to The potential atfor the-body method is approximated by That is, The approximated integral is of the first order and the fundamental function is singular whenever. The-body method applies a softening parameter technique forwith, to avoid the singularity.

To improve the calculation of (11), approximate the surface densityonin (8) linearly by whereandandare constant in the cell. The error of the discretization is. If a higher order accuracy is required, additional terms in the Taylor expansion in (13) can be considered.

Let By substituting (13) into (9), the potential can be approximated by where

Fortunately, the integrals of (14) can be obtained analytically with the help of the following simple formulae: The valueis then equal to where the notation. The calculations ofandare split into two parts by the identity and, respectively. This leads to

The forms of,, andin (16) are a type of convolution and the computational complexity can be reduced towith the help of FFT in theanddirections. In thedirection, the computational complexity is. This implies that the total complexity of this method isandfor the planeand 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, which can be reduced to polar coordinates (). 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 than-body calculations and the method in [17].

The support of the surface density is compact, contained in a region for some number. 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 integer, define,,,,,,, andwhere. The effect ofis 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 is discretized into thecells,for. For, the cellsdo not cover the origin and extra cellsshould be included. To simplify notation, denote,.

The potential functionof (1) under the assumptionin spherical coordinates can be expressed as where

The surface densityonin (20) can be linearly approximated by whereandandare constants in the cell. Equation (22) is the Taylor expansion in two dimensions and it follows that the error of the discretization is.

In this section, let Similarly for,, andon the central region . The termin (24) is for the formulation of a convolution type. According to (20) and (22), where Equation (27) can then be rewritten as

Define, whereis the dimensionless radius. The evaluation of (23), (24), and (25) can be obtained with the help of the following simple integrals:

The definitions ofandlead to Calculate the value of the integral 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 valueas follows: where the notation. Similarly,

#### 3. Convergence Results

Let us denoteas the area of the regioninand discuss the convergence of the proposed method. The integrals of over a region and, for some constant, are whereis a constant. Furthermore, if, then this integral is.

Lemma 1. *
If two compact disk densities and satisfy then their corresponding potentials and , respectively, satisfy**
whereis a positive constant and depends on the supports ofand.*

*Proof. *Let the supports ofandbe contained in a finite region. Following (5) and (35), it yields that
Sinceis 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 surfaceand truncated density,

The potential calculated numerically is denoted byand the analytic solution is. We estimate

We deduce the following theorem by aforementioned argument and Lemma 1.

Theorem 2. *If the truncated errors of the Taylor series of the disk surface densityare, whereis the mesh size in one dimension and, then the errors between the numerical and analytic potential also are .*

#### 4. Numerical Results

This study verifies that the proposed method achieved is of second-order accuracy through the following two examples: adisk [9] and a nonaxisymmetric diskconsisting of two superposeddisks. Thedisk has the surface density whereandis a given constant. The corresponding potential is where

##### 4.1. Order of Accuracy

This study investigates the order of accuracy of the proposed method in norms,, andto demonstrate the convergence in total variation, energy, and pointwise senses, respectively. These norms for a functionon a domainare defined as Without loss of generality for studying the order of accuracy, set the computational domain, and . Setfor adisk simulation andfor adisk simulation, wherewithand. Figure 1 presents the profiles of the surface density, the potential onplane and the residual between analytic and numerical solutions for. In Table 1, columnrepresents the error between the analytic and numerical solutions in-norm,, and. Columnin Table 1 represents the order of accuracy as measured byfortoand similarly for. These errors show that the proposed method achieves nearly second-order accuracy. More precisely, the order of convergence is approximatelyor.

Continue to use thedisk as an example and a unit diskas the computational domain to investigate the potential in polar coordinates. Set the value. Figure 2 presents the profiles of the surface density, potential, and the difference between analytic and numerical solutions for. 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.

##### 4.2. A Comparison Study

Start with And introduce a softening parameterto approximate For polar coordinates [1], the value ofcan be approximated by whereand. Note that when, 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 the-body method using the simulations of theanddisks in the previous section. For Cartesian coordinates, choose the softening parameters as the mesh size. Table 3 shows the errors for the disks and that the accuracy when using the softening parameter approach for theanddisks is of first order. This comparison confirms that the proposed method is more accurate and has a higher order of accuracy.

This 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 or The approach in [17] transforms the polar coordinateinto the coordinateby settingor. The potential-density pair in terms of the reduced surface density and reduced potential is given in [17], and it is whereis real and positive and is defined as The bounded unit diskcan be transformed to the unbounded domainand apply this method to thedisk using the polar coordinates. In this special case, computeand truncate where the valueis to approximately. The truncation process produces a hole in the unit disk and can introduce significant errors at the origin. Given a positive integer 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 of. The small window embedded within the panel zooms in on the residuals between the numerical and analytic solutions to the interval. 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.

##### 4.3. Potentials on -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 Examplein Section 4.1 shows that the order of accuracy for the-axis is of nearly second order (see Table 4), and Figure 4 shows the profiles, the potential on-axis and the residual between analytic and numerical solutions.

#### 5. Conclusion

This study presents an improved method of the-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 , which is the same order as the-body method. The proposed method also achieves second-order accuracy.

#### Conflict of Interests

The author declares that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

The author is grateful to anonymous referees for their valuable comments. This work was supported in part by the Ministry of Science and Technology, Taiwan, under the Grant NSC 101-2115-M-030-004.

#### References

- J. Binney and S. Tremaine,
*Galactic Dynamics*, Princeton Series in Astrophysics, 2nd edition, 2008. - C. Yuan and D. C. C. Yen, “Evolution of self-gravitating gas disks under the influence of a rotating bar potential,”
*Journal of the Korean Astronomical Society*, vol. 38, no. 2, pp. 197–201, 2005. View at Google Scholar - H. Zhang, C. Yuan, D. N. C. Lin, and D. C. C. Yen, “On the orbital evolution of a Jovian planet embedded in a self-gravitating disk,”
*Astrophysical Journal Letters*, vol. 676, no. 1, pp. 639–650, 2008. View at Publisher · View at Google Scholar · View at Scopus - G. D. Smith,
*Numerical Solution of Partial Differential Equations: Finite Difference Methods*, Oxford University Press, 1986. View at MathSciNet - J. C. Strikwerda,
*Finite Difference Schemes and Partial Differential Equations*, SIAM, 2004. View at Publisher · View at Google Scholar · View at MathSciNet - D. Y. Kwak and J. S. Lee, “Multigrid analysis for higher order finite difference scheme,”
*Journal of Numerical Mathematics*, vol. 12, no. 4, pp. 285–296, 2004. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet · View at Scopus - J. Zhang, “Fast and high accuracy multigrid solution of the three-dimensional Poisson equation,”
*Journal of Computational Physics*, vol. 143, no. 2, pp. 449–461, 1998. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet · View at Scopus - D. Li and F. J. Hickernell, “Trigonometric spectral collocation methods on lattices,” in
*Recent Advances in Scientific Computing and Partial Differential Equations*, S. Y. Cheng, C.-W. Shu, and T. Tang, Eds., vol. 330 of*Contemporary Mathematics*, pp. 121–132, American Mathematical Society, Providence, RI, USA, 2003. View at Publisher · View at Google Scholar · View at MathSciNet - E. Schulz, “Potential-density pairs for a family of finite disks,”
*Astrophysical Journal*, vol. 693, no. 2, pp. 1310–1315, 2009. View at Publisher · View at Google Scholar · View at Scopus - Z. Gimbutas, L. Greengard, and M. Minion, “Coulomb interactions on planar structures: inverting the square root of the Laplacian,”
*SIAM Journal on Scientific Computing*, vol. 22, no. 6, pp. 2093–2108, 2001. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus - J. Carrier, L. Greengard, and V. Rokhlin, “A fast adaptive multipole algorithm for particle simulations,”
*Society for Industrial and Applied Mathematics*, vol. 9, no. 4, pp. 669–686, 1988. View at Publisher · View at Google Scholar · View at MathSciNet - L. Greengard,
*The Rapid Evaluation of Potential Fields in Particle Systems*, ACM Distinguished Dissertations, MIT Press, Cambridge, Mass, USA, 1988. View at MathSciNet - L. Greengard and V. Rokhlin, “A new version of the fast multipole method for the Laplace equation in three dimensions,”
*Acta Numerica*, vol. 6, pp. 229–269, 1997. View at Google Scholar - A. P. Boss, “Protostellar formation in rotating interstellar clouds. I. Numerical methods and tests,”
*The Astrophysical Journal*, vol. 236, pp. 619–627, 1980. View at Google Scholar - L. D. G. Sigalotti, “Protostellar collapse and fragmentation: describing and testing a second-order accurate radiation hydrodynamic code,”
*The Astrophysical Journal Supplement Series*, vol. 116, no. 1, pp. 75–101, 1998. View at Publisher · View at Google Scholar · View at Scopus - L. D. G. Sigalotti and J. Klapp, “Protostellar collapse models of prolate molecular cloud cores,”
*Astronomy and Astrophysics*, vol. 378, no. 1, pp. 165–179, 2001. View at Publisher · View at Google Scholar · View at Scopus - A. J. Kalnajs, “Dynamics of flat galaxies-I,”
*The Astrophysical Journal*, vol. 166, pp. 275–293, 1971. View at Publisher · View at Google Scholar · View at MathSciNet - C. Yen, R. E. Taam, and K. H. a. Yeh, “Self-gravitational force calculation of infinitesimally thin gaseous disks,”
*Journal of Computational Physics*, vol. 231, no. 24, pp. 8246–8261, 2012. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus