## An iterative Monte Carlo method to solve nonlinear second-order differential equations

Mart´ ın Ch´vez-P´ez , Enrique Gonz´lez-Tovar a a a a a,* , and Guillermo Iv´n Guerrero-Garc´ ıa a b a Instituto de F´ ısica, Universidad Aut´noma de San Luis Potos´ ı, o ´ lvaro Obreg´n 64, 78000 San Luis Potos´ ı, S.L.P., M´xico A o e b Facultad de Ciencias, Universidad Aut´noma de San Luis Potos´ ı, o Av. Chapultepec 1570, Privadas del Pedregal, 78295 San Luis Potos´ ı, S.L.P., M´xico e

* Corresponding Author: enrique.gonzalez@uaslp.mx

## 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., ∇ 2 φ /vector ( r ) = -g . 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 g also depends on φ /vector ( r ) and/or on its derivatives, apart from an explicit dependence on /vector r . In this context, below, we introduce and exemplify an iterative Monte

Carlo method that allows to treat numerically ordinary differential equations of the form d y 2 ( x ) dx 2 = F ( x, y ( x , ) dy ( x ) dx ) , 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 y x ( ), with Dirichlet boundary conditions, is generally stated as:

$$E \left ( x, y \left ( x \right ), \frac { d y \left ( x \right ) } { d x }, \frac { d ^ { 2 } y \left ( x \right ) } { d x ^ { 2 } } \right ) = 0,$$

such as y a ( ) = y a and y b ( ) = y b , /owner a &lt; b .

If Eq. (1) can be put in the form

$$\frac { d ^ { 2 } y \left ( x \right ) } { d x ^ { 2 } } = F \left ( x, y \left ( x \right ), \frac { d y \left ( x \right ) } { d x } \right ),$$

the Iterative Monte Carlo Method (IMCM), discussed in the following, can be utilised to solve Eq. (2), for y a ( ) = y a and y b ( ) = y b . As a preliminary, let us, then, consider a simpler case of such equation, i.e.,

$$\frac { d ^ { 2 } y \left ( x \right ) } { d x ^ { 2 } } = F \left ( x \right )$$

(notice that the RHS term of the previous equation only contains an explicit dependence on the variable x ). The usual random walk-based Monte Carlo algorithm to solve Eq. (3) comprises the next steps [6]:

STEP 1: Given a uniform mesh { x 0 = a, x 1 , x 2 , ..., x N -1 , x N = b } , for each free point x i , such that i = 1 , ..., N -1, generate a large number, K , of random walks, which start at x i and end when hitting an absorbing boundary site, either x 0 or x N .

STEP 2: If the j -th walk arrives at the boundary after m ( i,j ) steps and has visited the sequence of locations { P ( i,j ) 0 = x , P i ( i,j ) 1 , P ( i,j ) 2 , ..., P ( i,j ) m ( i,j ) = x ( i,j ) end } , calculate the Monte Carlo estimator for y x ( i ) from

$$y ( x _ { i } ) = \frac { 1 } { K } \sum _ { j = 1 } ^ { K } y ( x _ { e n d } ^ { ( i, j ) } ) - \frac { h ^ { 2 } } { 2 K } \sum _ { j = 1 } ^ { K } \left \{ \sum _ { l = 0 } ^ { \mu ^ { ( i, j ) } } F ( P _ { l } ^ { ( i, j ) } ) \right \}, \quad \quad ( 4 ) \\ \text{here } \mu ^ { ( i, j ) } _ { \cdot \cdot \cdot } = m _ { \cdot \cdot \cdot \cdot } ^ { ( i, j ) } - 1 \text{ and } y ( x _ { e n d } ^ { ( i, j ) } ) \text{ being the value, } y _ { a } \text{ or } y _ { b }, \text{ of the function}.$$

where µ ( i,j ) = m ( i,j ) -1 and y x ( ( i,j ) end ) being the value, y a or y b , of the function y x ( ) at the particular absorbing boundary point, x 0 = a or x N = b , reached by the j -th random walk initiated at x i . Observe that the approximate values of y x ( ) for all the points in the mesh can be obtained for any arbitrary

sequence of the indices in the set 1 , ..., N -1, 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 y ( x ) and y ′ ( x ) = dy ( x ) dx . However, and given that, at the end, these two sought functions depend on x , the first version of our Iterative Monte Carlo Method ( IMCM ) proposes the use of an initial guess function, y [0] ( x ) , with a well-known and simple dependence on x in the interval [ x , x 0 N ] and that fulfills the associated boundary conditions y [0] ( x 0 ) = y 0 = y a and y [0] ( x N ) = y N = y b ; 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 y [1] ( x ) at x i :

$$\nu \, \text{www.y.co.jp} \, \text{www.w} \, y \, \quad \nu _ { 1 } \, \text{www.w} _ { 2 } \,. \\ y ^ { [ 1 ] } ( x _ { i } ) = \frac { 1 } { K } \sum _ { j = 1 } ^ { K } y ( x _ { e n d } ^ { ( i, j ) } ) - \\ \frac { h ^ { 2 } } { 2 K } \sum _ { j = 1 } ^ { K } \left \{ \sum _ { l = 0 } ^ { \mu ^ { ( i, j ) } } F \left ( P _ { l } ^ { ( i, j ) }, y ^ { [ 0 ] } ( P _ { l } ^ { ( i, j ) } ), y ^ { [ 0 ] ^ { \prime } } ( P _ { l } ^ { ( i, j ) } ) \right ) \right \}, \\ \text{or } i \, = \, 1, \dots, N \, - \, 1, \, \text{with } y ^ { [ 0 ] ^ { \prime } } ( x _ { i } ) \, \text{ being the numerical derivative of } y ^ { [ 0 ] } ( x ) \\ \text{at each node (which can be calculated through any finite difference formula).}$$

for i = 1 , ..., N -1, with y [0] ′ ( x i ) being the numerical derivative of y [0] ( x ) at each node (which can be calculated through any finite difference formula). After determining the values of y [1] ( x ) for all the free nodes, the corresponding numerical derivative has to be got to proceed with successive approximations of y x ( ), obtained via the iteration prescription

$$\nu \nu, \quad \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \  \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \ \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \quad \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds}\, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \, \text{sounds} \\ \text{complemented with a finite difference approximation of its derivative. } \text{ in a} \\ \text{similar way to the scheme associated to Eq. } \text{ (4), the prior iteration expression}$$

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 y [ s +1] ( x i ) 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, y [0] ( x ), in [ x , x 0 N ], that satisfies the boundary conditions of the problem and, then, initialise ˜( y x i ) = y [0] ( x i ), for i = 0 , ..., N .

STEP 2: For each free node, perform K absorbing random walks starting at x i to obtain and update immediately ˜( y x i ) employing the formula

$$\tilde { y } ( x _ { i } ) & = \frac { 1 } { K } \sum _ { j = 1 } ^ { K } y ( x _ { e n d } ^ { ( i, j ) } ) - \\ \frac { h ^ { 2 } } { 2 K } \sum _ { j = 1 } ^ { K } \left \{ \sum _ { l = 0 } ^ { \mu ( i, j ) } F \left ( P _ { l } ^ { ( i, j ) }, \tilde { y } ( P _ { l } ^ { ( i, j ) } ), \tilde { y } ^ { \prime } ( P _ { l } ^ { ( i, j ) } ) \right ) \right \}. \\ \text{EP 3: Right after calculating $\tilde{y}(x_{i}$)} \text{ in STEP 2, differentiate it numerically} \\ \text{produce $\tilde{y}^{\prime}(x_{i}$)}.$$

STEP 3: Right after calculating ˜( y x i ) in STEP 2, differentiate it numerically to produce ˜ ( y ′ x i ).

STEP 4: Repeat entirely the previous steps 2 and 3 to complete a cycle of T 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).

## 3 An example of application

Let us consider the following nonlinear second-order ODE of the type d y 2 ( x ) dx 2 = F ( x, y ( x )) [7]:

$$\frac { d ^ { 2 } y \left ( x \right ) } { d x ^ { 2 } } = \pi ^ { 2 } \sin \left ( \pi x \right ) + e ^ { \sin \left ( \pi x \right ) + x } - e ^ { - y \left ( x \right ) },$$

with y (0) = 0 and y (1) = -1. Eq. (8) has the exact solution

$$y ^ { e x a c t } \left ( x \right ) = - \sin \left ( \pi x \right ) - x,$$

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

$$F \left ( x, y \left ( x \right ) \right ) = \pi ^ { 2 } \sin \left ( \pi x \right ) + e ^ { \sin \left ( \pi x \right ) + x } - e ^ { - y \left ( x \right ) }, \text{ \quad \ \ } ( 1 0 )$$

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 [0 1], , with N = 100 equally-sized partitions, combined with K = 10 5 random walks, to estimate the value of y x ( ) at each non-boundary node of the grid via Eq. (7). The error in the s -th iteration, E ( s ) is calculated from

$$E ^ { ( s ) } & = \sqrt { \sum _ { i = 0 } ^ { N } \left [ y ^ { e x a c t } \left ( x _ { i } \right ) - \tilde { y } ^ { ( s ) } \left ( x _ { i } \right ) \right ] ^ { 2 } } \,, & ( 1 1 ) \\ \intertext { x \, \text{when the $a$ to $b$ to $m$} \, \text{the above one} \, \text{in} \, \text{or} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$}$$

with ˜ y ( s ) ( x ) being the s -th IMCM approximation. The corresponding errors for the IMCM results for iterations 1 and 10, plotted in Fig. 2, are 1 001 . × 10 0 and 2 787 . × 10 -2 , respectively.

## 4 Conclusions

In this paper, we have proposed an Iterative Monte Carlo Method to solve second-order ODEs of the general form d y 2 ( x ) dx 2 = F ( x, y ( x , ) dy ( x ) dx ) , 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´ zniakowski, Breaking Intractability, Scientific American 270 , Issue 1, 102 (1994). doi:10.1038/scientificamerican0194102
- [3] J. M. Hammersley and D. C. Handscomb, Monte Carlo Methods (Methuen &amp; Co. Ltd., 1964).
- [4] R. Courant, K. Friedrichs, and H. Lewy, ¨ ber U 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 American 220 , Issue 3, 66 (1969). doi:10.1038/scientificamerican036966
- [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 Communication 3 , 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 Physics 16 , 491 (2014). doi:10.4208/cicp.270713.280214a

Algorithm 1 Pseudocode for the Iterative Monte Carlo Method to solve a second-order ODE with Dirichlet boundary conditions.

```
Algorithm 1 Pseudocode for the Iterative Monte Carlo Method to solve a
second-order ODE with Dirichlet boundary conditions.
  1: Input data: x0, xN, N, K, T, j'[0] = {j0, j}[0(1),..., j'[0](N - 1), jN}
  2: Result: `j = {j}(0), j'(1),..., j'(N - 1), j'(N)}
  3: Define function:  F(a, b, c)
  4: Initialization: `j = j'[0], h = (xN - x0)/N
  5: x(0) = x0
  6: x(N) = xN
  7: Mesh and initial guess for the derivative at the free points:
  8: for l = 1 to N - 1 do
  9:     x(l) = x0 + (l * h)
 10:     j'(l) = (j'(l + 1) - j'(l - 1))/(2h)
 11: end for
 12: Iteration Loop:
 13: for s = 1 to T do
 14:     for i = 1 to N - 1 do
 15:         r = 0
 16:         for j = 1 to K do
 17:             sum = 0
 18:             index = i
 19:             p = x(index)
 20:             while (x0 < p < xN) do
 21:                 a = x(index)
 22:                 b = j'(index)
 23:                 c = j'(index)
 24:                 sum = sum + F(a, b, c)
 25:                 v = rnd()
 26:                 if (v < 0.5) then
 27:                     index = index - 1
 28:                 else
 29:                      index = index + 1
 30:                 end if
 31:                 p = x(index)
 32:             end while
 33:             r = r + ((j'(index) - sum)/K)
 34:        end for
 35:        j'(i) = r
 36:        j'(i) = (j'(i + 1) - j'(i - 1))/(2h)
 37:     end for
 38: end for
 39: Output: `j` and `j`
------------------------------------------------------------------------
Figure 1: Pseudocode for the Iterative Monte Carlo Method (successive dis-
placement version).
                                                      7
```

Figure 1: Pseudocode for the Iterative Monte Carlo Method (successive displacement version).

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00742v1-figure-e2798605-abb3-482f-ae8a-63b515c33a7b.jpg)

x

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).
Martín Chávez-Páez Enrique González-Tovar Guillermo Iván Guerrero-García
An iterative Monte Carlo method to solve nonlinear second-order differential equations
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.