An iterative Monte Carlo method to solve nonlinear second-order differential equations
Martín Chávez-Páez
Instituto de Física, Universidad Autónoma de San Luis Potosí, Álvaro Obregón 64, 78000 San Luis Potosí, S.L.P., México
Enrique González-Tovar
Instituto de Física, Universidad Autónoma de San Luis Potosí, Álvaro Obregón 64, 78000 San Luis Potosí, S.L.P., México
Corresponding Author: enrique.gonzalez@uaslp.mx
Guillermo Iván Guerrero-García
Facultad de Ciencias, Universidad Autónoma de San Luis Potosí, Av. Chapultepec 1570, Privadas del Pedregal, 78295 San Luis Potosí, S.L.P., México
Abstract
The Monte Carlo method is a thriving and mathematically beautiful numerical technique used extensively, nowadays, to deal with many demanding problems in diverse fields. Here, we present an iterative Monte Carlo algorithm to work out very general nonlinear second-order differential equations, with Dirichlet boundary conditions. An example of its usage is, also, reported.
1 Introduction
Since the conception of the Monte Carlo method (MCM) in its modern form [1], this tantalising mathematical idea has become a powerful and successful route to solve extremely difficult and even “intractable” [2] problems in many areas of knowledge [3]. In particular, the usage of the Monte Carlo probabilistic approach for the numerical solution of differential equations has a long history [4] and represents one of its most useful applications in science and technology. In this regard, the amazing relation discovered between Brownian motion and potential theory [5] has prompted the employment of the random walk-based MCM to solve the very fundamental Poisson equation [3, 6], i.e., . Despite the importance of this subject, to the best of our knowledge, there exists no general Monte Carlo algorithm to solve the Poisson equation when also depends on and/or on its derivatives, apart from an explicit dependence on . In this context, below, we introduce and exemplify an iterative Monte Carlo method that allows to treat numerically ordinary differential equations of the form , with Dirichlet boundary conditions. Notably, this reported probabilistic procedure can be readily expanded to deal with partial differential equations.
2 The method
A second-order ordinary differential equation (ODE) for , with Dirichlet boundary conditions, is generally stated as:
the Iterative Monte Carlo Method (IMCM), discussed in the following, can be utilised to solve Eq. (2), for and . As a preliminary, let us, then, consider a simpler case of such equation, i.e.,
(3)
(notice that the RHS term of the previous equation only contains an explicit dependence on the variable ). The usual random walk-based Monte Carlo algorithm to solve Eq. (3) comprises the next steps [6]:
STEP 1: Given a uniform mesh , for each free point , such that , generate a large number, , of random walks, which start at and end when hitting an absorbing boundary site, either or .
STEP 2: If the -th walk arrives at the boundary after steps and has visited the sequence of locations , calculate the Monte Carlo estimator for from
(4)
where and being the value, or , of the function at the particular absorbing boundary point, or , reached by the -th random walk initiated at . Observe that the approximate values of for all the points in the mesh can be obtained for any arbitrary sequence of the indices in the set , and, besides, that the computer implementation of the procedure can be easily parallelised.
Now, regarding the more general ODE stated by Eq. (2), in theory, the typical random walk methodology reviewed above can not be applied to work out such ODE, since the RHS of Eq. (2) depends, precisely, on the unknown functions and . However, and given that, at the end, these two sought functions depend on , the first version of our Iterative Monte Carlo Method (IMCM) proposes the use of an initial guess function, , with a well-known and simple dependence on in the interval and that fulfills the associated boundary conditions and ; this trial function will allow us to evaluate the last term in the RHS of Eq. (4) and, consequently, to set up the ensuing practicable estimator for the value of the next guess function at :
(5)
for , with being the numerical derivative of at each node (which can be calculated through any finite difference formula). After determining the values of for all the free nodes, the corresponding numerical derivative has to be got to proceed with successive approximations of , obtained via the iteration prescription
(6)
complemented with a finite difference approximation of its derivative. In a similar way to the scheme associated to Eq. (4), the prior iteration expression can be used to get for all the free points following any order of the mesh indices. Additionally, this first version of the IMCM, given by Eq. (6), can be likewise parallelised.
On the other hand, based on the idea of successive displacement we can try now to improve our IMCM as follows:
STEP 1: Advance an initial guess function, , in , that satisfies the boundary conditions of the problem and, then, initialise , for .
STEP 2: For each free node, perform absorbing random walks starting at to obtain and update immediately employing the formula
(7)
STEP 3: Right after calculating in STEP 2, differentiate it numerically to produce .
STEP 4: Repeat entirely the previous steps 2 and 3 to complete a cycle of iterations.
This last stochastic and successive displacement process represents a refined, and seemingly hastened, version of the IMCM introduced here (the pseudocode of the corresponding algorithm can be consulted in Fig. 1).
Figure 1: Pseudocode for the Iterative Monte Carlo Method (successive displacement version).
3 An example of application
Let us consider the following nonlinear second-order ODE of the type [7]:
which it will be utilised below to asses the performance of our approximate numerical IMCM results.
Coding the enhanced IMCM procedure described at the end of Section 2 (see Fig. 1), such that
(10)
we obtained the data included in Fig. 2. Therein, we show the exact solution given by Eq. (9) (red open symbols), the linear initial guess (green interrupted line), the IMCM results after a first iteration with Eq. (7) (blue continuous line), and the IMCM prediction after ten iterations with Eq. (7) (black continuous line). In each Monte Carlo iteration we have employed a fixed grid in , with equally-sized partitions, combined with random walks, to estimate the value of at each non-boundary node of the grid via Eq. (7). The error in the -th iteration, is calculated from
(11)
with being the -th IMCM approximation. The corresponding errors for the IMCM results for iterations 1 and 10, plotted in Fig. 2, are and , respectively.
Figure 2: Solution to the ODE specified in Eq. (8): exact solution given by Eq. (9) (red open symbols); linear initial guess (green interrupted line); IMCM numerical approximate results after one iteration with Eq. (7) (blue continuous line); IMCM numerical approximate results after ten iterations with Eq. (7) (black continuous line).
4 Conclusions
In this paper, we have proposed an Iterative Monte Carlo Method to solve second-order ODEs of the general form , with Dirichlet boundary conditions. Our algorithm is based on the classic random walk approach, but it is enriched with an iterative process starting from a trial or guess function, which allows the evaluation of the RHS of the ODE. To illustrate its performance, we have also included an example of application. The presented IMCM can be extended straightforwardly to second-order partial differential equations.
References
[1] N. G. Cooper, N. G., Los Alamos Science Number 15: Special Issue on Stanislaw Ulam 1909-1984 (1987). doi:10.2172/1054744
[2] J. F. Traub and H. Woźniakowski, Breaking Intractability, Scientific American270, Issue 1, 102 (1994). doi:10.1038/scientificamerican0194-102
[3] J. M. Hammersley and D. C. Handscomb, Monte Carlo Methods (Methuen & Co. Ltd., 1964).
[4] R. Courant, K. Friedrichs, and H. Lewy, Über die partiellen Differenzengleichungen der mathematischen Physik, Math. Ann.100, 32 (1928). doi.org/10.1007/BF01448839 (English translation in IBM J. Res. Develop.11, 215 (1967). doi.org/10.1147/rd.112.0215).
[5] R. Hersh and R. J. Griego, Brownian Motion and Potential Theory, Scientific American220, Issue 3, 66 (1969). doi:10.1038/scientificamerican0369-66
[6] M. N. O. Sadiku, R. C. Garcia, S. M. Musa, and S. R. Nelatury, Markov Chain Monte Carlo Solution of Poisson’s Equation, International Journal on Recent and Innovation Trends in Computing and Communication3, 106 (2015). doi:10.17762/ijritcc.v3i1.3770
[7] P. Yin, Y. Huang, and H. Liu, An Iterative Discontinuous Galerkin Method for Solving the Nonlinear Poisson Boltzmann Equation, Communications in Computational Physics16, 491 (2014). doi:10.4208/cicp.270713.280214a