Markov Chain Monte Carlo Solution of Laplace’s Equation in Axisymmetric Homogeneous Domain

Adebowale E. Shadare, Matthew N. O. Sadiku, Sarhan M. Musa

Department of Electrical/Computer Engineering, Prairie View A&M University, Prairie View, TX, USA.

DOI: 10.4236/ojmsi.2019.74012                

Abstract

With increasing complexity of today’s electromagnetic problems, the need and opportunity to reduce domain sizes, memory requirement, computational time and possibility of errors abound for symmetric domains. With several competing computational methods in recent times, methods with little or no iterations are generally preferred as they tend to consume less computer memory resources and time. This paper presents the application of simple and efficient Markov Chain Monte Carlo (MCMC) method to the Laplace’s equation in axisymmetric homogeneous domains. Two cases of axisymmetric homogeneous problems are considered. Simulation results for analytical, finite difference and MCMC solutions are reported. The results obtained from the MCMC method agree with analytical and finite difference solutions. However, the MCMC method has the advantage that its implementation is simple and fast.

Keywords

Laplace’s Equation, Axisymmetric Problem, Inhomogeneous Dirichlet Boundary Conditions, Markov Chain Monte Carlo (MCMC)

1. Introduction

Most real-world EM problems are difficult to solve using analytical methods and in most cases, analytical solutions are outright intractable [1]. Today, with increasing advancement in computer technologies and with increasing system complexities, the need for continuous development of computational techniques to address contemporary EM problems is as ever, critical. EM problems in nature are essentially three-dimensional (3D) requiring, in most cases, enormous time and computer memory to solve. However, approximations in 1D and 2D are often a smart way of solving 3D problems provided that approximations can be made without loss of accuracy. In axisymmetric systems, the cylindrical coordinate, is essentially independent of, reducing the complexity of the problem to two-dimensional problems in plane of the three-dimensional domain [2] - [10]. The resulting two-dimensional approximation can be treated like a Cartesian coordinate problem with either z or constant interface.

Several methods such as the Method of Lines [11], the Finite Element method [12] [13], Finite difference method [14] and the Boundary Integral Equation methods [15] [16] have all been applied in the modeling and analysis of axisymmetric problems. Since the connection between Brownian motion and the potential theory was established [17] and the application of probabilistic potential theory to electrical engineering related problems [18], several Monte Carlo techniques such as fixed random walk, floating random walk and Exodus method [9] [19] - [26] have evolved dramatically. However, conventional Monte Carlo techniques are limited in application because they are unsuitable for whole-field computation. They only allow single-point calculations. The shrinking boundary and the inscribed figure methods later proposed for whole-field calculations are not significantly superior to the classical Monte Carlo methods [27] [28]. To address this gap, Markov Chains for whole-field computations was proposed by Andrey Markov [29] [30]. The applications of MCMC to rectangular and axisymmetric problems are presented in [29] [31].

However, to the best of the authors’ knowledge, the solutions of Laplace’s equation in axisymmetric homogeneous domains with MCMC with inhomogeneous Dirichlet boundary condition are yet to be reported in the literature. Thus, this paper presents the MCMC solution of Laplace’s equation in axisymmetric homogeneous region.

2. Axisymmetric Problem Formulation

When it is necessary and convenient, electromagnetic problems in cylindrical coordinates may be approximated to axisymmetric solution region. Suppose a cylindrical coordinate system is as shown in Figure 1. It is immediately clear that the corresponding axisymmetric approximation in meshed region is as shown in Figure 2. The voltages, and imposed in Figure 1 satisfy the Dirichlet boundary condition, justifying axisymmetric approximation shown in Figure 2.

The Laplace’s equation in the axisymmetric region R depicted in Figure 2 is given as

(1)

The corresponding finite difference equivalence of Equation (1) for region using square grid is given as

(2)

Figure 1. Cylindrical geometry.

Figure 2. Axisymmetric Solution Region.

The transition probabilities in the Equation (2) are given as [14]

;;

Similarly, the finite difference equivalence of Equation (1) for region is given as [14]

(3)

The corresponding transition probabilities in the Equation (3) are given as [14]

;;

3. Absorbing Markov Chain

A Markov chain is a mathematical model that represents a sequence of random variables, such that the probability distribution of depends on the probability distribution of [14] [32] [33]. The process has a memoryless property—remembering only the most recent past. This paper considers discrete-state, discrete-time Markov chains where the Markov chain represents the random walk and the states represent the grid nodes. The transition probability is the probability that a randomly walking particle starting at node i will move to node j. The is expressed as

(4)

The transition probability P is defined as

; (5)

The sum of each row elements in P matrix is 1. This shows that the matrix P is stochastic as described in Equation (5).

The size of the transition matrix P is, where

(6)

From the Equation (6), are the free (non-absorbing) nodes which represent the nodes in the solution region excluding those on the boundary., on the other hand, are fixed (absorbing) nodes which are nodes on the boundary.

The transition matrix P in which absorbing nodes and the non-absorbing nodes are numbered first and last respectively is given by [14]

(7)

where

R is a matrix which describes the probabilities of moving from free nodes to fixed ones;

Q is a matrix which describes the probabilities of moving from one fixed node to another;

I is a identity matrix which describes the transitions between fixed nodes (and);

0 is a null matrix which shows that no transitions exist from fixed to free nodes.

To solve Laplace’s equation in the region R, the elements of matrix Q for nodes in the region are obtained from the Equation (2) as

(8)

Equation (8) describes the probabilities of moving from one node to another in the region.

Similarly, the elements of matrix Q for nodes at the line of symmetry, , are obtained from the Equation (3) as

(9)

The same Equation (8) and Equation (9) apply to except that j is a fixed node.

For any absorbing Markov chains, the fundamental matrix, describes the average number of times the randomly walking particles originating from node i passes through the node j before being absorbed. The fundamental matrix N with size is given as

(10)

The absorption probability matrix B which describes the probabilities that a randomly walking particle starting at free node i will be absorbed at a fixed node j is given as

(11)

The size of matrix B is. The matrix B is stochastic and it is given as

(12)

(13)

where with size and with size are free and fixed (prescribed) nodes potentials respectively.

Finally, the potential at any free node i is given as

(14)

where are the prescribed potentials,.

4. Simulation Results

In this section, simulation results are presented for the two cases of homogeneous axisymmetric problems considered. The inhomogeneous Dirichlet boundary condition with different levels of complexity was enforced for the two cases presented. At the line of symmetry, the Neumann boundary condition was imposed. Simulation results are reported for both cases.

Case I: Axisymmetric homogeneous Problem with Inhomogeneous Dirichlet Boundary Condition

Figure 3. Meshed Solution Region for Case I.

Suppose the axisymmetric cylinder given in Figure 3 has its ends at

and grounded respectively. The boundary condition imposed at is given as [14]:

(15)

In order to demonstrate the effectiveness of the MCMC method for axisymmetric homogeneous problems, the following simulation is carried out. The parameters used for the simulations are given in Table 1. Simulation is performed on MATLAB. With the length of the cylinder, and radius, given as dimensions of the axisymmetric region R, the solution region is discretized into grids using a step size of 0.025 m. This gives 3160 grid points excluding the nodes on the boundary. The inhomogeneous Dirichlet boundary condition imposed at is given in Equation (15). Other step sizes used are 0.05 m, 0.1 m and 0.2 m giving 780, 190 and 45 grid points respectively.

From the given parameters, the elements of matrix Q with the size 3160 × 3160 are formed based on the Equation (8) and Equation (9). Similarly, the elements of matrix R with size 3160 × 159 are formed from the Equation (8) and Equation (9) except that j is a fixed node. Then the identity matrix I with size 3160 × 3160 is determined.

Based on the foregoing, the fundamental matrix N with size 3160 × 3160 and the absorption matrix B with size are determined. Finally, the free node potentials are then determined from. The prescribed potential, has a size 159 × 1 and it is determined from the nodes on the boundary of the solution region. The simulation results from MCMC method are validated with the analytical solution given in [14]:

Table 1. Parameters for homogeneous axisymmetric problems.

(16)

where

is a modified Bessel function of the order zero.

A solution to the present problem is also obtained with the finite difference method and comparison was made with the MCMC and analytical solutions. Simulation results are presented in Figure 4. The results for the potential distributions along (middle of the grid nodes), , and are presented in Figures 4(a)-(c). At, the potential distribution at the line of symmetry is reported. Also, the potential distribution along as well as the surface and contour plots for all the grid nodes are reported in Figures 4(d)-(f). The MCMC simulations were repeated with step sizes of, and respectively and the results reported as shown in Figure 4. The smaller the step size, the closer the results are to the exact solution. In Table 2, the MCMC, analytical and finite difference solutions for some selected grid points are compared. As evident, the MCMC solutions agree well with the results obtained with the analytical solution and the finite difference method (FDM). The computation time for MCMC and FDM for this problem is respectively 2.1717 seconds and 1.0486 seconds. Hence, the computation time cannot be used as a basis for comparison.

Case II: Axisymmetric Homogeneous Problem with Inhomogeneous Dirichlet Boundary Condition.

In this section, another case of Laplace’s equation with inhomogeneous Dirichlet boundary conditions is presented. Suppose the cylinder given in Figure 5 has its ends grounded and given as 100 V respectively. The boundary at is described as in the Equation (15) where voltage is 100 V. The Neumann boundary condition is imposed at. Since the Laplace’s equation is a linear homogeneous equation, the problem in Figure 5 is simplified into Figure 5(a) and Figure 5(b). An analytical solution to the solution Region I is given as [34]:

Figure 4. Potential distribution along (a); (b) line of symmetry,; (c) ; (d); (e) surface plot; (f) contour plot for Laplace’s equation in Homogeneous Axisymmetric domain with Inhomogeneous Dirichlet boundary condition: Case I.

Table 2. Results comparison for analytical, FDM and MCMC for Laplace’s equation with inhomogeneous Dirichlet boundary condition: Case I.

(17)

where a and L are given in Table 1. are the roots of; and are Bessel functions of the first kind, order zero and one, respectively.

Similarly, recall that the analytical solution for region II depicted in Figure 5(b) is already provided in the Equation (16). Since the Laplace’s equation can be linearized, the analytical solution to the problem depicted in Figure 5 is given

(a)(b)(c)

Figure 5. (a) Meshed Solution Region for Case II; (b) Solution Region I; (c) Solution Region II.

as:

(18)

where

is a modified Bessel function of the order zero.

With the step size of 0.025 m, the MCMC results are reported in the Figure 6. The potential distributions along, (symmetry line) and are given in the Figures 6(a)-(c). Similarly, the potential distribution along, the surface plot and the contour plot for all the grid nodes are reported in the Figures 6(d)-(f) respectively. The results for selected grid points for MCMC and FDM using step size of 0.025 m agree perfectly with the analytical solution as shown in the Table 3. Computation time for MCMC and FDM for the problem is 1.8602 seconds and 0.8658 seconds respectively. The reduction in step size increases the accuracy of the solutions.

Figure 6. Potential distribution along (a); (b) line of symmetry,; (c) ; (d); (e) surface plot; (f) contour plot for Laplace’s equation in Axisymmetric Homogeneous Domain with Inhomogeneous Dirichlet Boundary Condition: Case II.

Table 3. Results comparison for analytical, FDM and MCMC for Laplace’s equation with inhomogeneous Dirichlet boundary condition: Case II.

5. Conclusion

This paper presents a comprehensive application of MCMC method to Laplace’s equation in homogeneous axisymmetric domains. Two broad cases of homogeneous axisymmetric problems were investigated. For Case I, the MCMC method was used to solve Laplace’s equation with inhomogeneous Dirichlet boundary condition in which the top and bottom boundaries were characterized by the same potential. For case II, the MCMC method was investigated with the Laplace’s equation using inhomogeneous Dirichlet boundary condition in which the top and bottom boundaries (prescribed potentials) were at different potentials. Also, Neumann boundary condition was imposed at the line of symmetry. Several plots were reported. The MCMC results were compared with the analytical and finite difference solutions. The proposed MCMC method agreed with the analytical and finite difference solutions for all reported cases. However, the difference in computation time between MCMC and FDM is in the order of seconds for the problems considered and thus cannot be used as a basis for comparison.

Nomenclature

MCMC: Markov Chain Monte Carlo;

1D: One dimension;

2D: Two dimension;

3D: Three dimension;

EM: Electromagnetics;

FDM: Finite difference method;

P: Transition probability;

N: Fundamental matrix;

B: Absorption probability matrix.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Sadiku, M.N.O. and Nelatury, S.R. (2016) Analytical Techniques in Electromagnetics. CRC Press, Boca Raton, FL.
[2]	Louai, F.Z., Said, N.N. and Drid, S. (2006) Numerical Analysis of Electromagnetic Axisymmetric Problems Using Element Free Galerkin Method. Journal of Electrical Engineering, 57, 99-104.
[3]	Hayati, A.N., Ahmadi, M.M. and Sadrnejad, S.A. (2012) Analysis of Axisymmetric Problems by Element-Free Galerkin Method. International Journal of Modeling and Optimization, 2, 712-717. https://doi.org/10.7763/IJMO.2012.V2.217
[4]	Soares, R.D., Moreira, F.J.S., Mesquita, R.C., Lowther, D.A. and Lima, N.Z. (2014) A Modified Meshless Local Petrov-Galerkin Applied to Electromagnetic Axisymmetric Problems. IEEE Transactions on Magnetics, 50, Article No. 7012604. https://doi.org/10.1109/TMAG.2013.2284472
[5]	Freeman, E.M. and Lowther, D.A. (1989) An Open Boundary Technique for Axisymmetric and Three Dimensional Magnetic and Electric Field Problems. IEEE Transactions on Magnetics, 25, 4135-4137. https://doi.org/10.1109/20.42546
[6]	Melissen, J.B.M. and Simkin, J. (1990) A New Coordinate Transform for the Finite Element Solution of Axisymmetric Problems in Magnetostatics. IEEE Transactions on Magnetics, 26, 391-394. https://doi.org/10.1109/20.106336
[7]	Boglietti, A., Chiampi, M., Chiarabaglio, D. and Tartaglia, M. (1990) Finite Element Approximation in Axisymmetrical Domains. IEEE Transactions on Magnetics, 26, 395-398. https://doi.org/10.1109/20.106337
[8]	Yoneta, A., Tsuchimoto, M. and Honma, T. (1990) An Analysis of Axisymmetric Modified Helmholtz Equation by Using Boundary Element Method. IEEE Transactions on Magnetics, 26, 1015-1018. https://doi.org/10.1109/20.106492
[9]	Sadiku, M.N.O. (1993) Monte Carlo Solution of Axisymmetric Potential Problems. IEEE Transactions on Industry Applications, 29, 1042-1046. https://doi.org/10.1109/28.259710
[10]	OH, M. (2010) Efficient Solution Techniques for Axisymmetric Problems. University of Florida, Gainesville, FL.
[11]	Sadiku, M.N.O. and Garcia, R.C. (2000) Method of Lines Solution of Axisymmetric Problems. Proceedings of IEEE SoutheastCon 2000, Nashville, TN, 9 April 2000, 527-530.
[12]	Gratkowski, S., Todaka, T., Enokizono, M. and Sikora, R. (2000) Asymptotic Boundary Conditions for the Finite Element Modeling of Axisymmetric Electrical Field Problems. IEEE Transactions on Magnetics, 36, 717-721. https://doi.org/10.1109/20.877549
[13]	Chen, Q., Konrad, A. and Baronijan, S. (1994) Asymptotic Boundary Conditions for Axisymmetric Finite Element Electrostatic Analysis. IEEE Transactions on Magnetics, 30, 4335-4337. https://doi.org/10.1109/20.334079
[14]	Sadiku, M.N.O. (2019) Computational Electromagnetics with MATLAB. 4th Edition, CRC Press, Boca Raton, FL. https://doi.org/10.1201/9781315151250
[15]	Ehrich, M., Kuhlmann, J. and Netzler, D. (1997) High Accuracy Integration of Boundary Integral Equations Describing Axisymmetric Field Problems. Proceedings of Asia Pacific Microwave Conference, Hong Kong, 2-5 December 1997, 462-464. https://doi.org/10.1109/APMC.1997.659423
[16]	Krahenbuhl, L. and Nicolas, A. (1983) Axisymmetric Formulation for Boundary Integral Equation Methods in Scalar Potential Problems. IEEE Transactions on Magnetics, 19, 2364-2366. https://doi.org/10.1109/TMAG.1983.1062872
[17]	Kakutani, S. (1944) Two Dimensional Brownian Motion and Harmonic Function. Proceedings of Imperial Academy (Tokyo), 20, 706-714. https://doi.org/10.3792/pia/1195572706
[18]	Bevensee, R.M. (1973) Probabilistic Potential Theory Applied to Electrical Engineering Problems. Proceedings of IEEE, 61, 423-427. https://doi.org/10.1109/PROC.1973.9056
[19]	Garcia, R.C. and Sadiku, M.N.O. (1996) Monte Carlo Fixed-Radius Floating Random Walk Solution for Potential Problems. Proceedings of SoutheastCon, Tampa, FL, 11-14 April 1996, 88-91. https://doi.org/10.1109/SECON.1996.510032
[20]	Sadiku, M.N.O. and Garcia, R.C. (1993) Monte Carlo Floating Random Walk Solution of Poisson’s Equation. Proceedings of SoutheastCon, Charlotte, NC, 4-7 April 1993, 4.
[21]	Gu, K. and Sadiku, M.N.O. (1995) A Triangular Mesh Random Walk Method for Dirichlet Problem. Journal of Franklin Institute, 332, 569-578. https://doi.org/10.1016/0016-0032(95)00081-X
[22]	Sadiku, M.N.O. and Gu, K. (1994) Floating Random Walk Method on Axisymmetric Potential Problems. International Symposium on Electromagnetic Compatibility, Japan, 1994, 659-662.
[23]	Sadiku, M.N.O., Ajose, S.O. and Fu, Z. (1994) Applying the Exodus Method to Solve Poisson’s Equation. IEEE Transactions on Microwave Theory and Techniques, 42, 661-666. https://doi.org/10.1109/22.285073
[24]	Sadiku, M.N.O. and Gu, K. (1996) A New Monte Carlo Method for Neumann Problems. Proceedings of SoutheastCon, Tampa, FL, 11-14 April 1996, 92-95. https://doi.org/10.1109/SECON.1996.510033
[25]	Godoy, E., Boccardo, V. and Durán, M. (2017) A Dirichlet-to-Neumann Finite Element Method for Axisymmetric Elastostatics in a Semi-Infinite Domain. Journal of Computational Physics, 328, 1-26. https://doi.org/10.1016/j.jcp.2016.09.066
[26]	Beasley, M.D.R., et al., (1979) Comparative Study of Three Methods for Computing Electric Fields. Proceedings of Institution of Electrical Engineers, 126, 126-134.
[27]	Zinsmeiter, G.E. and Sawyerr, J.A. (1974) A Method for Improving the Efficiency of Monte Carlo Calculation of Heat Conduction Problems. Transactions of the ASME, Journal of Heat Transfer, 96, 246-248. https://doi.org/10.1115/1.3450172
[28]	Zinsmeiter, G.E. and Pan, S.S. (1976) A Modification of the Monte Carlo Method. International Journal for Numerical Methods in Engineering, 10, 1057-1064. https://doi.org/10.1002/nme.1620100507
[29]	Sadiku, M.N.O. and Garcia, R.C. (1998) Whole Field Computation Using Monte Carlo Method. International Journal of Numerical Modelling: Electronic Networks, Devices and Fields, 10, 303-312. https://doi.org/10.1002/(SICI)1099-1204(199709/10)10:5<303::AID-JNM281>3.0.CO;2-R
[30]	Fusco, V.F. and Linden, P.A. (1988) A Markov Chain Approach for Static Field Analysis. Microwave and Optical Technology Letter, 1, 216-220. https://doi.org/10.1002/mop.4650010610
[31]	Garcia, R.C. and Sadiku, M.N.O. (1998) Neuro-Monte Carlo Solution of Electrostatic Problems. Journal of Franklin Institute, 335, 53-69. https://doi.org/10.1016/S0016-0032(96)00115-9
[32]	Sadiku, M.N.O. (1990) Monte Carlo Methods in an Introductory Electromagnetic Course. IEEE Transactions on Education, 33, 73-80. https://doi.org/10.1109/13.53630
[33]	Sadiku, M.N.O. (2009) Monte Carlo Methods for Electromagnetics. CRC Press, Boca Raton, FL. https://doi.org/10.1201/9781439800720
[34]	Gu, K. (1996) Monte Carlo Solution for Potential and Waveguide Problems. Ph.D. Dissertation, Temple University, Philadelphia, PA.