Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                

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., 2ϕ(r)=gsuperscript2italic-ϕ𝑟𝑔\nabla^{2}\phi\left(\vec{r}\,\right)=-g∇ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_ϕ ( over→ start_ARG italic_r end_ARG ) = - italic_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𝑔gitalic_g also depends on ϕ(r)italic-ϕ𝑟\phi\left(\vec{r}\,\right)italic_ϕ ( over→ start_ARG italic_r end_ARG ) and/or on its derivatives, apart from an explicit dependence on r𝑟\vec{r}over→ start_ARG italic_r end_ARG. In this context, below, we introduce and exemplify an iterative Monte Carlo method that allows to treat numerically ordinary differential equations of the form d2y(x)dx2=F(x,y(x),dy(x)dx)superscript𝑑2𝑦𝑥𝑑superscript𝑥2𝐹𝑥𝑦𝑥𝑑𝑦𝑥𝑑𝑥\dfrac{d^{2}y\left(x\right)}{dx^{2}}=F\left(x,y\left(x\right),\dfrac{dy\left(x% \right)}{dx}\right)divide start_ARG italic_d start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_y ( italic_x ) end_ARG start_ARG italic_d italic_x start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG = italic_F ( italic_x , italic_y ( italic_x ) , divide start_ARG italic_d italic_y ( italic_x ) end_ARG start_ARG italic_d italic_x end_ARG ), 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)𝑦𝑥y(x)italic_y ( italic_x ), with Dirichlet boundary conditions, is generally stated as:

E(x,y(x),dy(x)dx,d2y(x)dx2)=0,𝐸𝑥𝑦𝑥𝑑𝑦𝑥𝑑𝑥superscript𝑑2𝑦𝑥𝑑superscript𝑥20E\left(x,y\left(x\right),\dfrac{dy\left(x\right)}{dx},\dfrac{d^{2}y\left(x% \right)}{dx^{2}}\right)=0,italic_E ( italic_x , italic_y ( italic_x ) , divide start_ARG italic_d italic_y ( italic_x ) end_ARG start_ARG italic_d italic_x end_ARG , divide start_ARG italic_d start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_y ( italic_x ) end_ARG start_ARG italic_d italic_x start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ) = 0 , (1)

such as y(a)=ya𝑦𝑎subscript𝑦𝑎y(a)=y_{a}italic_y ( italic_a ) = italic_y start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT and y(b)=yb𝑦𝑏subscript𝑦𝑏y(b)=y_{b}italic_y ( italic_b ) = italic_y start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT, a<bcontainsabsent𝑎𝑏\ni\,a<b∋ italic_a < italic_b.

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

d2y(x)dx2=F(x,y(x),dy(x)dx),superscript𝑑2𝑦𝑥𝑑superscript𝑥2𝐹𝑥𝑦𝑥𝑑𝑦𝑥𝑑𝑥\dfrac{d^{2}y\left(x\right)}{dx^{2}}=F\left(x,y\left(x\right),\dfrac{dy\left(x% \right)}{dx}\right),divide start_ARG italic_d start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_y ( italic_x ) end_ARG start_ARG italic_d italic_x start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG = italic_F ( italic_x , italic_y ( italic_x ) , divide start_ARG italic_d italic_y ( italic_x ) end_ARG start_ARG italic_d italic_x end_ARG ) , (2)

the Iterative Monte Carlo Method (IMCM), discussed in the following, can be utilised to solve Eq. (2), for y(a)=ya𝑦𝑎subscript𝑦𝑎y(a)=y_{a}italic_y ( italic_a ) = italic_y start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT and y(b)=yb𝑦𝑏subscript𝑦𝑏y(b)=y_{b}italic_y ( italic_b ) = italic_y start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT. As a preliminary, let us, then, consider a simpler case of such equation, i.e.,

d2y(x)dx2=F(x)superscript𝑑2𝑦𝑥𝑑superscript𝑥2𝐹𝑥\dfrac{d^{2}y\left(x\right)}{dx^{2}}=F\left(x\right)divide start_ARG italic_d start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_y ( italic_x ) end_ARG start_ARG italic_d italic_x start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG = italic_F ( italic_x ) (3)

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

STEP 1: Given a uniform mesh {x0=a,x1,x2,,xN1,xN=b}formulae-sequencesubscript𝑥0𝑎subscript𝑥1subscript𝑥2subscript𝑥𝑁1subscript𝑥𝑁𝑏\left\{x_{0}=a,x_{1},x_{2},...,x_{N-1},x_{N}=b\right\}{ italic_x start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = italic_a , italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , … , italic_x start_POSTSUBSCRIPT italic_N - 1 end_POSTSUBSCRIPT , italic_x start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT = italic_b }, for each free point xisubscript𝑥𝑖x_{i}italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT, such that i=1,,N1𝑖1𝑁1i=1,...,N-1italic_i = 1 , … , italic_N - 1, generate a large number, K𝐾Kitalic_K, of random walks, which start at xisubscript𝑥𝑖x_{i}italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT and end when hitting an absorbing boundary site, either x0subscript𝑥0x_{0}italic_x start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT or xNsubscript𝑥𝑁x_{N}italic_x start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT.

STEP 2: If the j𝑗jitalic_j-th walk arrives at the boundary after m(i,j)superscript𝑚𝑖𝑗m^{(i,j)}italic_m start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT steps and has visited the sequence of locations {P0(i,j)=xi,P1(i,j),P2(i,j),,Pm(i,j)(i,j)=xend(i,j)}formulae-sequencesuperscriptsubscript𝑃0𝑖𝑗subscript𝑥𝑖superscriptsubscript𝑃1𝑖𝑗superscriptsubscript𝑃2𝑖𝑗superscriptsubscript𝑃superscript𝑚𝑖𝑗𝑖𝑗superscriptsubscript𝑥𝑒𝑛𝑑𝑖𝑗\left\{P_{0}^{(i,j)}=x_{i},P_{1}^{(i,j)},P_{2}^{(i,j)},...,P_{m^{(i,j)}}^{(i,j% )}=x_{end}^{(i,j)}\right\}{ italic_P start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT = italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT , italic_P start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT , italic_P start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT , … , italic_P start_POSTSUBSCRIPT italic_m start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT = italic_x start_POSTSUBSCRIPT italic_e italic_n italic_d end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT }, calculate the Monte Carlo estimator for y(xi)𝑦subscript𝑥𝑖y(x_{i})italic_y ( italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) from

y(xi)=1Kj=1Ky(xend(i,j))h22Kj=1K{l=0μ(i,j)F(Pl(i,j))},𝑦subscript𝑥𝑖1𝐾superscriptsubscript𝑗1𝐾𝑦superscriptsubscript𝑥𝑒𝑛𝑑𝑖𝑗superscript22𝐾superscriptsubscript𝑗1𝐾superscriptsubscript𝑙0superscript𝜇𝑖𝑗𝐹superscriptsubscript𝑃𝑙𝑖𝑗y(x_{i})=\dfrac{1}{K}\sum_{j=1}^{K}y(x_{end}^{(i,j)})-\dfrac{h^{2}}{2K}\sum_{j% =1}^{K}\left\{\sum_{l=0}^{\mu^{(i,j)}}F(P_{l}^{(i,j)})\right\},italic_y ( italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) = divide start_ARG 1 end_ARG start_ARG italic_K end_ARG ∑ start_POSTSUBSCRIPT italic_j = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_K end_POSTSUPERSCRIPT italic_y ( italic_x start_POSTSUBSCRIPT italic_e italic_n italic_d end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT ) - divide start_ARG italic_h start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_K end_ARG ∑ start_POSTSUBSCRIPT italic_j = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_K end_POSTSUPERSCRIPT { ∑ start_POSTSUBSCRIPT italic_l = 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_μ start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT end_POSTSUPERSCRIPT italic_F ( italic_P start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT ) } , (4)

where μ(i,j)=m(i,j)1superscript𝜇𝑖𝑗superscript𝑚𝑖𝑗1\mu^{(i,j)}=m^{(i,j)}-1italic_μ start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT = italic_m start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT - 1 and y(xend(i,j))𝑦superscriptsubscript𝑥𝑒𝑛𝑑𝑖𝑗y(x_{end}^{(i,j)})italic_y ( italic_x start_POSTSUBSCRIPT italic_e italic_n italic_d end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT ) being the value, yasubscript𝑦𝑎y_{a}italic_y start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT or ybsubscript𝑦𝑏y_{b}italic_y start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT, of the function y(x)𝑦𝑥y(x)italic_y ( italic_x ) at the particular absorbing boundary point, x0=asubscript𝑥0𝑎x_{0}=aitalic_x start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = italic_a or xN=bsubscript𝑥𝑁𝑏x_{N}=bitalic_x start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT = italic_b, reached by the j𝑗jitalic_j-th random walk initiated at xisubscript𝑥𝑖x_{i}italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT. Observe that the approximate values of y(x)𝑦𝑥y(x)italic_y ( italic_x ) for all the points in the mesh can be obtained for any arbitrary sequence of the indices in the set 1,,N11𝑁1{1,...,N-1}1 , … , italic_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)𝑦𝑥y\left(x\right)italic_y ( italic_x ) and y(x)=dy(x)dxsuperscript𝑦𝑥𝑑𝑦𝑥𝑑𝑥y^{\prime}(x)=\dfrac{dy\left(x\right)}{dx}italic_y start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ( italic_x ) = divide start_ARG italic_d italic_y ( italic_x ) end_ARG start_ARG italic_d italic_x end_ARG. However, and given that, at the end, these two sought functions depend on x𝑥xitalic_x, the first version of our Iterative Monte Carlo Method (IMCM) proposes the use of an initial guess function, y[0](x)superscript𝑦delimited-[]0𝑥y^{[0]}(x)italic_y start_POSTSUPERSCRIPT [ 0 ] end_POSTSUPERSCRIPT ( italic_x ), with a well-known and simple dependence on x𝑥xitalic_x in the interval [x0,xN]subscript𝑥0subscript𝑥𝑁[x_{0},x_{N}][ italic_x start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT , italic_x start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT ] and that fulfills the associated boundary conditions y[0](x0)=y0=yasuperscript𝑦delimited-[]0subscript𝑥0subscript𝑦0subscript𝑦𝑎y^{[0]}(x_{0})=y_{0}=y_{a}italic_y start_POSTSUPERSCRIPT [ 0 ] end_POSTSUPERSCRIPT ( italic_x start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) = italic_y start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = italic_y start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT and y[0](xN)=yN=ybsuperscript𝑦delimited-[]0subscript𝑥𝑁subscript𝑦𝑁subscript𝑦𝑏y^{[0]}(x_{N})=y_{N}=y_{b}italic_y start_POSTSUPERSCRIPT [ 0 ] end_POSTSUPERSCRIPT ( italic_x start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT ) = italic_y start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT = italic_y start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT; 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)superscript𝑦delimited-[]1𝑥y^{[1]}(x)italic_y start_POSTSUPERSCRIPT [ 1 ] end_POSTSUPERSCRIPT ( italic_x ) at xisubscript𝑥𝑖x_{i}italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT:

y[1](xi)=1Kj=1Ky(xend(i,j))h22Kj=1K{l=0μ(i,j)F(Pl(i,j),y[0](Pl(i,j)),y[0](Pl(i,j)))},superscript𝑦delimited-[]1subscript𝑥𝑖1𝐾superscriptsubscript𝑗1𝐾𝑦superscriptsubscript𝑥𝑒𝑛𝑑𝑖𝑗superscript22𝐾superscriptsubscript𝑗1𝐾superscriptsubscript𝑙0superscript𝜇𝑖𝑗𝐹superscriptsubscript𝑃𝑙𝑖𝑗superscript𝑦delimited-[]0superscriptsubscript𝑃𝑙𝑖𝑗superscriptsuperscript𝑦delimited-[]0superscriptsubscript𝑃𝑙𝑖𝑗\begin{multlined}y^{[1]}(x_{i})=\dfrac{1}{K}\sum_{j=1}^{K}y(x_{end}^{(i,j)})\,% -\\ \dfrac{h^{2}}{2K}\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\},% \end{multlined}y^{[1]}(x_{i})=\dfrac{1}{K}\sum_{j=1}^{K}y(x_{end}^{(i,j)})\,-% \\ \dfrac{h^{2}}{2K}\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\},start_ROW start_CELL italic_y start_POSTSUPERSCRIPT [ 1 ] end_POSTSUPERSCRIPT ( italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) = divide start_ARG 1 end_ARG start_ARG italic_K end_ARG ∑ start_POSTSUBSCRIPT italic_j = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_K end_POSTSUPERSCRIPT italic_y ( italic_x start_POSTSUBSCRIPT italic_e italic_n italic_d end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT ) - end_CELL end_ROW start_ROW start_CELL divide start_ARG italic_h start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_K end_ARG ∑ start_POSTSUBSCRIPT italic_j = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_K end_POSTSUPERSCRIPT { ∑ start_POSTSUBSCRIPT italic_l = 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_μ start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT end_POSTSUPERSCRIPT italic_F ( italic_P start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT , italic_y start_POSTSUPERSCRIPT [ 0 ] end_POSTSUPERSCRIPT ( italic_P start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT ) , italic_y start_POSTSUPERSCRIPT [ 0 ] end_POSTSUPERSCRIPT start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ( italic_P start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT ) ) } , end_CELL end_ROW (5)

for i=1,,N1𝑖1𝑁1i=1,...,N-1italic_i = 1 , … , italic_N - 1, with y[0](xi)superscriptsuperscript𝑦delimited-[]0subscript𝑥𝑖{y^{[0]}}^{\prime}(x_{i})italic_y start_POSTSUPERSCRIPT [ 0 ] end_POSTSUPERSCRIPT start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ( italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) being the numerical derivative of y[0](x)superscript𝑦delimited-[]0𝑥y^{[0]}(x)italic_y start_POSTSUPERSCRIPT [ 0 ] end_POSTSUPERSCRIPT ( italic_x ) at each node (which can be calculated through any finite difference formula). After determining the values of y[1](x)superscript𝑦delimited-[]1𝑥y^{[1]}(x)italic_y start_POSTSUPERSCRIPT [ 1 ] end_POSTSUPERSCRIPT ( italic_x ) for all the free nodes, the corresponding numerical derivative has to be got to proceed with successive approximations of y(x)𝑦𝑥y(x)italic_y ( italic_x ), obtained via the iteration prescription

y[s+1](xi)=1Kj=1Ky(xend(i,j))h22Kj=1K{l=0μ(i,j)F(Pl(i,j),y[s](Pl(i,j)),y[s](Pl(i,j)))},superscript𝑦delimited-[]𝑠1subscript𝑥𝑖1𝐾superscriptsubscript𝑗1𝐾𝑦superscriptsubscript𝑥𝑒𝑛𝑑𝑖𝑗superscript22𝐾superscriptsubscript𝑗1𝐾superscriptsubscript𝑙0superscript𝜇𝑖𝑗𝐹superscriptsubscript𝑃𝑙𝑖𝑗superscript𝑦delimited-[]𝑠superscriptsubscript𝑃𝑙𝑖𝑗superscriptsuperscript𝑦delimited-[]𝑠superscriptsubscript𝑃𝑙𝑖𝑗\begin{multlined}y^{[s+1]}(x_{i})=\dfrac{1}{K}\sum_{j=1}^{K}y(x_{end}^{(i,j)})% \,-\\ \dfrac{h^{2}}{2K}\sum_{j=1}^{K}\left\{\sum_{l=0}^{\mu^{(i,j)}}F\left(P_{l}^{(i% ,j)},y^{[s]}(P_{l}^{(i,j)}),{y^{[s]}}^{\prime}(P_{l}^{(i,j)})\right)\right\},% \end{multlined}y^{[s+1]}(x_{i})=\dfrac{1}{K}\sum_{j=1}^{K}y(x_{end}^{(i,j)})\,% -\\ \dfrac{h^{2}}{2K}\sum_{j=1}^{K}\left\{\sum_{l=0}^{\mu^{(i,j)}}F\left(P_{l}^{(i% ,j)},y^{[s]}(P_{l}^{(i,j)}),{y^{[s]}}^{\prime}(P_{l}^{(i,j)})\right)\right\},start_ROW start_CELL italic_y start_POSTSUPERSCRIPT [ italic_s + 1 ] end_POSTSUPERSCRIPT ( italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) = divide start_ARG 1 end_ARG start_ARG italic_K end_ARG ∑ start_POSTSUBSCRIPT italic_j = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_K end_POSTSUPERSCRIPT italic_y ( italic_x start_POSTSUBSCRIPT italic_e italic_n italic_d end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT ) - end_CELL end_ROW start_ROW start_CELL divide start_ARG italic_h start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_K end_ARG ∑ start_POSTSUBSCRIPT italic_j = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_K end_POSTSUPERSCRIPT { ∑ start_POSTSUBSCRIPT italic_l = 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_μ start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT end_POSTSUPERSCRIPT italic_F ( italic_P start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT , italic_y start_POSTSUPERSCRIPT [ italic_s ] end_POSTSUPERSCRIPT ( italic_P start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT ) , italic_y start_POSTSUPERSCRIPT [ italic_s ] end_POSTSUPERSCRIPT start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ( italic_P start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT ) ) } , end_CELL end_ROW (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 y[s+1](xi)superscript𝑦delimited-[]𝑠1subscript𝑥𝑖y^{[s+1]}(x_{i})italic_y start_POSTSUPERSCRIPT [ italic_s + 1 ] end_POSTSUPERSCRIPT ( italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) 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)superscript𝑦delimited-[]0𝑥y^{[0]}(x)italic_y start_POSTSUPERSCRIPT [ 0 ] end_POSTSUPERSCRIPT ( italic_x ), in [x0,xN]subscript𝑥0subscript𝑥𝑁[x_{0},x_{N}][ italic_x start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT , italic_x start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT ], that satisfies the boundary conditions of the problem and, then, initialise y~(xi)=y[0](xi)~𝑦subscript𝑥𝑖superscript𝑦delimited-[]0subscript𝑥𝑖\tilde{y}(x_{i})=y^{[0]}(x_{i})over~ start_ARG italic_y end_ARG ( italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) = italic_y start_POSTSUPERSCRIPT [ 0 ] end_POSTSUPERSCRIPT ( italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ), for i=0,,N𝑖0𝑁i=0,...,Nitalic_i = 0 , … , italic_N.

STEP 2: For each free node, perform K𝐾Kitalic_K absorbing random walks starting at xisubscript𝑥𝑖x_{i}italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT to obtain and update immediately y~(xi)~𝑦subscript𝑥𝑖\tilde{y}(x_{i})over~ start_ARG italic_y end_ARG ( italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) employing the formula

y~(xi)=1Kj=1Ky(xend(i,j))h22Kj=1K{l=0μ(i,j)F(Pl(i,j),y~(Pl(i,j)),y~(Pl(i,j)))}.~𝑦subscript𝑥𝑖1𝐾superscriptsubscript𝑗1𝐾𝑦superscriptsubscript𝑥𝑒𝑛𝑑𝑖𝑗superscript22𝐾superscriptsubscript𝑗1𝐾superscriptsubscript𝑙0superscript𝜇𝑖𝑗𝐹superscriptsubscript𝑃𝑙𝑖𝑗~𝑦superscriptsubscript𝑃𝑙𝑖𝑗superscript~𝑦superscriptsubscript𝑃𝑙𝑖𝑗\begin{multlined}\tilde{y}(x_{i})=\dfrac{1}{K}\sum_{j=1}^{K}y(x_{end}^{(i,j)})% \,-\\ \dfrac{h^{2}}{2K}\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\}% .\end{multlined}\tilde{y}(x_{i})=\dfrac{1}{K}\sum_{j=1}^{K}y(x_{end}^{(i,j)})% \,-\\ \dfrac{h^{2}}{2K}\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\}.start_ROW start_CELL over~ start_ARG italic_y end_ARG ( italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) = divide start_ARG 1 end_ARG start_ARG italic_K end_ARG ∑ start_POSTSUBSCRIPT italic_j = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_K end_POSTSUPERSCRIPT italic_y ( italic_x start_POSTSUBSCRIPT italic_e italic_n italic_d end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT ) - end_CELL end_ROW start_ROW start_CELL divide start_ARG italic_h start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_K end_ARG ∑ start_POSTSUBSCRIPT italic_j = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_K end_POSTSUPERSCRIPT { ∑ start_POSTSUBSCRIPT italic_l = 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_μ start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT end_POSTSUPERSCRIPT italic_F ( italic_P start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT , over~ start_ARG italic_y end_ARG ( italic_P start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT ) , over~ start_ARG italic_y end_ARG start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ( italic_P start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_i , italic_j ) end_POSTSUPERSCRIPT ) ) } . end_CELL end_ROW (7)

STEP 3: Right after calculating y~(xi)~𝑦subscript𝑥𝑖\tilde{y}(x_{i})over~ start_ARG italic_y end_ARG ( italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) in STEP 2, differentiate it numerically to produce y~(xi)superscript~𝑦subscript𝑥𝑖\tilde{y}^{\prime}(x_{i})over~ start_ARG italic_y end_ARG start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ( italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ).

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

Refer to caption
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 d2y(x)dx2=F(x,y(x))superscript𝑑2𝑦𝑥𝑑superscript𝑥2𝐹𝑥𝑦𝑥\dfrac{d^{2}y\left(x\right)}{dx^{2}}=F\left(x,y\left(x\right)\right)divide start_ARG italic_d start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_y ( italic_x ) end_ARG start_ARG italic_d italic_x start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG = italic_F ( italic_x , italic_y ( italic_x ) ) [7]:

d2y(x)dx2=π2sin(πx)+esin(πx)+xey(x),superscript𝑑2𝑦𝑥𝑑superscript𝑥2superscript𝜋2𝜋𝑥superscript𝑒𝜋𝑥𝑥superscript𝑒𝑦𝑥\dfrac{d^{2}y\left(x\right)}{dx^{2}}=\pi^{2}\sin\left(\pi x\right)+e^{\sin% \left(\pi x\right)+x}-e^{-y\left(x\right)},divide start_ARG italic_d start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_y ( italic_x ) end_ARG start_ARG italic_d italic_x start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG = italic_π start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_sin ( italic_π italic_x ) + italic_e start_POSTSUPERSCRIPT roman_sin ( italic_π italic_x ) + italic_x end_POSTSUPERSCRIPT - italic_e start_POSTSUPERSCRIPT - italic_y ( italic_x ) end_POSTSUPERSCRIPT , (8)

with y(0)=0𝑦00y\left(0\right)=0italic_y ( 0 ) = 0 and y(1)=1𝑦11y\left(1\right)=-1italic_y ( 1 ) = - 1. Eq. (8) has the exact solution

yexact(x)=sin(πx)x,superscript𝑦𝑒𝑥𝑎𝑐𝑡𝑥𝜋𝑥𝑥y^{exact}\left(x\right)=-\sin\left(\pi x\right)-x,italic_y start_POSTSUPERSCRIPT italic_e italic_x italic_a italic_c italic_t end_POSTSUPERSCRIPT ( italic_x ) = - roman_sin ( italic_π italic_x ) - italic_x , (9)

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(x,y(x))=π2sin(πx)+esin(πx)+xey(x),𝐹𝑥𝑦𝑥superscript𝜋2𝜋𝑥superscript𝑒𝜋𝑥𝑥superscript𝑒𝑦𝑥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)},italic_F ( italic_x , italic_y ( italic_x ) ) = italic_π start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_sin ( italic_π italic_x ) + italic_e start_POSTSUPERSCRIPT roman_sin ( italic_π italic_x ) + italic_x end_POSTSUPERSCRIPT - italic_e start_POSTSUPERSCRIPT - italic_y ( italic_x ) end_POSTSUPERSCRIPT , (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 [0,1]01\left[0,1\right][ 0 , 1 ], with N=100𝑁100N=100italic_N = 100 equally-sized partitions, combined with K=105𝐾superscript105K=10^{5}italic_K = 10 start_POSTSUPERSCRIPT 5 end_POSTSUPERSCRIPT random walks, to estimate the value of y(x)𝑦𝑥y(x)italic_y ( italic_x ) at each non-boundary node of the grid via Eq. (7). The error in the s𝑠sitalic_s-th iteration, E(s)superscript𝐸𝑠E^{\left(s\right)}italic_E start_POSTSUPERSCRIPT ( italic_s ) end_POSTSUPERSCRIPT is calculated from

E(s)=i=0N[yexact(xi)y~(s)(xi)]2,superscript𝐸𝑠superscriptsubscript𝑖0𝑁superscriptdelimited-[]superscript𝑦𝑒𝑥𝑎𝑐𝑡subscript𝑥𝑖superscript~𝑦𝑠subscript𝑥𝑖2E^{\left(s\right)}=\sqrt{\sum_{i=0}^{N}\left[y^{exact}\left(x_{i}\right)-{% \tilde{y}}^{\left(s\right)}\left(x_{i}\right)\right]^{2}}\,,italic_E start_POSTSUPERSCRIPT ( italic_s ) end_POSTSUPERSCRIPT = square-root start_ARG ∑ start_POSTSUBSCRIPT italic_i = 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_N end_POSTSUPERSCRIPT [ italic_y start_POSTSUPERSCRIPT italic_e italic_x italic_a italic_c italic_t end_POSTSUPERSCRIPT ( italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) - over~ start_ARG italic_y end_ARG start_POSTSUPERSCRIPT ( italic_s ) end_POSTSUPERSCRIPT ( italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) ] start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG , (11)

with y~(s)(x)superscript~𝑦𝑠𝑥{\tilde{y}}^{\left(s\right)}\left(x\right)over~ start_ARG italic_y end_ARG start_POSTSUPERSCRIPT ( italic_s ) end_POSTSUPERSCRIPT ( italic_x ) being the s𝑠sitalic_s-th IMCM approximation. The corresponding errors for the IMCM results for iterations 1 and 10, plotted in Fig. 2, are 1.001×1001.001superscript1001.001\times 10^{0}1.001 × 10 start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT and 2.787×1022.787superscript1022.787\times 10^{-2}2.787 × 10 start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPT, respectively.

Refer to caption
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 d2y(x)dx2=F(x,y(x),dy(x)dx)superscript𝑑2𝑦𝑥𝑑superscript𝑥2𝐹𝑥𝑦𝑥𝑑𝑦𝑥𝑑𝑥\dfrac{d^{2}y\left(x\right)}{dx^{2}}=F\left(x,y\left(x\right),\dfrac{dy\left(x% \right)}{dx}\right)divide start_ARG italic_d start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_y ( italic_x ) end_ARG start_ARG italic_d italic_x start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG = italic_F ( italic_x , italic_y ( italic_x ) , divide start_ARG italic_d italic_y ( italic_x ) end_ARG start_ARG italic_d italic_x end_ARG ), 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 American 270, 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 American 220, 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 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