A numerical procedure to solve Poisson’s equation in spherical coordinates

Resumo

This paper sets out to present a numerical procedure that solves Poisson’s equation in a spherical coordinate system. To discretize this equation, integration techniques at the interfaces between different regions have been carried out allowing the calculation of both the potential and the corresponding field inside and outside a charge distribution. The Gauss-Seidel method is adopted to determine the potential in each region and the results, whenever compared with the analytical solutions found in the literature, come out very satisfactory, with errors less than 1% for distances of the order of 1 × 10⁻¹⁴ m and, for larger distances, they never reach 4%.

keywords
Charge distribution; Gauss-Seidel method; electrostatic potential; field; numerical solution

1. Introduction

Poisson’s equation is an elliptic partial differential equation with a known non-trivial source term. This equation has a wide application in several areas of Physics and Engineering, such as Electrodynamics, Mechanics, Fluid Dynamics and the study of topological deffects. In Mechanics, for example, it is used to study the gravitational potential of mass distributions. In this case, it is referred to as Poisson’s equation for gravity. Similarly, in Maxwellian Electromagnetism, it allows the computation of the scalar and vector potentials generated by charge and current distributions over space-time.

In the specific literature, several analytic solutions to Poisson’s equation may be found in different coordinate systems. The simplest cases are the one-and two-dimensional systems described in Cartesian coordinates. It is however important to mention that, in other coordinate systems, such as cylindrical or spherical coordinates, analytic solutions may not be found for generic source distributions. However, numerical solutions may always be attained with the help of some specific methods. Different numerical techniques to solve Poisson’s equation to obtain electrostatic potentials are found in the literature of Computational Physics and Applied Mathematics, which use Fourier series [¹], approximate analytic solutions [²] and finite difference discretization scheme [³]. In many situations, to simplify the problem, one assumes a homogeneous medium and the absence of electric charges, which reduces the problem to solving the Laplace equation.

This paper focus on solving Poisson’s equation for problems that involve a uniform spherical charge distribution. To do that, azimuthal symmetry with respect to the ϕ-coordinate is invoked, yielding the electrostatic potential as a function of the r- and θ-coordinates, V(r, θ). Next, by taking the gradient of the latter, one readily computes the electrostatic field over all the space, inside and outside the charge distribution.

In addition, we next contemplate a non-homogeneous medium, where more than one region is involved. The best method we have found to tackle this problem consists in using integration techniques, which means discretizing Poisson’s equation and integrating it numerically over an arbitrary volume element. As a step further, in this method, we can perform an integration on an interface region between the distinct media to have a complete answer.

This technique, known as point- (or interface-) centered scheme [⁴], allows an arbitrary number of regions to be considered, as long as they are physically acceptable. This is possible by using Uniqueness theorem [⁵], which states that the potential and its derivative are both continuous at the interface between two different media. (We should however notice that the presence of a superficial density of charges at the interface would prevent us from using the continuity of the derivative of the potential). As an outcome, the numerical procedure to solve Poisson’s equation we present here, namely, based on the use of integration techniques, opens up the possibility to carry out field calculation for different geometries.

2. Poisson’s Equation

The Poisson’s equation for the electrostatic potential in spherical coordinates with azimuthal symmetry, so that the potential is ϕ-independent [⁵], is written as

(1) $$\begin{array}{l}\begin{array}{l}\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial }{\partial r}V\left(r,\theta \right)\right)\hfill \\ +\frac{1}{r^2\sin \theta }\frac{\partial }{\partial \theta }\left(\sin \theta \frac{\partial }{\partial \theta }V\left(r,\theta \right)\right)=\frac{\rho \left(r\right)}{{\epsilon }_o},\hfill \end{array}\hfill \end{array}$$

where ρ(r) is the charge density distribution, ∊_o is the electric permittivity of the vacuum and V(r, θ) is the potential we are searching for. The regions of interest are bounded as follows: r ≥ 0 and 0° θ ≤ 180°. An analytic solution may be obtained to Eq. (1), if it is assumed that the potential is only r-dependent, that is, the potential is a function V(r). Thus, the angular dependence represents the projection of the electric potential radial solution for different angles within the domain [0°, 180°]. In this paper, an analytical solution is going to be derived in order to validate the numerical method that we shall implement further on.

2.1. The analytic solution

The Poisson’s equation written in spherical coordinates for the r-dependent electrostatic potential reads as follows below:

(2) $$\frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{d}{dr}V\left(r\right)\right)=-\frac{\rho \left(r\right)}{{\epsilon }_o}.$$

The purpose is to obtain solutions inside, (r ≤ R), and outside, (r > R), the charged sphere. The analytic solution to Eq. (2) is obtained by using differential equation techniques, which leads to the solutions in both regions. To do that, we take the charge density as given in what follows:

(3) $$\rho \left(r\right)=\{\begin{array}{lll}{\rho }_o\hfill & \text{if}\hfill & r\le R\hfill \\ 0,\hfill & \text{if}\hfill & r>R\hfill \end{array}$$

Inside the charged sphere, where the density is constant, the electric potential can be obtained by Eq. (2) upon a direct integration, which results in

(4) $$V_{in}\left(r\right)=-\frac{1}{6}\frac{{\rho }_o}{{\epsilon }_o}{r}^2-\frac{A}{r}+B,$$

where V_in(r) is the solution for r ≤ R. Similarly, for the external region, where the charge density is zero, the solution takes the form

(5) $$V_{out}\left(r\right)=C+\frac{D}{r}.$$

where V_out(r) is the potential for r > R. The coefficients of these solutions are fixed by imposing the following conditions:

1. The potential is continuous in R, such that

   (6a) $$V_{in}\left(R\right)=V_{out}\left(R\right).$$

2. The derivative of the potential is continuous in R,

   (6b) $${\frac{d}{dr}{V}_{in}\left(r\right)|}_{r=R}={\frac{d}{dr}{V}_{out}\left(r\right)|}_{r=R}.$$

3. The potential is zero at infinity, such that

   (6c) $$\underset{r\to \infty }{\lim }{V}_{out}\left(r\right)=0.$$

4. The derivative of the potential vanishes at r = 0,

   (6d) $${\frac{d}{dr}{V}_{in}\left(r\right)|}_{r=0}=0,$$

where have

(7) $$V\left(r\right)=\{\begin{array}{lll}\frac{1}{2}\frac{{\rho }_o}{{\epsilon }_o}\left(R^2-\frac{r^2}{3}\right),\hfill & \text{if}\hfill & r\le R\hfill \\ \frac{1}{3}\frac{{\rho }_o}{{\epsilon }_o}\frac{R^3}{r},\hfill & \text{if}\hfill & r>R\hfill \end{array}$$

The electric field E(r) can then be readily attained:

(8) $$\text{E}\left(r\right)=-\nabla V\left(r\right),$$

finally,

(9) $$\text{E}\left(r\right)=\{\begin{array}{lll}\frac{1}{3}\frac{{\rho }_o}{{\epsilon }_o}r\widehat{\text{r}},\hfill & \text{if}\hfill & r\le R\hfill \\ \frac{1}{3}\frac{{\rho }_o}{{\epsilon }_o}\frac{R^3}{r^2}\widehat{\text{r}},\hfill & \text{if}\hfill & r>R\hfill \end{array}$$

These analytic solutions shall be, later on, compared with the numerical solutions, which are going to be worked out in the next Section.

2.2. Getting the numerical solution

The numerical solution to the electric potential in spherical coordinates is now presented. For this, the Finite difference method [⁶], to discretize the Eq. (1) using techniques of integration at the interface, has been used. This procedure approximates the derivative of a function, f(u), for consecutive points [u, u + Δu] as given below:

(10) $$\frac{d}{du}f\left(u\right)=\frac{f\left(u+\Delta u\right)-f\left(u\right)}{\Delta u}.$$

This allows us to approximate the function f(u) in the interval with a truncation error, Θ(Δu), where u = [r, θ]. So, the smaller the interval, more negligeable will be the error, rendering the numerical solution considerably close to the analytical solution. In the limit in which Δu tends to zero in Eq. (10), we have the exact definition of the derivative of a function.

As this system involves different types of regions, the best form to approach the problem is by integrating numerically Eq. (1) in the volume of a sphere in the interval $\left[r_{i-\frac{1}{2}},r_{i+\frac{1}{2}}\right]$ and $\left[{\theta }_{j-\frac{1}{2}},{\theta }_{j+\frac{1}{2}}\right]$, as shown in Figure 1, such that

Figure 1
Different regions of the problem.

(11) $$\begin{array}{l}{\displaystyle {\int }_{{\theta }_{j-\frac{1}{2}}}^{{\theta }_{j+\frac{1}{2}}}{\displaystyle {\int }_{r_{i-\frac{1}{2}}}^{r_{i+\frac{1}{2}}}\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial }{\partial r}V\left(r,\theta \right)\right)2\pi r^2\sin \theta drd\theta }}\hfill \\ +{\displaystyle {\int }_{{\theta }_{j-\frac{1}{2}}}^{{\theta }_{j+\frac{1}{2}}}{\displaystyle {\int }_{r_{i-\frac{1}{2}}}^{r_{i+\frac{1}{2}}}\frac{1}{r^2\sin \theta }\frac{\partial }{\partial \theta }\left(\sin \theta \frac{\partial }{\partial \theta }V\left(r,\theta \right)\right)2\pi r^2\sin \theta drd\theta }}\hfill \\ =-\frac{1}{{\epsilon }_o}{\displaystyle {\int }_{{\theta }_{j-\frac{1}{2}}}^{{\theta }_{j+\frac{1}{2}}}{\displaystyle {\int }_{r_{i-\frac{1}{2}}}^{r_{i+\frac{1}{2}}}\rho \left(r\right)2\pi r^2\sin \theta drd\theta .}}\hfill \end{array}$$

Notice that the 2π factor appears in the equation as a result of the integral over the coordinate ϕ. That can be simplified by dividing the equation by 2π. Integrating this equation, we arrive at

(12) $$\begin{array}{l}\underset{i}{\underbrace{{{\displaystyle {\int }_{{\theta }_{j-\frac{1}{2}}}^{{\theta }_{j+\frac{1}{2}}}\left(r^2\frac{\partial }{\partial r}V\left(r,\theta \right)\right)}|}_{r_{i+\frac{1}{2}}}-{r^2\frac{\partial }{\partial r}V\left(r,\theta \right)|}_{r_{i-\frac{1}{2}}})\sin \theta d\theta }}\hfill \\ +\underset{ii}{\underbrace{{{\displaystyle {\int }_{r_{i-\frac{1}{2}}}^{r_{i+\frac{1}{2}}}\left(\sin \theta \frac{\partial }{\partial \theta }V\left(r,\theta \right)\right)}|}_{{\theta }_{j+\frac{1}{2}}}-\sin \theta {\frac{\partial }{\partial r}V\left(r,\theta \right)|}_{{\theta }_{j-\frac{1}{2}}})dr}}\hfill \\ =\underset{iii}{\underbrace{\frac{\cos {\theta }_{j+\frac{1}{2}}-\cos {\theta }_{j-\frac{1}{2}}}{{\epsilon }_o}{\displaystyle {\int }_{r_{i-\frac{1}{2}}}^{r_{i+\frac{1}{2}}}\rho \left(r\right)r^{2}dr.}}}\hfill \end{array}$$

where each term is described below.

1. The derivatives of the potential in the points $r_{i+\frac{1}{2}}$ and $r_{i-\frac{1}{2}}$ are given by

   (13a) $${r^2\frac{\partial }{\partial r}V\left(r,\theta \right)|}_{r_{i+\frac{1}{2}}}\simeq r_{i+\frac{1}{2}}^2\frac{V\left(r_{i+1},\theta \right)-V\left(r_i,\theta \right)}{\Delta r_i},$$

   and

   (13b) $${r^2\frac{\partial }{\partial r}V\left(r,\theta \right)|}_{r_{i-\frac{1}{2}}}\simeq r_{i-\frac{1}{2}}^2\frac{V\left(r_i,\theta \right)-V\left(r_{i-1},\theta \right)}{\Delta r_{i-1}}.$$

   We can now integrate these results. This yields:

   (14) $$\begin{array}{l}-\left(\cos {\theta }_{j+\frac{1}{2}}-\cos {\theta }_{j-\frac{1}{2}}\right)\hfill \\ \times (\frac{r_{i+\frac{1}{2}}^2}{\Delta r_i}V\left(r_{i+1},{\theta }_j\right)-\left(\frac{r_{i+\frac{1}{2}}^2}{\Delta r_i}+\frac{r_{i-\frac{1}{2}}^2}{\Delta r_{i-1}}\right)\hfill \\ \times V\left(r_i,{\theta }_j\right)+\frac{r_{i-\frac{1}{2}}^2}{\Delta r_{i-1}}V\left(r_{i-1},{\theta }_j\right)).\hfill \end{array}$$

2. Similarly, the derivatives of the potential in ${\theta }_{j+\frac{1}{2}}$ and ${\theta }_{j-\frac{1}{2}}$ are

   (15a) $${\sin \theta \frac{\partial }{\partial \theta }V\left(r,\theta \right)|}_{{\theta }_{j+\frac{1}{2}}}\simeq \sin {\theta }_{j+\frac{1}{2}}\frac{V\left(r,{\theta }_{j+1}\right)-V\left(r,{\theta }_j\right)}{\Delta {\theta }_j},$$

   and

   (15b) $${\sin \theta \frac{\partial }{\partial \theta }V\left(r,\theta \right)|}_{{\theta }_{j-\frac{1}{2}}}\simeq \sin {\theta }_{j-\frac{1}{2}}\frac{V\left(r,{\theta }_j\right)-V\left(r,{\theta }_{j-1}\right)}{\Delta {\theta }_{j-1}},$$

   integrating these terms

   (16) $$\begin{array}{l}\left(r_{i+\frac{1}{2}}-r_{i-\frac{1}{2}}\right)(\frac{\sin {\theta }_{j+\frac{1}{2}}}{\Delta {\theta }_j}V\left(r_{i,}{\theta }_{j+1}\right)\hfill \\ -\left(\frac{\sin {\theta }_{j+\frac{1}{2}}}{\Delta {\theta }_j}+\frac{\sin {\theta }_{j-\frac{1}{2}}}{\Delta {\theta }_{j-1}}\right)V\left(r_i,{\theta }_j\right)\hfill \\ +\frac{\sin {\theta }_{j-\frac{1}{2}}}{\Delta {\theta }_{j-1}}V\left(r_i,{\theta }_{j-1}\right)).\hfill \end{array}$$

3. The integral of the charge density involves the regions i − 1 and i. So, we split the integral to simplify the calculation,

   (17) $$\begin{array}{l}\frac{\cos {\theta }_{j+\frac{1}{2}}-\cos {\theta }_{j-\frac{1}{2}}}{{\epsilon }_o}\hfill \\ \times \left({\displaystyle {\int }_{r_{i-\frac{1}{2}}}^{r_i^{-}}\rho \left(r\right)r^{2}dr}+{\displaystyle {\int }_{r_i^{+}}^{r_{i+\frac{1}{2}}}\rho \left(r\right)r^{2}dr}\right).\hfill \end{array}$$

   The charge density depends on the region, so the integral of $\left[r_{i-\frac{1}{2}},r_i^{-}\right]$ is related to the region i − 1 and the integral of $\left[r_i^{+},r_{i+\frac{1}{2}}\right]$ is related to the region i, resulting in

   (18) $$\begin{array}{l}\frac{\cos {\theta }_{j+\frac{1}{2}}-\cos {\theta }_{j-\frac{1}{2}}}{3{\epsilon }_o}\hfill \\ \times \left({\rho }_{o,i-1}\left(r_i^3-r_{i-\frac{1}{2}}^3\right)+{\rho }_{o,i}\left(r_{i+\frac{1}{2}}^3-r_i^3\right)\right),\hfill \end{array}$$

   where ρ_o,i−1 and ρ_o,i are, respectively, the charge densities in the regions i − 1 and i and the radius $r_i=r_i^{+}=r_i^{-}$.

The reason for the use of integration techniques at the interface between different regions is that it already includes the continuity condition and derivative of the function at the interface. Combining the results obtained above in (i), (ii) and (iii):

(19) $$\begin{array}{l}-\left(\cos {\theta }_{j+\frac{1}{2}}-\cos {\theta }_{j-\frac{1}{2}}\right)\frac{r_{i-\frac{1}{2}}^2}{\Delta r_{i-1}}V\left(r_{i-1},{\theta }_j\right)\hfill \\ -\left(\cos {\theta }_{j+\frac{1}{2}}-\cos {\theta }_{j-\frac{1}{2}}\right)\frac{r_{i+\frac{1}{2}}^2}{\Delta r_{i-1}}V\left(r_{i+1},{\theta }_j\right)\hfill \\ +(\left(\cos {\theta }_{j+\frac{1}{2}}-\cos {\theta }_{j-\frac{1}{2}}\right)\left(\frac{r_{i+\frac{1}{2}}^2}{\Delta r_i}+\frac{r_{i-\frac{1}{2}}^2}{\Delta r_{i-1}}\right)\hfill \\ -\left(r_{i+\frac{1}{2}}-r_{i-\frac{1}{2}}\right)\left(\frac{\sin {\theta }_{j+\frac{1}{2}}}{\Delta {\theta }_j}+\frac{\sin {\theta }_{j-\frac{1}{2}}}{\Delta {\theta }_{j-1}}\right))V\left(r_i,{\theta }_j\right)\hfill \\ +\left(r_{i+\frac{1}{2}}-r_{i-\frac{1}{2}}\right)\frac{\sin {\theta }_{j-\frac{1}{2}}}{\Delta {\theta }_{j-1}}V\left(r_i,{\theta }_{j-1}\right)\hfill \\ +\left(r_{i+\frac{1}{2}}-r_{i-\frac{1}{2}}\right)\frac{\sin {\theta }_{j+\frac{1}{2}}}{\Delta {\theta }_{j-1}}V\left(r_i,{\theta }_{j+1}\right)\hfill \\ =\frac{\cos {\theta }_{j+\frac{1}{2}}\cos {\theta }_{j-\frac{1}{2}}}{3{\epsilon }_o}\left({\rho }_{o,i-1}\left(r_i^3-r_{i-\frac{1}{2}}^3\right)\right)\hfill \\ +{\rho }_{o,i}\left(r_{i+\frac{1}{2}}^3-r_i^3\right))\hfill \end{array}$$

which can be rewritten as follows,

(20) $$\begin{array}{l}{a}_{i-1,j}{V}_{i-1,j}+a_{i+1,j}{V}_{i+1,j}+a_{i,j}{V}_{i,j}+a_{i,j-1}{V}_{i,j-1}\hfill \\ +a_{i,j+1}{V}_{i,j+1}=s_{i,j},\hfill \end{array}$$

where 0 ≤ i ≤ I and 0 ≤ j ≤ J. I and J is the total number of regions. It is important to stress that we have I, J regions and I + 1, J + 1 points calculation. To simplify the notation, we have defined V_i,j ≡ V(r_i, θ_j). The mesh will be defined according to the size of the region of calculation, such that, Δr = R_I/I e Δθ = θ_J/J. The coefficients of Eq. (20) and the source term will be defined as follows,

(21) $$a_{i-1,j}\equiv -\left(\cos {\theta }_{j+\frac{1}{2}}-\cos {\theta }_{j-\frac{1}{2}}\right)\frac{r_{i-\frac{1}{2}}^2}{\Delta r_{i-1}},$$

(22) $$a_{i+1,j}\equiv -\left(\cos {\theta }_{j+\frac{1}{2}}-\cos {\theta }_{j-\frac{1}{2}}\right)\frac{r_{i+\frac{1}{2}}^2}{\Delta r_i}.$$

(23) $$\begin{array}{ll}{a}_{i,j}\equiv \hfill & \left(\cos {\theta }_{j+\frac{1}{2}}-\cos {\theta }_{j-\frac{1}{2}}\right)\left(\frac{r_{i+\frac{1}{2}}^2}{\Delta r_i}+\frac{r_{i-\frac{1}{2}}^2}{\Delta r_{i-1}}\right)\hfill \\ \hfill & -\left(r_{i+\frac{1}{2}}-r_{i-\frac{1}{2}}\right)\left(\frac{\sin {\theta }_{j+\frac{1}{2}}}{\Delta {\theta }_j}+\frac{\sin {\theta }_{j-\frac{1}{2}}}{\Delta {\theta }_{j-1}}\right).\hfill \end{array}$$

(24) $$a_{i,j-1}\equiv \left(r_{i+\frac{1}{2}}-r_{i-\frac{1}{2}}\right)\frac{\sin {\theta }_{j-\frac{1}{2}}}{\Delta {\theta }_{j-1}}.$$

(25) $$a_{i,j+1}\equiv \left(r_{i+\frac{1}{2}}-r_{i-\frac{1}{2}}\right)\frac{\sin {\theta }_{j+\frac{1}{2}}}{\Delta {\theta }_j}.$$

(26) $$\begin{array}{ll}{s}_{i,j}\equiv \hfill & \frac{\cos {\theta }_{j+\frac{1}{2}}-\cos {\theta }_{j-\frac{1}{2}}}{3{\epsilon }_o}\hfill \\ \hfill & \times \left({\rho }_{o,i-1}\left(r_i^3-r_{i-\frac{1}{2}}^3\right)+{\rho }_{o,i}\left(r_{i+\frac{1}{2}}^3-r_i^3\right)\right).\hfill \end{array}$$

In general, Eq. (20) is the numerical solution to Eq. (1) using the finite difference method, whose precision is directly related to the dimensions of the regions of calculation Δr_i = r_i+1 − r_i and Δθ_j = θ_j+1 − θ_j. The method proposed for numericalling solve this equation is known as Gauss-Seidel method [⁶] which is an iterative method developed to solve systems of linear equations. There are other methods of solution, such as the Jacobi method [⁶], that follows the same criteria of convergence as the Gauss-Seidel method, and non-iterative matrix methods, e.g., LU [⁶] factorization. The Gauss-Seidel method is generally used in problems involving diagonally dominant matrices, which is a necessary condition to ensure convergence.

The iteration of the Gauss-Seidel method to pentadiagonal matrices, as represented by Eq. (20), becomes

(27) $$\begin{array}{ll}{V}_{i,j}^t=\hfill & \frac{1}{a_{i,j}}(s_{i,j}-a_{i-1,j}{V}_{i-1,j}^t-a_{i,j-1}{V}_{i,j-1}^t\hfill \\ \hfill & -a_{i+1,j}{V}_{i+1,j}^{t-1}-a_{i,j+1}{V}_{i,j+1}^{t-1}),\hfill \end{array}$$

where t − 1 and t are, respectively, the previous and the current iteration. Notice that, in moving forward in the regions i, j, values of the potential previously calculated have been used, left $V_{i-1,j}^t$ and behind $V_{i,j-1}^t$ for previous regions, associated with the lower triangular matrix for the current iteration; on the other hand, values of the potential in the regions which have not yet been computed, right $V_{i,+1,j}^{t-1}$ and in front $V_{i,j+1}^{t-1}$, coming from the data initialization or previous iterations are associated with the upper triangular matrix. (For a better understanding of our notations left $V_{i-1,j}^t$, behind $V_{i,j-1}^t$, right $V_{i+1,j}^{t-1}$ and in front $V_{i,j+1}^{t-1}$, we refer the reader to the diagram depicted in Figure 1.) This is the great advantage of the Gauss-Seidel method when compared to the Jacobi method. The iterative process is over when one or more convergence criteria are reached. These criteria are obtained from the percentage relative errors of the potential between the current and previous iterations.

Once the discretized equation for the potential is known, the boundary conditions for this problem can be imposed.

2.3. Boundary condition

There are two possibilities to express the boundary conditions. The first one is to assume that the potential is zero on the boundary; the second one is to use a symmetry condition to simplify the problem. In this case, due to the azimuthal symmetry, the electric potential will be calculated in regions r ≥ 0. In the literature, a more general condition is known as Albedo boundary condition [⁷], given by

(28) $${\left[\alpha V+\beta \frac{d}{du}V\right]}_{u=u_n}=0,$$

where u = [r, θ] and u_n are the boundary regions, whereas u₀ and u_N are, respectively, left and right boundaries. The Eq. (28) shows that, if α = 0, the derivative of the potential is zero at the boundary point, representing the symmetry condition or singularity at the origin. Now, β = 0 implies that the potential is zero at a boundary point.

• Vanishing potential We are assuming that the potential is zero on the right boundary, i.e., V(r_I, θ_j) = 0. This is the simplest case, since that the potential is calculated in i = 0, …, I − 1 e j = 0, …, J.

• Condition of symmetry and singularity The symmetry conditions [⁸] are applied to the left boundary for the angle θ = 0° and right boundary for the angles θ = 90°, 180°. In this case, we impose that,

  (29) $$\begin{array}{ccc}{\frac{\partial }{\partial \theta }V\left(r,\theta \right)|}_{\theta ={\theta }_0}=0,& \forall & r\ge 0.\end{array}$$

  The singularity [³] is applied at the origin, r = 0. In this case, it is assumed that the derivative of the potential is zero, i.e.,

  (30) $$\begin{array}{ccc}{\frac{\partial }{\partial r}V\left(r,\theta \right)|}_{r=r_0}=0,& \forall & 0^{\circ }\le \theta \le {180}^{\circ }.\end{array}$$

  Eq. (19), in its simplified form given in Eq. (20), corresponds to the numerical solution of Eq. (1) for regions distant from the boundary. Now, to treat the condition of singularity and implement the symmetry, a new set of equations must be obtained to represent these regions we now wish to include. To meet our target, we believe it is more pedagogical to split the procedure in six steps, which altogether describe how we must work to include these new regions. The mentioned six steps are cast below:

  • The potential V(0, 0) is determined by integrating Eq. (1) over the interval $\left[r_0,r_{0+\frac{1}{2}}\right]$ and $\left[{\theta }_0,{\theta }_{0+\frac{1}{2}}\right]$. Thus, we solve the problem of the singularity at the origin, as well as the left symmetry for a null angle.

  • The potential V(0, θ_j) is determined by integrating Eq. (1) over the interval $\left[r_0,r_{0+\frac{1}{2}}\right]$ and $\left[{\theta }_{j-\frac{1}{2}},{\theta }_{j+\frac{1}{2}}\right]$. In this case, we treat the singularity at the origin.

  • The potential V(0, θ_J) is determined by integrating Eq. (1) over the interval $\left[r_0,r_{0+\frac{1}{2}}\right]$ and $\left[{\theta }_{J-\frac{1}{2}},{\theta }_J\right]$. As shown above, we solve the problem of the singularity and right symmetry for the angles 90° or 180°.

  • (iv)

    The potential V(r_i, 0) is determined by integrating Eq. (1) over the interval $\left[r_{i-\frac{1}{2}},r_{i+\frac{1}{2}}\right]$ and $\left[{\theta }_0,{\theta }_{0+\frac{1}{2}}\right]$. Thus, we treat the problem of left symmetry for a null angle.

  • (v)

    The potential V(r_i, θ_j) is given by Eq. (19).

  • (vi)

    The potential V(r_i, θ_J) is determined by integrating Eq. (1) over the interval $\left[r_{i-\frac{1}{2}},r_{i+\frac{1}{2}}\right]$ and $\left[{\theta }_{J-\frac{1}{2}},{\theta }_J\right]$. The right symmetry for the angles 90° or 180° has been treated in a similar way.

3. Results

The comparison between the results obtained by the numerical and analytic solutions to the Poisson’s equation shall be now presented. In these calculations, it has been assumed Δr = 10⁻¹⁶ m and Δθ = 0.5°. The potential has been calculated in the domain of the electron radius. The electric charge density has been obtained by means of the following relation,

(31) $${\rho }_o=\frac{Q}{\frac{4\pi R_e^3}{3}},$$

where Q = 1.6 × 10⁻¹⁹C is the electron charge and R_e ≃ 3.0 × 10⁻¹⁵ m the electron radius. The convergence criteria of ∈ = 10⁻⁶ for the calculation of the electric potential has been used.

The numerical method works with finite boundary conditions. Eq. (7) was obtained by using the condition that the potential outside the charged sphere is zero at infinity. So, to compare the numerical results with this analytic solution, we must assume a finite physical boundary for analytic solution and recalculate the coefficients of the solution. Notice that the solutions for the electric potential inside and outside the charged sphere, given by Eqs. (4) and (5), remain unchanged. The result obtained with the numerical solution for θ = 180° is cast in Figure 2.

The comparison between the results obtained by the analytic and numerical solutions for the electric potential is depicted in Figure 3. From the solution to the electric potential, the numerical solution to the electric field is readily attained. This has been calculated by approximating the derivative of Eq. (8) by finite differences, such that

Figure 2
Electric potential in spherical coordinates.

Figure 3
Comparison of results for the electric potential.

(32) $$\begin{array}{ll}\text{E}\left(r_i,{\theta }_j\right)=\hfill & -\frac{V\left(r_i,{\theta }_j\right)-V\left(r_{i-1},{\theta }_j\right)}{\Delta r_{i-1}}\widehat{\text{r}}\hfill \\ \hfill & -\frac{1}{r_i}\frac{V\left(r_i,{\theta }_j\right)-V\left(r_{i,}{\theta }_{j-1}\right)}{\Delta {\theta }_{j-1}}\widehat{\theta }.\hfill \end{array}$$

Once again, a comparison between the analytic and numerical solutions to the radial component was established in terms of the electric field (Figure 4).

Figures 3 and 4 show a qualitative comparison between the solutions. To quantify such a comparison, the percentage relative errors have been calculated using the following equation,

(33) $$Error\left(\%\right)=\left|\frac{R_a-R_n}{R_a}\right|\times 100\%,$$

where R_a and R_n represent the analytic and numerical results, respectively. For a better assessment, the resulting potential and radial component of the electric field and their relative errors are summarized in Table 1. The maximum percentage relative error obtained in the calculation of the electric potential was 3.9%.

Figure 4
Comparison of results for the electric field.

4. Concluding Comments

In this contribution, our major effort consisted in setting up a formulation to solve numerically Poisson’s equation in spherical coordinates; continuity conditions on the function and its derivative at the interface between non-homogeneous media have been imposed. By following this procedure, the electrostatic potential and field have been determined in the inner and external regions, and also compared with the analytic solution of the problem to give confidence that the method is reliable. For this, a program has been developed that allows to define spatial meshes with different sizes in each region. The accuracy of the method is directly related to the sizes of the meshes of the regions appearing in the calculations. We may therefore refine the spatial meshes to obtain more accurate results, since the angular mesh expands the radial solution at different angles within the angular domain, not affecting the radial solution of the problem. Fortran 90 has been adopted to develop this code.

By integrating numerically the Poisson’s equation, the singularity problem and symmetry have been treated and the electric potential has been written down in all radial and angular regions. Both qualitative and quantitative comparisons between the numerical and analytic solutions show a very satisfactory result of the numerical method, with relative errors less than 1% for distances of the order of 1 × 10⁻¹⁴ m from the electric charge and always less than 4% for larger distances. We conclude by stating that this method is very helpful for the determination of the electrostatic potential in different non-homogeneous regions. Therefore, with the validation we have ascertained in the previous Sections, the method can be reliably applied to solve problems of a greater complexity.

Table 1
Results for the electric potential and electric field.

r (fm)	Va(r)(V)	Vn(r)(V)	Error(%)	Ea(r)(V/m)	En(r)(V/m)	Error(%)
0	6.962 × 105	6.952 × 105	0.1436	0	0	0
1.50	6.363 × 105	6.353 × 105	0.1571	7.989 × 1019	7.722 × 1019	3.3420
10.0	1.210 × 105	1.201 × 105	0.7438	1.438 × 1019	1.452 × 1019	0.9735
30.0	2.511 × 104	2.451 × 104	2.3894	1.598 × 1018	1.585 × 1018	0.8135
50.0	5.935 × 103	5.719 × 103	3.6394	5.752 × 1017	5.578 × 1017	3.0250
60.0	1.141 × 103	1.097 × 103	3.8562	3.994 × 1017	3.847 × 1017	3.6805

• fm = 10⁻¹⁵ m.
• _a(r) is the analytic potential.
• _n(r) is the numerical potential for θ = 0°.

As a natural follow-up of the method and the results presented in this contribution, it would be advisable to consider the case of a non-spherically symmetric charge density, by contemplating the situation of a θ-dependent charge distribution. We shall endeavour this study and we shall be presenting our results in a forthcoming paper.

References

• ^[1]
  C.S. Tan, Journal of Computational Physics 59, 81 (1985).
• ^[2]
  J.P. Boyd and C. Zhou, Journal of Computational Physics 228, 4702 (2009).
• ^[3]
  O.D. Momoh, M.N.O. Sadiku and C.M. Akujuobi, in: Proceedings of the PIERS Conference, (PIERS, Cambridge, 2010).
• ^[4]
  S. Nakamura, Computational Methods in Engineering and Science (John Wiley & Sons, New York, 1977).
• ^[5]
  D.J. Griffiths, Introduction to Electrodynamics (Prentice Hall, New Jersey, 1999).
• ^[6]
  R.L. Burden and J.D. Faires, Numerical Analysis (Cengage Learning, Boston, 2011).
• ^[7]
  W.M. Stacey, Nuclear Reactor Physics (Wiley-VCH, Weinheim, 2007).
• ^[8]
  O.D. Momoh and M.N.O. Sadiku, Int. J. of Modelling and Simulations 2, 196 (2009).

Datas de Publicação

• Publicação nesta coleção
  19 Maio 2021
• Data do Fascículo
  2021

Histórico

• Recebido
  16 Jan 2021
• Revisado
  09 Abr 2021
• Aceito
  19 Abr 2021