<figure><img src="image_1.png" /></figure>

HAL
open science

# Monte-Carlo methods in Numerical Analysis

Charbel Tannous

► To cite this version:

Charbel Tannous. Monte-Carlo methods in Numerical Analysis. Master. Analyse Numérique, UBO Brest, France. 2018, pp.21. hal-04670584

HAL Id: hal-04670584

https://hal.science/hal-04670584v1

Submitted on 12 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Monte-Carlo methods in Numerical Analysis

C. Tannous

<figure><img src="image_3.png" /></figure>

Université de Brest, Lab-STICC, CNRS-UMR 6285, F-29200 Brest, FRANCE

(Dated: August 12, 2024)

Monte-Carlo method is a powerful and versatile numerical approach used in many fields of Science and Engineering. It is described in this work as an efficient tool of Numerical Analysis tailored to deal with Integration, evaluation of Mathematical Constants, finding roots of Vector functions, solving Linear Algebra Systems, Partial Differential equations...

PACS numbers: 87.55.K-, 87.55.Gh, 41.20.Cv, 02.30.Jr, 05.40.-a

Keywords: Monte-Carlo methods, simulation, Laplace equation, Partial differential equations, statistical physics

## Contents

<table><tr><td>I. Introduction</td><td>1</td></tr><tr><td>A. Advantages of Monte-Carlo simulation</td><td>2</td></tr><tr><td>B. Limitations of Monte-Carlo simulation</td><td>3</td></tr><tr><td>C. Speeding up Monte-Carlo with Variance reduction</td><td>3</td></tr><tr><td>II. Random number generation</td><td>5</td></tr><tr><td>III. Monte-Carlo quadrature</td><td>6</td></tr><tr><td>IV. Evaluation of mathematical constants</td><td>7</td></tr><tr><td>A. Evaluation of π</td><td>7</td></tr><tr><td>B. Evaluation of φ the Golden mean</td><td>9</td></tr><tr><td>C. Evaluation of Euler e</td><td>9</td></tr><tr><td>V. Solution of 1D diffusion equation</td><td>10</td></tr><tr><td>A. Exact solution of 1D diffusion equation</td><td>10</td></tr><tr><td>VI. Solution of a 2D Laplace equation</td><td>11</td></tr><tr><td>A. Exact solution of 2D Laplace equation</td><td>12</td></tr><tr><td>VII. Linear systems and large matrix inversion</td><td>16</td></tr><tr><td>A. Monte-Carlo method in Statistical Physics</td><td>16</td></tr><tr><td>B. Deterministic Numerical Quadrature</td><td>19</td></tr><tr><td>References</td><td>20</td></tr></table>

## I. INTRODUCTION

Monte-Carlo methods constitute a powerful approach that make use of random numbers to solve numerically a problem.

An obvious reaction to the use of random numbers is that our World is inherently deterministic, then how come a random approach might resolve a deterministic picture. The answer to this thorny critique is that statistical means and tools [1, 2] are used widely in Science, Engineering, Social Sciences... to provide a tangible answer.

A famous example is Quantum Mechanics where deterministic physical quantities are obtained from probability amplitudes (Wave functions) after defining state vectors and operators, performing special averages with adequate mathematical tools [3, 4] such as Hilbert Space, Dirac Algebra, Schrödinger equation, Heisenberg representation...

Historically, the first large scale Monte-Carlo work carried out, dates back to the middle 20th century. This work pertained to studies of neutron multiplication, scattering, propagation and eventual absorption in matter... Nicholas
Metropolis, Stanislas Ulam, John Von Neumann and Enrico Fermi [5] were the first to propose and employ the Monte-Carlo method as a viable numerical technique for solving practical problems of Nuclear Physics.

A few years later, Metropolis et al. [6] developed a powerful algorithm to deal with many other problems occurring particularly in Quantum, Statistical, Atomic and Nuclear Physics.

Several earlier isolated attempts and perhaps not fully developed instances were developed, akin to a Monte-Carlo approach disguised into some form. An example is an experiment performed in middle 18th century, consisting of throwing a needle randomly on a wooden board onto which parallel lines are drawn, and evaluating $\pi$ value from the number of times, the needle intersects a line (cf. Buffon method IV A).

In fact, early descriptions of several modern Monte-Carlo techniques appear in a paper by Lord Kelvin [7], written almost 125 years ago and pertaining to a discussion of the Boltzmann equation.

Later, in the 1920s, Karl Pearson [5] used random numbers (RN) for solving complex problems in probability theory and statistics that did not possess analytic solutions. After publication of another paper on RN by R. A. Fisher and F. Yates [5], Pearson and his students used RN to obtain distributions of many complex statistics [5].

## A. Advantages of Monte-Carlo simulation

* Versatility [8]

Monte-Carlo (MC) methods can be used in many branches of Physics and Engineering such as Mechanics, Electromagnetism, Statistical Physics, Quantum Mechanics, Atomic and Nuclear Physics, Field Theory, High Energy Physics [8]... In Appendix A, application of MC to Statistical Physics is described in detail. In numerical analysis, MC methods are used in Quadrature, Linear Algebra, root finding of vector functions, Optimization, Ordinary and Partial Differential Equations...

* Curse of Dimension resolution [9]

In Quadrature problems, error is extremely dependent on system dimension (cf. Table 1). As shown in Table 1,

<table><caption>**Table 1:** Integration algorithms and associated errors in 1D and dD with N subdivisions (in integration) or steps (in Monte-Carlo). In section IC, it is shown why Monte-Carlo error is ~1/√N. In Gauss method, m is the Legendre polynomial degree [9, 10] used for integration. Simpson 1/3 and 3/8 methods are described in Appendix B.</caption><thead><tr><th>Integration</th><th>Error in 1D</th><th>Error in dD</th></tr></thead><tbody><tr><td>Trapezoidal</td><td>1/N<sup>2</sup></td><td>1/N<sup>2/d</sup></td></tr><tr><td>1/3 and 3/8 Simpson</td><td>1/N<sup>4</sup></td><td>1/N<sup>4/d</sup></td></tr><tr><td>Gauss</td><td>1/N<sup>2m-1</sup></td><td>1/N<sup>2m-1/d</sup></td></tr><tr><td>Monte-Carlo</td><td>1/√N</td><td>1/√N</td></tr></tbody></table>

when system dimension d grows, error increases since one has to use more meshing points as d increases whereas in Monte-Carlo, it is independent of d.

* Limited computational resources:

Since error is given by the number of steps, it suffices to let a computer with limited resources run for a time long enough to achieve the required steps. Thus runtime length compensates lack of resources.

* Applications to Partial Differential Equations (PDE):

Application of MC methods to a 3D PDE satisfied by some function $V(x, y, z)$ defined over a 3D domain $D$ with boundary conditions over $\partial D$, an MC run is based on random walks over $D$. A random walk (RW) entails generating RN such that a displacement: $\delta_x \sim aU(0, 1)$, $\delta_y \sim bU(0, 1)$, $\delta_z \sim cU(0, 1)$ where a, b, c are typical scales along $x$, $y$, $z$ directions and $U(0, 1)$ is the uniform RNG (cf. Section II).

A typical MC algorithm is then:

1. For each $(x, y, z)$ in $D$ start a RW.

2. Generate a uniform RN $0 < r < 1$, with $r \sim U(0, 1)$ and move to next point.

3. If the RW hits the boundary $\partial D$, terminate it.

4. Evaluate $V_b$ on boundary $\partial D$ and go to step 1 to begin another RW from $(x, y, z)$.

5. After N RW, determine: $V(x, y, z) = \frac{1}{N} \sum_{k=1}^{N} V_b(k)$
## B. Limitations of Monte-Carlo simulation

* **Convergence error:**

Consider the following example adapted from Ref. [11]. Given a sequence $(x_k)$ defined by:

$$ 
\begin{aligned}
x_0 &= 1.5100050721319 \\
x_{k+1} &= \frac{3x_k^4 - 20x_k^3 + 35x_k^2 - 24}{4x_k^3 - 30x_k^2 + 70x_k - 50}
\end{aligned} 
\quad (1) 
$$

The sequence converges numerically to either 1, 2, 3, or 4 depending on machine precision. Below are given $x_k$ values starting with IEEE Quadruple precision (32 digit accuracy) using rounded decimal arithmetic of different precisions and proceeding downwards to reach 2 digit only precision. Actually IEEE double precision and higher converges to 3 whereas single precision converges to 2.

32 digits 3.000000000000000000000000000000000000000000

24 digits 3.00000000000000000000000000000000

20 digits 3.000000000000000000000000

16 digits 3.000000000000000000

12 digits 1.999999999990

10 digits 4.0000000000

8 digits 1.9999980

6 digits 1.99990

4 digits 2.097

2 digits 1.0

* **Statistical and other errors [8, 12]**

Several difficulties might arise in many different situations with different algorithms. All computers operate with limited word length and hence limited precision proper to numerical values of any variable. Truncation and round-off errors may, in some cases, lead to serious problems. In addition there are statistical errors which arise as an inherent feature of the simulation algorithm originating from a finite number of terms characterizing the generated statistical sample. These errors must be estimated and then a policy decision must be made, related to CPU time, available RAM or storage in order to reduce statistical errors. Other errors may be systematic.

## C. Speeding up Monte-Carlo with Variance reduction

The mean value $\mu$ of a function $h(X)$, with $X$ a random variable is given by:

$$ \mu = \int_{-\infty}^{+\infty} h(x)f(x)dx, \quad (2) $$

$f(x)$ is the probability density function of the random variable $X$. In Monte-Carlo, a random sequence $\{x_i\}, i = 1, 2, \dots, N$ is sampled from $f(x)$ to yield:

$$ \langle h \rangle_N = \frac{1}{N} \sum_{i=1}^{N} h(x_i). \quad (3) $$

In the limit $N \to \infty$, $\langle h \rangle_N \to \mu$. Moreover, by Central Limit Theorem [2], when $N \to \infty$, the probability density of $\langle h \rangle_N$ tends to a Gaussian (cf. eq. 5) with mean $\mu$ and variance $\sigma^2/N$, where

$$ \sigma^2 = \int_{-\infty}^{+\infty} [h(x) - \mu]^2 f(x)dx \quad (4) $$

Monte-Carlo estimate of $\mu$ is given by $\langle h \rangle_N \pm \sigma/\sqrt{N}$. This implies that $\langle h \rangle_N$ is expected to lie within $\pm \sigma/\sqrt{N}$ around $\mu$ (for $N$ sufficiently large), with a probability $p$ given by:

$$ p = \frac{\sqrt{N}}{\sigma\sqrt{2\pi}} \int_{\mu-(\sigma/\sqrt{N})}^{\mu+(\sigma/\sqrt{N})} \exp\left[-\frac{N(x-\mu)^2}{2\sigma^2}\right] dx \approx \frac{1}{\sqrt{2\pi}} \int_{-1}^{+1} \exp\left[-\frac{x^2}{2}\right] dx = 0.6826895 \quad (5) $$
Incidentally $\sigma$ is not known. Thus it is approximated by its Monte-Carlo estimate $S_N$:

$$S_N^2 = \frac{1}{N} \sum_{i=1}^{N} h^2(x_i) - \left[ \frac{1}{N} \sum_{i=1}^{N} h(x_i) \right]^2. \qquad (6)$$

$\pm S_N/\sqrt{N}$ is the statistical error. Notice that sample size $N$ must be large enough for the above estimates of error and mean to be acceptable. Another valuable test is to check whether sample mean has reached a Gaussian distribution.

The error is of order $1/\sqrt{N}$ with $N$ the number of MC runs. While an MC method allows to avoid dimension curse, it is extremely slow otherwise, since for a reasonable accuracy such as $10^{-6}$, the number of steps $N$ should be about $10^{12}$.

The remedy is to resort to variance reduction techniques [13] that minimize variance without altering the mean. The basic principle behind these techniques is called importance sampling [8, 12] and the procedure terminology is called "importance-based".

Importance sampling helps select from sample space, regions providing larger contributions, instead of those leading to small irrelevant values to the Monte-Carlo sum. Suppose $\{x_i\}, i = 1, 2, \dots, N$ are sampled from an importance density $g(x)$ instead of $f(x)$. In order to preserve the mean, we define a modified function $H(x)$ given by:

$$H(x) = h(x) \frac{f(x)}{g(x)}. \qquad (7)$$

When $H$ mean is evaluated over the importance density $g(x)$ it is equal to the mean of the original function $h(x)$ over density $f(x)$:

$$\mu(H) = \int_{-\infty}^{+\infty} H(x)g(x)dx = \int_{-\infty}^{+\infty} h(x)f(x)dx = \mu(h) \qquad (8)$$

Thus $\{x_i\}, i = 1, 2, \dots, N$ sampled from $g(x)$, allow to calculate $f(x_i)/g(x_i)$. The importance-based estimate of $\mu$ is given by:

$$\langle H \rangle_N = \frac{1}{N} \sum_{i=1}^{N} h(x_i) \frac{f(x_i)}{g(x_i)}. \qquad (9)$$

When $N \to \infty$, $\langle H \rangle_N \to \mu$.

It is now possible to calculate the statistical error associated with $\langle H \rangle_N$. The second moment [1, 2] is a valuable measure, since mean is preserved under importance sampling [12].

The importance-based second moment [13] using $g(x)$ is given by:

$$\langle H^2 \rangle_g = \int_{-\infty}^{+\infty} H^2(x)g(x)dx = \int_{-\infty}^{+\infty} \left[ \frac{f(x)}{g(x)} \right] h^2(x)f(x)dx \qquad (10)$$

The ordinary second moment [13] using $f(x)$ is given by:

$$\langle h^2 \rangle_f = \int_{-\infty}^{+\infty} h^2(x)f(x)dx \qquad (11)$$

Thus choosing an adequate importance function $g(x)$ and ensuring that, on average, the ratio $[f(x)/g(x)] \ll 1$, it is possible to make $\langle H^2 \rangle_g \ll \langle h^2 \rangle_f$. Consequently, sampling from an importance density helps estimate the mean with a much better statistical accuracy for a given sample size, since the corresponding variance has been reduced.

In statistical physics, importance sampling means that configurations yielding small contribution to the MC average are not considered [8, 12]. The best example, in that respect, is Metropolis Algorithm (cf. Appendix A) that evaluates an energy change $\Delta E$ of a magnetic spin system, resulting from a single spin state being altered. If $\Delta E$ corresponding to a Boltzmann factor $\exp(-\Delta E/k_B T)$ ($k_B$ is Boltzmann constant and $T$ the spin system temperature) is small, then the spin state is not altered.

In conclusion, a general advice given by Binder et al. [12] is reproduced below with some of the steps one should consider in order to perform a successful MC simulation:
* Work with small domains and short runs.

This is useful for obtaining results quickly and checking algorithm validity. This also allows us to search rather rapidly through a wide range of parameter space to determine ranges with physically interesting behavior.

* Test the RN generator

Find some limiting cases where accurate, or exact values of certain properties can be calculated, and compare results of your algorithm with different RN sequences and/or different RN generators.

* Look at systematic variations with system size and run length

Use a wide range of sizes and run lengths and then use scaling forms to analyze data.

* Calculate error bars

Search for and estimate both statistical and systematic errors. This enables to evaluate the correctness of conclusions drawn from data evaluated.

* Make a few very long runs

Do this to ensure that there is not some hidden time scale that is much longer than anticipated.

## II. RANDOM NUMBER GENERATION

Monte-Carlo methods are based on fast and efficient Random Number Generators (RNG) in order to obtain tangible simulation results.

These are in fact Pseudo-random numbers based on a numerical algorithm making them "apparently random" since produced by a deterministic rule [14, 15]. Pseudo-random numbers are also reproducible since their generation starts from a seed that determines all subsequent "random looking" numbers and this is very important for debugging a simulation.

* Uniform RNG:

The first type of RNG is Uniform and the usual method of producing numbers that are distributed uniformly over some interval is a simple multiplicative congruential algorithm:

$$I_{j+1} = aI_j \pmod{m} \qquad (12)$$

with $a = 7^5 = 16807$ and $m = 2^{31} - 1 = 2147483647$.

There exist a more general linear congruential generator, which generates a sequence of integers $I_1, I_2, I_3$, each between 0 and $m-1$ (a large number) by the recurrence relation

$$I_{j+1} = aI_j + c \pmod{m} \qquad (13)$$

$m$ is called the modulus, and $a$ and $c$ are positive integers called the multiplier and the increment, respectively. The recurrence has a period that is obviously no greater than $m$.

If $m, a$ and $c$ are properly chosen (cf. Table 2), then the period will be of maximal length, i.e., of length $m$. It has the disadvantage that it is not free of sequential correlation over successive calls. If $k$ RN at a time are used to plot points in $k$ dimensional space (with each coordinate between 0 and 1), then the points will not tend to fill up the $k$-dimensional space, but rather will lie on $(k-1)$-dimensional planes. There will be at most about $m^{1/k}$ such planes as in the IBM failed RANDU generator case [14, 15]. Another type of RNG is based on shift registers such as R250 and R1279:

$$X_n = X_{n-103} \cdot XOR \cdot X_{n-250}, \quad X_n = X_{n-1063} \cdot XOR \cdot X_{n-1279} \qquad (14)$$

where XOR is the Boolean Exclusive Or logic operation [16]. For a complete description of RNG, several reviews exist such as Ref. [12–14].

* Non-uniform RNG:

When a non-uniform RNG is required, one may use the transformation method based on probability conservation. Thus going from a probability density $p(x)$ to another one $g(y)$ where $y = y(x)$ might be done by invoking probability conservation that is expressed by: $|p(x)dx| = |g(y)dy|$ to yield: $g(y) = p(x)|\frac{dx}{dy}|$.

To generate realizations of a random variable $X$, an acceptance-rejection method makes use of realizations of another random variable $Y$ having probability density $g_Y$ similar to the probability density of $X$, $p_X$. The basic idea is that selective subsamples from samples from one distribution are statistically equivalent to samples
<table><caption>**Table 2:** Constants for Linear congruential RNG. Adapted from Ref. [10].</caption><thead><tr><th>Overflow at</th><th>m</th><th>a</th><th>c</th><th>Overflow at</th><th>m</th><th>a</th><th>c</th></tr></thead><tbody><tr><td rowspan="2">2<sup>20</sup></td><td>6075</td><td>106</td><td>1283</td><td rowspan="4">2<sup>24</sup></td><td>14406</td><td>967</td><td>3041</td></tr><tr><td>7875</td><td>211</td><td>1663</td><td>29282</td><td>419</td><td>6173</td></tr><tr><td rowspan="2">2<sup>21</sup></td><td>7875</td><td>421</td><td>1663</td><td>12960</td><td>1741</td><td>2731</td></tr><tr><td>6075</td><td>1366</td><td>1283</td><td>14000</td><td>1541</td><td>2957</td></tr><tr><td rowspan="3">2<sup>22</sup></td><td>6655</td><td>936</td><td>1399</td><td>21870</td><td>1291</td><td>4621</td></tr><tr><td>11979</td><td>430</td><td>2531</td><td>31104</td><td>625</td><td>6571</td></tr><tr><td>11979</td><td>430</td><td>2531</td><td>139968</td><td>205</td><td>29573</td></tr><tr><td>2<sup>23</sup></td><td></td><td></td><td></td><td>2<sup>25</sup></td><td></td><td></td><td></td></tr></tbody></table>

from a different distribution. The acceptance-rejection technique is one of the most important methods in RN generation.

In order to generate a variate from a distribution [13] with density $p_X$, a random variable $Y$ is chosen so that we can easily generate realizations of it and so that its density $g_Y$ can be scaled to upper bound $p_X$ using some constant $C$, such that $Cg_Y(x) \ge p_X(x) \quad \forall x$. The density $g_Y$ is an upper-bounding density, and $Cg_Y$ is an upper-bounding function.

### Quality of RNG:

Properties of RNG have been carefully examined using several mathematical tests [12, 13, 15]. A few simple examples of such tests are:

- Uniformity test: Break up the interval between zero and one into a large number of small bins and after generating a large amount of RN, check for uniformity in the number of entries in each bin.

- Overlapping M-tuple test: Check the statistical properties of the number of times, sequences of M digit numbers appear in the RN sequence.

- Parking lot test: Plot points in an m-dimensional space where the m-coordinates of each point are determined by m-successive calls to the RNG. If there are regular structures, the algorithm should be discarded.

## III. MONTE-CARLO QUADRATURE

In deterministic quadrature (as described in the Appendix B), a function $f(x)$ is defined over $[a, b]$ divided into a set $x_i$, with $i = 0, \dots, N$, where $x_0 = a$ and $x_N = b$. The $x_i$ may be uniformly or non-uniformly distributed depending on the function properties. For instance when the function varies steeply over some sub-interval $[c, d]$ it is preferable to use a finer meshing...

In MC the meshing is random, that is $x_i$ are random variables that are picked from some probability distribution and the function $f(x)$ is thus sampled in a random fashion instead of a deterministic one [10]. The naive way to evaluate an integral in a random fashion is to express it as:

$$ \int_a^b dxf(x) \approx \frac{1}{N} \sum_{i=1}^N f(x_i) \qquad (15) $$

However there is a major paradigm shift in MC, in the way quadrature is performed since random points $(x_i, y_i)$ are generated in order to evaluate the area under $f(x)$ as shown in Fig. 1.

In order to perform integration of $f(x)$ shown in Fig. 1 function values are bound over $[a, b]$ and random points $x_i, y_i$ are generated such that $x_i \sim a + (b - a)U(a, b)$ and $y_i \sim U(0, c)$ where $U(a, b)$ (resp. $U(0, c)$) is a uniform distribution over the interval $[a, b]$ (resp. $[0, c]$). Since $y_i$ is uniformly distributed over $[0, c]$, its position is checked with the condition that $y_i \le f(x_i)$.

Then for $i = 1...N$, a counter $n$ stores the number of times, a point $(x_i, y_i) \in D$.

The probability of being in $D$ is the ratio of the surface below $f(x)$ i.e. the integral $\int_a^b dxf(x)$ and the rectangular domain $[a, b] \times [0, c]$ surface.
<figure><img src="image_4.png" /><figcaption>**Fig.1:** Function and random meshing at $x_i$, with $i = 0, \dots, N$, where $x_0 = a$ and $x_N = b$. $f(x)$ is bound by $c$ over $[a, b]$. The surface of domain $\mathcal{D}$ is bound by $[a, b]$ and the function $f(x)$ over $[a, b]$ is the integral $\int_a^b xf(x)dx$. Integration is done by throwing random points $(x_i, y_i)$ and counting the number of times they fall into $\mathcal{D}$.</figcaption></figure>

Moreover the probability of being in $\mathcal{D}$ is the total number $n$ of being in $\mathcal{D}$ over the total number of trials $N$. Thus
$$ \frac{n}{N} = \frac{S}{(b-a)c} $$
and the integral value is $S = \frac{n(b-a)c}{N}$.

# IV. EVALUATION OF MATHEMATICAL CONSTANTS

## A. Evaluation of $\pi$

### Archimedes hexagons

Historically, Archimedes developed two geometrical methods to evaluate $\pi$: From perimeter evaluation of a unit radius circle (2$\pi$) or from its surface ($\pi$). He proceeded by dividing the circle into 6 hexagons (having a unit side) and evaluated the area of 96 hexagons. Hexagons are drawn inside the unit circle and outside in order to bound the area giving $\pi$. Archimedes found the bounds [9]: $3 + \frac{10}{71} < \pi < 3 + \frac{1}{7}$ or $3.1408 < \pi < 3.1428$ whereas the exact value is $\pi = 3.1415926535...$.

While in an MC method, a half-circle of radius 1 (cf. Fig. 2 is considered and its surface $\pi/2$ is evaluated by picking $(x - i, y_i)$ points with $x_i \sim U(0, 2)$ and $y_i \sim U(0, 1)$ counting the number $n$ of times the point falls within $\mathcal{D}$. Thus the surface ratio $(\pi/2)/2$ is given by $n/N$ where $N$ is the total number of trials. Then $\pi \approx 4n/N$.

<figure><img src="image_5.png" /><figcaption>**Fig.2:** Half circle with radius 1 inscribed in $\mathcal{D} = [0, 2] \times [0, 1]$ domain.</figcaption></figure>

### Buffon needle

Another historical way to evaluate $\pi$ was developed in 1733 by Buffon, a French mathematician named Georges-Louis Leclerc aliased as Comte de Buffon [17].

Buffon used to throw "baguettes" (French bread), in public, over a pavement ruled by parallel lines drawn with chalk and count the number of times the "baguette" crosses a chalk line in order to evaluate $\pi$. After being severely criticized of wasting bread, he transformed his method into a more gentile and mundane method consisting of throwing
a needle randomly, on a wooden board ruled by parallel lines, and inferring the value of $\pi$ from the number of times the needle intersects a line.

Considering a needle of length $\ell < a$ where $a$ is the distance between parallel lines (cf. Fig. 3), we get the condition $y < \frac{\ell}{2} \sin \theta$ for the needle to intersect a line where $y$ is distance between needle center and crossed line (cf. Fig. 3).

<figure><img src="image_6.png" /><figcaption>**Fig.3:** Needle with length $\ell < a$ intersecting a line. Note the dashed symmetry lines at mid-distance between two neighboring parallel lines.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig.4:** (Color on-line) Monte-Carlo setup for the Buffon problem. The red curve is $\frac{\ell}{2} \sin \theta$ and the domain $\mathcal{D}$ is the area where $y \le \frac{\ell}{2} \sin \theta$. Avoiding redundancy consists of working with degrees and replacing $\pi$ with $180^{\circ}$. The probability that the needle intersects a line is given by the ratio of the $\mathcal{D}$ area equal to $\ell$ and the rectangular domain given by $[0, \pi] \times [0, a/2]$ with area equal to $a\pi/2$ thus $p = \frac{2\ell}{a\pi}$</figcaption></figure>

In fact, the MC method might lead to a redundancy issue, since the condition $y \le \frac{\ell}{2} \sin \theta$ requires generating $y \sim U(0, a/2)$ and $\theta \sim U(0, \pi)$. Thus one needs to know $\pi$ in order to evaluate $\pi$!

Note that $y \sim U(0, a/2)$ and not $y \sim U(0, a)$ since there is a symmetry line at mid-distance between two neighboring parallel lines.

The probability that the needle intersects a line is $p = \frac{2\ell}{a\pi}$ revealing the MC method to evaluate $\pi$ by generating points $(y_i, \theta_i)$ and evaluating the number of times, the needle cuts a parallel line. For $N$ trials and $n$ hits we have: $\pi \approx \frac{2\ell N}{an}$.

One way to avoid this redundancy pitfall, is to work in degrees instead of radians. Thus the condition $y \le \frac{\ell}{2} \sin \theta$ is satisfied by generating $y \sim U(0, a/2)$ and $\theta \sim U(0, 180^{\circ})$ and count the number the point $(\theta_i, y_i)$ falls in $\mathcal{D}$ (cf. Fig. 4).

Another way to avoid this pitfall, is to determine $\pi$ with a fixed-point algorithm since the probability $p$ that the needle intersects a line being $p = \frac{2\ell}{a\pi}$ can be expressed by the double integral:

$$p = \frac{1}{\pi a/4} \int_{0}^{\frac{\ell}{a p}} d\theta \int_{0}^{\ell \sin(\theta)/2} dy = \frac{2\ell}{a\pi} \qquad (16)$$

where the upper $\theta$ value has been redefined to half its value as well as the normalizing rectangular surface. Taking $\ell = a = 2$ gives $p = 2/\pi$ which can be substituted in the integral to become:

$$p \int_{0}^{1/p} d\theta \int_{0}^{\sin(\theta)} dy = p \qquad (17)$$
Performing integration yields the identity: $p[1 - \cos(1/p)] = p$.

This implies that the function $g(p) = p[1 - \cos(1/p)]$ is a decreasing function over the interval $[0.4, 1]$ containing the fixed point given by $p^* = 2/\pi \approx 0.63$.

An MC algorithm that finds the fixed point $p^* = 2/\pi$ starts by evaluating the probability $p$ from a surface ratio via generation of $y \sim U(0, a/2)$ and $\theta \sim U(0, 1/p)$. The MC surface ratio evaluation is within a loop that performs iterations in order to find $p^*$.

## B. Evaluation of $\phi$ the Golden mean

The Golden mean is as fascinating as $\pi$ and its use is well known in painting, music, acoustics and general art [18]... It is determined from a relation between three points $ABC$ lying on a straight line with $AB=\phi$ and $BC=1$, and $AB/BC = AC/AB$ yielding $\phi/1 = (1+\phi)/\phi$ or $\phi^2 = 1+\phi$. Geometrically, the Golden mean gives the isosceles Golden triangle with an apex angle of 36° and leading to the trigonometric relation $\phi = 2\cos(\pi/5)$ (cf. Fig. 5).

<figure><img src="image_8.png" /><figcaption>**Fig.5:** (Color on-line) Golden triangle $ABC$ with an apex angle of 36° (in black). A smaller triangle $BDC$ is obtained by drawing $BD$, the bisector of $\widehat{ABC}$ a 72° (in red) angle. Triangle $BDC$ is also Golden and geometrically similar to triangle $ABC$ since both are isosceles with an apex angle of 36° or $\pi/5$. Moreover triangle $ABD$ is isosceles with a base angle of 36°. By projecting $BD$ over $AB$ we have: $a+b = 2a\cos(\pi/5)$. $ABC$ being geometrically similar to $BDC$, we have: $AB/BC = BC/DC$ or $(a+b)/a = a/b$. Combining both relations leads to $\phi = a/b = 2\cos(\pi/5)$.</figcaption></figure>

MC determination of $\phi$ is done by considering the integral:

$$ \int_{0}^{\pi/5} dx \sin x = -\cos(\pi/5) + 1 = -\phi/2 + 1 \quad (18) $$

and evaluating it by MC quadrature to get $\phi$.

## C. Evaluation of Euler $e$

A handy integral containing $e$ as a result may be used as:

$$ \int_{0}^{\pi/2} \cos(a \tan x) dx = \frac{\pi}{2} e^{-a}, \quad a \ge 0 \quad (19) $$

Evaluating the above integral by MC quadrature for a given $a$ will provide an approximation to $e^a$.
V. SOLUTION OF 1D DIFFUSION EQUATION

Let us proceed to solve, by MC, the 1D diffusion equation:

$$ \frac{\partial C}{\partial t} = \frac{\partial^2 C}{\partial x^2} \qquad (20) $$

where $C(x, t)$ is the density of particles at position $x$ and instant $t$, starting from initial conditions (at $t = 0$) that are strongly localized such that:

$$ C(x, t = 0) = \delta(x) \qquad (21) $$

where $\delta(x)$ is Dirac delta function [19].

<figure><img src="image_9.png" /><figcaption>Fig.6: (Color on-line) Exact and Monte-Carlo diffusion profiles. At $t = 0$ concentration $C(x, 0)$ (in red) is $n_0\delta(x)$ that can be approximated by a Gaussian $\delta(x) \approx \frac{1}{\sigma_0\sqrt{2\pi}}\exp(-\frac{x^2}{2\sigma_0^2})$ centered at origin with an initial number of particles $n_0 = 100$ with a standard deviation $\sigma_0 = 1$. Time $t$ provided by the number of steps, delivers the profile $C(x, t)$. (Left) $C(x, t)$ is shown after 100 RW steps and 50 MC runs (in blue) displaying discrepancies with the analytic profile (in black). (Right) $C(x, t)$ is shown after 100 RW steps and 5000 MC runs (in blue) proving that when the number of MC runs increases, $C(x, t)$ obtained by MC becomes indistinguishable from the analytic solution (in black).</figcaption></figure>

A. Exact solution of 1D diffusion equation

Diffusion equation: $\frac{\partial P}{\partial t} = D \frac{\partial^2 P}{\partial x^2}$ is a parabolic PDE (Partial Differential Equation) [19] with $P(x, t)$ the probability density (particle concentration) and $D$ the diffusion constant, has the analytic solution [19]:

$$ P(x, t) = \frac{1}{\sqrt{4\pi Dt}} \exp\left(-\frac{(x - \langle x \rangle)^2}{4Dt}\right) \qquad (22) $$

that starts at $t = 0$ with an initial concentration $P(x, 0) = \delta(x)$ with boundary conditions $P(\pm\infty, t) = 0$. $\langle x \rangle$ is mean displacement $x$ and $P(x, t)$ is a Gaussian $\frac{1}{\sigma\sqrt{2\pi}}\exp(-\frac{(x-\langle x \rangle)^2}{2\sigma^2})$.

This yields $\sigma = \sqrt{2Dt}$ linking time and Gaussian standard deviation. Note that particle density $C(x, t)$ is $n_0P(x, t)$ where $n_0$ is the number of particles. In the spatially discrete case diffusion originates from a RW over a space-time described with $x_i = i\ell$ and temporally by $t_j = j\tau$. $\ell$ and $\tau$ are spatial and temporal steps respectively. If $p$ is the probability of moving in the right direction, then $q = 1 - p$ is the probability of moving left.

If $N$ is the total number of steps with $n_1$ to the right and $n_2$ to the left (then $N = n_1 + n_2$), the probability of having made $n_1$ steps is obtained from inspection of words such as: $ppppp...qqqq = p^{n_1}q^{n_2}$ in order to obtain all possible combinations: $\frac{N!}{n_1!n_2!}p^{n_1}q^{n_2}$.

This originates from the Combinatorial identity: $(p+q)^N = 1 = \sum_{n=0}^{N} \frac{N!}{n!(N-n)!}p^nq^{N-n}$.

Thus the probability of being at a position $x = ml$ after $N$ steps (i.e. in a time $t = N\tau$) is $P_N(m) = \frac{N!}{n_1!n_2!}p^{n_1}q^{n_2}$ with $n_1 = \frac{1}{2}(N+m)$ and $n_2 = \frac{1}{2}(N-m)$.
Mean displacement is given by $\langle x \rangle = N(p-q)\ell$ whereas mean square displacement $\langle x^2 \rangle = [N(p-q)\ell]^2 + 4pqN\ell^2$.

Standard deviation of displacement is given by: $\sigma = \sqrt{\langle x^2 \rangle - \langle x \rangle^2} = \ell\sqrt{N}$.

In the left-right equiprobable case ($p = q = \frac{1}{2}$) when position is $x = m\ell$, probability to move forward in time is given by equal left-right contributions $P(x, t + \tau) = \frac{1}{2}P(x + \ell, t) + \frac{1}{2}P(x - \ell, t)$.

Applying this condition to $P(x, t - \tau)$, yields:

$$ \frac{P(x,t) - P(x,t-\tau)}{\tau} = \frac{P(x+\ell,t-\tau) + P(x-\ell,t-\tau) - 2P(x,t-\tau)}{2\tau} \quad (23) $$

that is equivalent to:

$$ \frac{\partial P}{\partial t} \approx \frac{\ell^2}{2\tau} \frac{\partial^2 P}{\partial x^2} \quad \text{since} \quad \frac{P(x+\ell, t-\tau) + P(x-\ell, t-\tau) - 2P(x,t-\tau)}{\ell^2} \approx \frac{\partial^2 P}{\partial x^2} \quad (24) $$

This gives the diffusion coefficient value as $D = \frac{\ell^2}{2\tau}$ and consequently $\sigma = \sqrt{2Dt} = \sqrt{2\frac{\ell^2}{2\tau}N\tau} = \ell\sqrt{N}$, validating a previously obtained result.

Continuous probability $P(x, t)$ where $t = N\tau$ and $x = m\ell$ is obtained from discrete $P_N(m)$ when $N \to \infty$ resulting in $P(x, t) = \frac{1}{\sqrt{\pi Dt}} \exp(-\frac{(x-\langle x \rangle)^2}{4Dt})$.

Note that discretization leading to this probability is twice as much as the continuum solution given by $P(x, t) = \frac{1}{\sqrt{4\pi Dt}} \exp(-\frac{(x-\langle x \rangle)^2}{4Dt})$. This originates from the fact spatial discretization excludes alternatively odd and even values of position $x$.

**Note:** Extension to 1D Schrödinger equation: Schrödinger equation is a diffusion equation in imaginary time, thus the correspondence $t \to it$ with $i = \sqrt{-1}$.

## VI. SOLUTION OF A 2D LAPLACE EQUATION

Consider the 2D Laplace equation:

$$ \Delta V(x, y) = 0, \quad \text{or} \quad \frac{\partial^2 V}{\partial x^2} + \frac{\partial^2 V}{\partial y^2} = 0 \quad (25) $$

where $V(x, y)$ is an electrostatic potential at $(x, y)$ over a domain $\partial D$ defined in an $x, y$ plane: $[0, a] \times [0, b]$.

Dirichlet boundary conditions (cf. Fig. 7) are: $V(0, y) = V(a, y) = V(x, b) = 0, V(x, 0) = V_0$.

<figure><img src="image_10.png" /><figcaption>**Fig.7:** (Color on-line) Rectangular domain (in cyan) $D = [0, a] \times [0, b]$ and boundary $\partial D$ (in red). Boundary conditions are defined on $\partial D$ such that $V(0, y) = V(a, y) = V(x, b) = 0, V(x, 0) = V_0$.</figcaption></figure>

$V(x, y)$ is obtained by a Finite Difference method (FDM) (cf. Fig. 8). FDM consists of a discretization of Laplace operator and imposing boundary conditions. This results into a linear system $AX = B$ where $X$ is a vector that represents potential discretized values and $B$ is a vector representing boundary conditions. System $AX = B$ can be
solved iteratively by Jacobi or Gauss-Seidel methods [10].

In MC, discretization of Laplace operator leads to:

$$V(x, y) = \frac{1}{4} [V(x + \delta_x, y) + V(x - \delta_x, y) + V(x, y + \delta_y) + V(x, y - \delta_y)] \quad (26)$$

where $\delta_x, \delta_y$ are displacements originating from discretization. This may be interpreted in a probabilistic manner akin to RW and Diffusion equation. While in the Diffusion equation, we have space-time steps $\frac{1}{2}$, in Laplace case, we have only spatial steps $\delta_x, \delta_y$ with probability $\frac{1}{4}$.

In order to speed up the MC procedure, we make the following steps [20] to calculate $V$ at every point $(x, y)$ in domain $D$:

1. For each $(x, y)$ in $D$ start a RW.

2. Generate a uniform RN $r \sim U(0, 1)$, and move to next point.

3. If the next point is not on boundary, repeat step 2.

4. If the RW hits the boundary $\partial D$, terminate it.

5. Evaluate $V_b$ on boundary $\partial D$, go to step 1 and begin another RW.

6. After $N$ RW, determine $V(x, y) = \frac{1}{N} \sum_{k=1}^{N} V_b(k)$

Note the additional step 3 (above) that speeds up the MC procedure and that was not stated in the general algorithm of section I A.

With this algorithm, when the number of runs (10 and 100 MC runs in Fig. 8 center) $V(x, y)$ profile looks (despite fluctuations) like the FDM case. As MC runs increase (1000 in lowest part of Fig. 8), $V(x, y)$ profile gets smoother and closer to the FDM solution.

In order to validate further the solution, we proceed to obtain an analytic solution, as described below.

## A. Exact solution of 2D Laplace equation

Laplace equation that appears in Electromagnetism and Gravity problems is an elliptic PDE [10, 19].

Supposing separability of $x, y$ coordinates in 2D Laplace equation yields: $V(x, y) = X(x)Y(y)$. Substituting in Laplace PDE gives: $X'' + XY'' = 0$.

Dividing by $XY$ we get $\frac{X''}{X} + \frac{Y''}{Y} = 0$. Taking $\frac{X''}{X} = -\lambda^2$, then $\frac{Y''}{Y} = \lambda^2$.

Both equations are rewritten as: $X'' + \lambda^2 X = 0$ and $Y'' - \lambda^2 Y = 0$ with solutions given by:

$$X = A \cos \lambda x + B \sin \lambda x, \quad Y = C \cosh \lambda y + D \sinh \lambda y \quad (27)$$

Applying boundary conditions gives: $V(0, y) = V(a, y) = 0$ to $X(x)$ at $x = 0$ and $x = a$ gives $X = B \sin \lambda x$ with $\lambda = \frac{m\pi}{a}$ where $m = 1, 2, 3...$.

Boundary conditions $V(x, b) = 0, V(x, 0) = V_0$ applied to $Y(x)$ at $y = 0$ and $y = b$ yield:

$$Y = V_0 \cosh \left( \frac{m\pi}{a} y \right) - \frac{V_0}{\tanh\left(\frac{m\pi b}{a}\right)} \sinh \left( \frac{m\pi}{a} y \right) \quad (28)$$

The resulting potential becomes:

$$V(x, y) = V_0 \sum_{m=1}^{\infty} B_m \sin \left( \frac{m\pi}{a} x \right) \left[ \cosh \left( \frac{m\pi}{a} y \right) - \frac{\sinh \left( \frac{m\pi}{a} y \right)}{\tanh \left( \frac{m\pi b}{a} \right)} \right] \quad (29)$$

Clearly boundary conditions are satisfied along lines $x = 0, x = a$ as well as along $y = b$.

Accounting for boundary condition along line $y = 0$ gives:

$$V(x, y = 0) = V_0 \sum_{m=1}^{\infty} B_m \sin \left( \frac{m\pi}{a} x \right) = V_0 \quad (30)$$

which is equivalent to:

$$\sum_{m=1}^{\infty} B_m \sin \left( \frac{m\pi}{a} x \right) = 1 \quad (31)$$
<figure><img src="image_11.png" /><figcaption>**Fig.8:** (Color on-line) (Top graph) $V(x, y)$ obtained by FDM. Lower graphs display MC results after 10, 100 and 1000 steps successively. Boundary conditions are: $V(0, y) = V(a, y) = V(x, b) = 0$, $V(x, 0) = 100$ with $a = b = 50$. For 1000 MC steps $V(x, y)$ is very close to the smooth FDM solution.</figcaption></figure>
Multiplying by $\sin(\frac{n\pi}{a}x)$ and integrating over interval $[0, a]$,
results in the fact $m \neq n$ yields zero (terms orthogonality $\sin(\frac{m\pi}{a}x)$ and $\sin(\frac{n\pi}{a}x)$).
When $m = n$ we get:

$$B_n \int_0^a \sin^2\left(\frac{n\pi}{a}x\right) dx = \int_0^a \sin\left(\frac{n\pi}{a}x\right) dx \quad (32)$$

This yields $B_n \frac{a}{2} = \frac{[1-(-1)^n]}{n\pi/a}$. Thus $B_n = 0$ when $n$ is even and $B_n = \frac{4}{n\pi}$ when $n$ is odd.
Consequently, the potential is given by the Fourier series:

$$V(x, y) = \frac{4V_0}{\pi} \sum_{m=1}^{\infty} \frac{1}{m} \sin\left(\frac{m\pi}{a}x\right) \left[ \cosh\left(\frac{m\pi}{a}y\right) - \frac{\sinh\left(\frac{m\pi}{a}y\right)}{\tanh\left(\frac{m\pi b}{a}\right)} \right]; \quad m = 1, 3, 5, \dots \quad (33)$$

The first Fourier component ($m = 1$) is graphed in Fig.9.

<figure><img src="image_12.png" /><figcaption>Fig.9: (Color on-line) First Fourier ($m = 1$) component of $V(x, y)$ given by $\frac{4V_0}{\pi} \sin(\frac{\pi}{a}x) \left[ \cosh(\frac{\pi}{a}y) - \frac{\sinh(\frac{\pi}{a}y)}{\tanh(\frac{\pi b}{a})} \right]$ with $V_0 = 100, a = b = 50$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig.10: (Color on-line) Normalized $V(x, 0)/V_0$ potential evaluated with several Fourier terms, summing over $m = 1, 3, 5, 7, 9$. For $m = 1$, $V(x, 0)/V_0$ at $x = a/2 = 25$ is 1.27. As $m$ increases, smoothness is observed and progress toward 1 is steady as predicted by Fourier theory.</figcaption></figure>

Comparison of Fig.9 to FDM result of Fig.8 infers that the first Fourier component overshoots $V_0 = 100$ in $V(x, 0)$ profile since its value is $\frac{4V_0}{\pi} \sim 127$. This is due to convergence slowness of Fourier components that depend weakly
<figure><img src="image_14.png" /><figcaption>Fig.11: (Color on-line) $V(x, y)/V_0$ potential evaluated with several Fourier terms: 25 terms in top part of figure and 50 terms at lower part. Fluctuations and overshoots of $x$ are visible with 25 terms. In the 50 term case, the FDM solution is retrieved (cf top of Fig.8).</figcaption></figure>

on $m$ as $\frac{1}{m}$. It is expected that $V(x, 0)$ will ultimately to evolve with $m$ toward FDM result of Fig.8.

Indeed $V(x, 0) = \frac{4V_0}{\pi} \sum_{m=1}^{\infty} \frac{1}{m} \sin(\frac{m\pi}{a}x)$ and $B_m$ is evaluated from $\sum_{m=1}^{\infty} B_m \sin(\frac{m\pi}{a}x) = 1$.

In order to check this issue, $\frac{V(x, 0)}{V_0} = \frac{4}{\pi} \sum_{m=1}^{\infty} \frac{1}{m} \sin(\frac{m\pi}{a}x)$ is drawn (cf. Fig. 10) for several values of $m$ in order to evaluate Fourier series convergence speed. Indeed Fig.10 displays smoothness provided by Fourier series, validating the fact the potential might be represented by an infinite Fourier series.

Fourier convergence is tested when $m \to \infty$, $\tanh(\frac{m\pi b}{a}) \sim 1$ and $y$ term becomes:

$\cosh(\frac{m\pi}{a}y) - \sinh(\frac{m\pi}{a}y) = \exp(-\frac{m\pi}{a}y)$ yielding the potential:

$$V(x, y) = \frac{4V_0}{\pi} \sum_{m=1}^{\infty} \frac{1}{m} \sin\left(\frac{m\pi}{a}x\right) \exp\left(-\frac{m\pi}{a}y\right) \quad (34)$$

confirming convergence of $V(x, y)$ when $m \to \infty$ since it is bound by the series $\frac{4V_0}{\pi} \sum_{m=1}^{\infty} \exp(-\frac{m\pi}{a}y)$ approximated by $\frac{4V_0}{\pi} \frac{1}{[1-\exp(-\frac{m\pi}{a})]}$. The final confirmation is obtained numerically, by increasing the number of terms. In Fig.11 the number of terms is 25 then increased to 50. Obviously, agreement with Fig.8 top part is achieved by considering 50 Fourier terms.
VII. LINEAR SYSTEMS AND LARGE MATRIX INVERSION

The computation of the inverse of a sparse matrix or some of its elements is one of the main topics in numerical
Science and Engineering. A sparse matrix has typically a small number of non-zero elements.

A large sparse matrix can be usually inverted from its factors obtained by using special adapted techniques [21]. Although it can be tackled by reordering rows and columns, these techniques suffer from large amount of data required by LU factorization [10] and their algorithms usually demand large random access memory (RAM) as well.

Monte-Carlo algorithms are used to estimate empirically the expectation value of random variables. Using them to invert matrices consist basically in obtaining and averaging several realized values of a variable whose expectation is the required inverse [21]. Although Monte-Carlo algorithms bear a statistical error, they can be used to obtain estimates of the inverse when the matrix is too large or complicated to be inverted by using conventional or sparse algorithms [21]. These algorithms have high parallel efficiency, they do not demand a great amount of RAM and the processing time is proportional to the number of non-zero elements of the matrix to be inverted.

The most common Monte-Carlo algorithms to invert matrices [22, 23] are based on finite discrete Markov chains (cf. Appendix A. In these algorithms several random trajectories are generated to solve a linear system of equations for each matrix column. Alternatively, there are Monte-Carlo methods based on random variables whose expectation is the whole inverse instead one of its columns [22, 23].

The problem is expressed in the following way: Find the solution vector $\boldsymbol{x}$ of $A\boldsymbol{x} = \boldsymbol{b}$ where $\boldsymbol{b}$ is an $N$ dimensional vector and $A$ is $N \times N$ matrix amounts to introducing a matrix $H$ such that $I - H = A$ where $I$ is the identity matrix to result in $\boldsymbol{x} = H\boldsymbol{x} + \boldsymbol{b}$ where $H$ is an $N \times N$ matrix.

It is straightforward to transform $\boldsymbol{x} = H\boldsymbol{x} + \boldsymbol{b}$ into an iterative process in order to obtain a numerical solution through:

$$ \boldsymbol{x}^{(k+1)} = H\boldsymbol{x}^{(k)} + \boldsymbol{b} \qquad (35) $$

This iteration converges when the modulus $||H|| < 1$.
Then, starting from $\boldsymbol{x}^0 = 0$, we have:

$$ \boldsymbol{x}^{(k+1)} = (I + H + H^2 + \dots H^k)\boldsymbol{b}, \quad \lim_{k \to \infty} \boldsymbol{x}^{(k)} = \lim_{k \to \infty} H^k \boldsymbol{b} = (I - H)^{-1} \boldsymbol{b} = A^{-1} \boldsymbol{b} = \boldsymbol{x} \qquad (36) $$

Thus MC estimators are given by the Von Neumann [22, 23] series representation:

$$ x_j = b_j + \sum_{i_1=1}^{N} h_{ji_1} b_{i_1} + \sum_{i_2=1}^{N} \sum_{i_1=1}^{N} h_{ji_2} h_{i_2 i_1} b_{i_1} + \dots \qquad (37) $$

where $x_j$ is the j-th component of the solution vector $\boldsymbol{x}$. The first step in estimating $x_j$ is to introduce a discrete Markov chain (cf. Appendix A) with state space $\{1, \dots, N+1\}$, a transition matrix $P = (p_{ij})$, where $p_{ij}$ denotes the probability of transition from state $i$ to $j$, and a vector $\boldsymbol{a} = (a_j)$. Then we introduce the $N$ -dimensional vector $\boldsymbol{p} = (p_j)$, where $p_j$ is the probability that the RW (cf. Appendix A) is absorbed in state $j$, once it visits it.

Inverse matrix $C = A^{-1}$ elements $C_{rr'}$ are calculated as [24]:

$$ C_{rr'} \approx \frac{1}{N} \sum_{i=1}^{N} \left[ \sum_{(j|k_j=r')}^{N} W_j \right]_i, \quad W_j = W_{j-1} \frac{h_{k_{j-1}k_j}}{p_{k_{j-1}k_j}} \qquad (38) $$

In eq. 38, the discrete Markov chain [22, 23] is represented by a set of indices $k_j$ with transition probabilities $p_{k_{j-1}k_j}$. To compute a matrix element $C_{rr'}$, the Markov chain [24] starts with a row index $k_0 = r$, and the only $W_j$ terms included into $C_{rr'}$ sum in eq. 38 are those for which $k_j = r'$ as indicated by the logical notation $(j|k_j = r')$.
In order to compute one entire row, $N$ Markov chains are used for evaluating averages and obtaining the result.

APPENDIX A: MONTE-CARLO METHOD IN STATISTICAL PHYSICS

In a Monte-Carlo simulation, time dependence of a model is followed by an algorithm depending on a sequence of RN generated during simulation. A very large number of different problems fall into this category: in percolation an empty lattice is gradually filled with particles by placing a particle on the lattice randomly with each clock tick. A percolating cluster reaches from one boundary of a macroscopic system to the opposite one. The properties of such
objects are of interest in Science and Engineering problems such as conductivity of random mixtures, flow through porous rocks (oil mining), behavior of dilute magnets [12]...

Another example is diffusion limited aggregation (DLA) where a particle makes a RW in space, until it encounters a seed object and sticks to it producing its growth [8, 25].

In statistical physics, and specially magnetic spin systems, the most used algorithm is Metropolis. In the classic, Metropolis method [6] configurations are generated from a previous state using a transition probability which depends on the energy difference between initial and final states. The sequence of states produced follows a time ordered path, however time in this case is a non-deterministic Monte-Carlo time. Time-dependent behavior is described by a Master equation [12]:

$$ \frac{\partial P_n(t)}{\partial t} = \sum_{n \neq m} [P_n(t)W_{n \to m} - P_m(t)W_{m \to n}] \quad (A1) $$

where $P_n(t)$ is the probability of the system being in state $n$ at time $t$, and $W_{n \to m}$ is the transition rate for $n$ to $m$. At equilibrium $\frac{\partial P_n(t)}{\partial t} = 0$ and the two RHS terms must be equal. The result is known as Detailed Balance:

$$ P_n(t)W_{n \to m} = P_m(t)W_{m \to n} \quad (A2) $$

The probability of the $n$-th state occurring in a classical system is given by

$$ P_n(t) = e^{E_n/k_B T} / Z \quad (A3) $$

where $Z$ is the partition function. This probability is usually not exactly known because of the denominator. It is possible to avoid this difficulty by generating a Markov chain of states, i.e. generate each new state directly from the preceding state.

The concept of Markov chains is a central theme in Monte-Carlo simulations. A stochastic process [2] at discrete times $t_1, t_2, t_3, \dots$, is defined for a system with a finite set of possible states $S_1, S_2, S_3, \dots$, where $X_t$ is the system state at time $t$.

The conditional probability that $X_{t_n} = S_{i_n}$,

$$ P(X_{t_n} = S_{i_n} | X_{t_{n-1}} = S_{i_{n-1}}, X_{t_{n-2}} = S_{i_{n-2}}, \dots, X_{t_1} = S_{i_1}) \quad (A4) $$

given that at the previous time, the system state $X_{t_{n-1}}$ was in state $S_{i_{n-1}}$, etc. Such a memoryless process is called a Markov process [12], if this conditional probability is independent of all states but the immediate predecessor, i.e. $P(X_{t_n} = S_{i_n} | X_{t_{n-1}} = S_{i_{n-1}})$. The corresponding sequence of states $\{X_t\}$ is called a Markov chain, and the above conditional probability can be interpreted as the transition probability to move from state $S_i$ to state $S_j$,

$$ W_{ij} = W(S_i \to S_j) = P(X_{t_n} = S_j | X_{t_{n-1}} = S_i) \quad (A5) $$

We further require that

$$ W_{ij} \ge 0, \quad \sum_j W_{ij} = 1 \quad (A6) $$

as usual for transition probabilities. The total probability $P(X_{t_n} = S_j)$ that at time $t_n$ the system is in state $S_j$ as $P(X_{t_n} = S_j) = P(X_{t_n} = S_j | X_{t_{n-1}} = S_i) = W_{ij}P(X_{t_n} = S_j)$.

Master equation [12] considers the change of this probability with time $t$ (treating time as a continuous rather than a discrete variable).

Let us apply these concepts to system with states $m$ and $n$. If we produce the $n$-th state from the $m$-th state, the relative probability is the ratio of the individual probabilities and the denominator cancels. As a result, only the energy difference between the two states is needed [12]:

$$ \Delta E = E_n - E_m \quad (A7) $$

Any transition rate which satisfies detailed balance is acceptable. The first choice of rate which was used in statistical physics is the Metropolis form [6]

$$ W_{n \to m} = \tau_0^{-1} e^{-\Delta E / k_B T} \quad \text{if } \Delta E > 0; \quad W_{m \to n} = \tau_0^{-1} \quad \text{if } \Delta E < 0 \quad (A8) $$

The simplest case is the single spin-flip sampling for the simple Ising model [12] where $\tau_0$ is the time required to attempt a spin-flip. The way the Metropolis algorithm [6] is implemented can be described by a simple procedure: Metropolis importance sampling Monte-Carlo scheme [6]:
1. Choose an initial state

2. Choose a site i

3. Calculate the energy change $\Delta E$ which results if the spin at site i is flipped

4. Generate $r$, a uniform RN such that $r \sim U(0, 1)$

5. If $r < \exp(-\Delta E/k_B T)$, flip the spin

6. Go to next site and apply step (3)

Given a magnetic system with a number of spins, the physical properties are determined and added to the statistical average [25]. Note that $r \sim U(0, 1)$ and successive RN such as $r$ should be uncorrelated. States are generated with a probability proportional to above step (4) of Metropolis algorithm, once the number of states is large enough so that initial probability transients are negligible.

Then, the desired averages of a variable, simply become arithmetic averages over the entire generated sample of states.

As an example, the 2D Ising model has energy given by: $E = -J \sum_{i,j} s_i s_j$ where $J$ is Heisenberg interaction constant [3] between next neighboring spins $s_i, s_j$ ($i, j$ are 2D coordinates, corresponding each to two integer numbers along $x, y$ axes) and the sum $\sum_{i,j}$ is restricted to these spins. An Ising spin variable $s_i = \pm 1$ meaning that the Ising model is a simple representation of a strongly uniaxial anisotropic magnetic system. Total magnetization is given by $M = \sum_i s_i$.

The specific heat at constant volume [26] is given by:

$$C_V = \frac{\langle E^2 \rangle - \langle E \rangle^2}{k_B T^2} \qquad (A9)$$

where averages that are obtained with Fluctuation theory of the Canonical Ensemble [27], are identified with MC averages. MC results for average magnetization $\langle M(T) \rangle$ and specific heat $C_V(T)$ versus normalized temperature $k_B T/J$ for a 2D Ising model in Fig. 12 containing 20x20 spins. While $\langle M(T) \rangle$ is relatively smooth, $C_V(T)$ displays strong fluctuations specifically around a critical temperature $T_c$ where $\langle M(T_c) \rangle = 0$ signaling phase transition from ferromagnetic to paramagnetic state.

<figure><img src="image_15.png" /><figcaption>**Fig.12:** (Color on-line) Average magnetization $\langle M(T) \rangle$ and specific heat $C_V(T)$ versus normalized temperature $k_B T/J$ for a 2D Ising model containing 20x20 spins. The critical temperature $T_c \approx 2.3J/k_B$ where $\langle M(T_c) \rangle \approx 0$. Periodic Boundary Conditions are used such that during an MC run, if a RW leads outside the 2D domain, it is back folded in. Adapted from Giordano [25].</figcaption></figure>

Fig. 12 agrees with Onsager [12] exact solution of the 2D Ising model where magnetization and specific heat behave respectively near $T_c$ as $M \sim (T_c - T)^\beta$, $C \sim \frac{1}{|T - T_c|^\alpha}$ where $\alpha = 0, \beta = 1/8$. $\alpha = 0$ means a logarithmic divergence of the specific heat around $T_c$ as observed in Fig. 12. Onsager normalized critical temperature [25] $k_B T_c/J = 2/\ln(1 + \sqrt{2})$ corresponds to $T_c = 2.269J/k_B$ [8] agreeing with MC result $T_c \approx 2.3J/k_B$ observed in Fig. 12.
APPENDIX B: DETERMINISTIC NUMERICAL QUADRATURE

Given a function $f(x)$ defined over $[a, b]$ and known at equally spaced points $x_i$, for $i = 0, \dots, n$, where $x_0 = a$ and $x_n = b$, we qualify quadrature as local (simple) if one is interested in a given interval such as $[a, a + nh]$ where $h$ is the integration step (usually "small" since typically $h = \frac{b-a}{n}$). When the entire $[a, b]$ interval is considered for integration, quadrature is qualified as global (composite).

Newton-Cotes formalism [28] is used to derive an integration formula like:

$$ \int_{a}^{a+nh} f(x)dx = \frac{h[(1+\Delta)^n - 1]}{\ln(1+\Delta)}f(a), \quad n \neq 0 \qquad (B1) $$

where $\Delta$ is the finite difference operator [28] yielding $\Delta f(a) = f(a+h) - f(a)$

<figure><img src="image_16.png" /><figcaption>Fig.13: Deterministic integration showing a periodic set of $x_i$ coordinates and the corresponding function values that allow to determine the surface by various Newton-Cotes algorithms.</figcaption></figure>

Replacing $\Delta$ by $x$ and performing Taylor expansion around $x = 0$ yields:

$$ \frac{(1+x)^n - 1}{\ln(1+x)} \approx n + \frac{n^2 x}{2} + \frac{(2n^3 - 3n^2)x^2}{12} + \frac{(n^4 - 4n^3 + 4n^2)x^3}{24} + \frac{(6n^5 - 45n^4 + 110n^3 - 90n^2)x^4}{720} + \frac{(2n^6 - 24n^5 + 105n^4 - 200n^3 + 144n^2)x^5}{1440} + \frac{(12n^7 - 210n^6 + 1428n^5 - 4725n^4 + 7672n^3 - 5040n^2)x^6}{60480} \qquad (B2) $$

Performing the Taylor expansion of $\frac{(1+x)^n - 1}{\ln(1+x)}$ to order $n+1$, taking $n$ from 1 to 7, Table 3 is built.

<table><caption>Table 3: Taylor expansions corresponding to Newton-Cotes formulas. Adapted from Ref. [28]</caption><thead><tr><th>Order n</th><th>Expansion to order n + 1</th></tr></thead><tbody><tr><td>n = 1</td><td>1 + x/2 + x<sup>2</sup>/12 + ...</td></tr><tr><td>n = 2</td><td>2 + 2x + x<sup>3</sup>/3 + ...</td></tr><tr><td>n = 3</td><td>3 + 9x/2 + 9x<sup>2</sup>/4 + 3x<sup>3</sup>/8 + 3x<sup>4</sup>/80 + ...</td></tr><tr><td>n = 4</td><td>4 + 8x + 20x<sup>2</sup>/3 + 8x<sup>3</sup>/3 + 14x<sup>4</sup>/45 + ...</td></tr><tr><td>n = 5</td><td>5 + 25x/2 + 175x<sup>2</sup>/12 + 75x<sup>3</sup>/8 + 425x<sup>4</sup>/144 + 95x<sup>5</sup>/288 + 275x<sup>6</sup>/12096 + ...</td></tr><tr><td>n = 6</td><td>6 + 18x + 27x<sup>2</sup> + 24x<sup>3</sup> + 123x<sup>4</sup>/10 + 33x<sup>5</sup>/10 + 41x<sup>6</sup>/140 + ...</td></tr><tr><td>n = 7</td><td>7 + 49x/2 + 539x<sup>2</sup>/12 + 1225x<sup>3</sup>/24 + 26117x<sup>4</sup>/720 + 2499x<sup>5</sup>/160 + 30919x<sup>6</sup>/8640 + 5257x<sup>7</sup>/17280 - 8183x<sup>8</sup>/518400 ...</td></tr></tbody></table>

Newton-Cotes formula of degree $n \neq 0$ is:

$$ \int_{a}^{b} f(x) dx \approx \sum_{i=0}^{n} w_{i} f(x_{i}) \qquad (B3) $$
<table><caption>**Table 4:** Local corresponds to the simple rule [29] whereas global corresponds to the composite rule. $\xi$ belongs to the $[a, a + nh]$ interval in the local case and $[a, b]$ in the global case. Function derivatives of order $k$ are given by $f^{(k)}(\xi)$. Adapted from Ref. [28]</caption><thead><tr><th>n</th><th>local interval</th><th>Local error</th><th>Global error</th><th>Rule</th></tr></thead><tbody><tr><td>0</td><td>[a, a + h]</td><td>h<sup>2</sup> <sup>1</sup>&frasl;<sub>2</sub> f<sup>(1)</sup>(ξ)</td><td>(b - a)h<sup>1</sup>&frasl;<sub>2</sub> f<sup>(1)</sup>(ξ) ≈ <sup>1</sup>&frasl;<sub>N</sub></td><td>Rectangular/Euler</td></tr><tr><td>1</td><td>[a, a + h]</td><td>h<sup>3</sup> <sup>1</sup>&frasl;<sub>12</sub> f<sup>(2)</sup>(ξ)</td><td>(b - a)h<sup>2</sup> <sup>1</sup>&frasl;<sub>12</sub> f<sup>(2)</sup>(ξ) ≈ <sup>1</sup>&frasl;<sub>N<sup>2</sup></sub></td><td>Trapezoidal</td></tr><tr><td>2</td><td>[a, a + 2h]</td><td>h<sup>5</sup> <sup>1</sup>&frasl;<sub>90</sub> f<sup>(4)</sup>(ξ)</td><td>(b - a)h<sup>4</sup> <sup>1</sup>&frasl;<sub>2 90</sub> f<sup>(4)</sup>(ξ) ≈ <sup>1</sup>&frasl;<sub>N<sup>4</sup></sub></td><td>1/3 Simpson</td></tr><tr><td>3</td><td>[a, a + 3h]</td><td>h<sup>5</sup> <sup>3</sup>&frasl;<sub>80</sub> f<sup>(4)</sup>(ξ)</td><td>(b - a)h<sup>4</sup> <sup>3</sup>&frasl;<sub>3 80</sub> f<sup>(4)</sup>(ξ) ≈ <sup>1</sup>&frasl;<sub>N<sup>4</sup></sub></td><td>3/8 Simpson</td></tr><tr><td>4</td><td>[a, a + 4h]</td><td>h<sup>7</sup> <sup>8</sup>&frasl;<sub>945</sub> f<sup>(6)</sup>(ξ)</td><td>(b - a)h<sup>6</sup> <sup>1</sup>&frasl;<sub>4 945</sub> f<sup>(6)</sup>(ξ) ≈ <sup>1</sup>&frasl;<sub>N<sup>6</sup></sub></td><td>Boole</td></tr><tr><td>5</td><td>[a, a + 5h]</td><td>h<sup>7</sup> <sup>275</sup>&frasl;<sub>12096</sub> f<sup>(6)</sup>(ξ)</td><td>(b - a)h<sup>6</sup> <sup>1</sup>&frasl;<sub>5 12096</sub> f<sup>(6)</sup>(ξ) ≈ <sup>1</sup>&frasl;<sub>N<sup>6</sup></sub></td><td>Milne</td></tr><tr><td>6</td><td>[a, a + 6h]</td><td>h<sup>9</sup> <sup>9</sup>&frasl;<sub>1400</sub> f<sup>(8)</sup>(ξ)</td><td>(b - a)h<sup>8</sup> <sup>1</sup>&frasl;<sub>6 1400</sub> f<sup>(8)</sup>(ξ) ≈ <sup>1</sup>&frasl;<sub>N<sup>8</sup></sub></td><td>Weddle</td></tr></tbody></table>

where $x_i = i.h + x_0$, with $h$ the step size equal to $(x_n - x_0)/n = (b - a)/n$ and $w_i$ the weights. The weights are obtained typically from Lagrange [28, 29] basis polynomials. They depend only on the $x_i$ and not on function $f$. Lagrange interpolation polynomial corresponding to $n + 1$ points $a \le x_0 < x_1 < ... < x_{n-1} < x_n \le b$ is:

$$L(x) = \sum_{k=0}^{n} f(x_k) \ell_k(x), \quad \ell_k(x) = \prod_{j=0, j \neq k}^{n} \left( \frac{x - x_j}{x_k - x_j} \right).$$

Hence, given $L(x)$ the Lagrange interpolation polynomial [29] for the $n + 1$ data points $(x_0, f(x_0)), \dots, (x_n, f(x_n))$, we get:

$$\int_a^b f(x) dx \approx \int_a^b L(x) dx = \int_a^b \left( \sum_{i=0}^n f(x_i) \ell_i(x) \right) dx \quad (B4)$$

Using finite difference approximation: $f(x) = f(a) + (x - a)f'(\xi)$, the local error is:

$$\epsilon = \int_a^b [f(x) - f(a)] dx = \int_a^b (x - a)f'(\xi) dx = f'(\xi) \int_a^b (x - a) dx \quad (B5)$$

yielding: $\epsilon = \frac{(b-a)^2}{2} f'(\xi) = \frac{h^2}{2} f'(\xi)$ where $\xi$ belongs to the $[a, b]$ integration interval [29]. The global error over the full interval $[a, b]$ is $(b - a)h^{1/2} f'(\xi)$. Table 4 summarizes local and global errors according to Newton-Cotes [28] formulas.

[1] W. Feller, *An Introduction to Probability Theory and its Applications*, Vols. I and II, John Wiley, New York (1968).

[2] A. Papoulis, *Probability Theory, Random Variables and Stochastic Processes*, McGraw Hill (1965).

[3] L. D. Landau and E. M. Lifshitz, *Quantum Mechanics*, Vol. 3, Pergamon, Oxford (1965).

[4] E. Merzbacher, *Quantum Mechanics*, Third edition, Wiley, New-York (1998).

[5] N. Metropolis and S. M. Ulam, *The Monte Carlo Method*, J. Amer. Statist. Assoc., **44**, 333 (1949). J. Von Neumann, *Various techniques used in connection with random digits*, U. S. National Bur. Stand. Appl. Math. Ser., No. 12, page 36, (1951); J. M. Hammersley and D. C. Handscomb, *Monte-Carlo Methods*, Chapman and Hall, London (1964).

[6] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller and E. Teller, *J. Chem. Phys.* **21**, 1087 (1953).

[7] Lord Kelvin, *Nineteenth century clouds over the dynamical theory of heat and light*, Phil. Mag. (6), **2**, 1 (1901).

[8] D. W. Heermann, *Computer Simulation Methods in Theoretical Physics* Second Edition, Springer-Verlag, Berlin (1990). N. G. van Kampen, *Stochastic Processes in Physics and Chemistry*, North Holland (1981). C. W. Gardiner, *A Handbook of Stochastic Methods for Physics, Chemistry and Natural sciences*, Springer (1985). J. B. Kogut, *The lattice gauge-theory approach to quantum chromodynamics*, *Reviews of Modern Physics*, **55**, 775, (1983).

[9] J. R. Rice, *Numerical Methods, Software and Analysis*, Second Edition, Academic Press, New-York (1993).

[10] W. H. Press, W. T. Vetterling, S. A. Teukolsky and B. P. Flannery, *Numerical Recipes in C: The Art of Scientific Computing* Third Edition, Cambridge University Press, New-York (2007).

[11] J.-C. Bajard, D. Michelucci, J.-M. Moreau, J.-M. Muller, Introduction to the Special Issue on "Real Numbers and Computers", *Journal of Universal Computer Science* 1:7, page 438, Jul 28, (1995).

[12] D. P. Landau and K. Binder, *A Guide to Monte-Carlo Simulations in Statistical Physics*, 2nd Edition, Cambridge (2000).
[13] J. E. Gentle, *Random number generation and Monte-Carlo methods*, Second Edition, Springer Science Business Media, Inc., New-York (2003).

[14] L. Devroye, *Non-uniform random variate generation*, Springer-Verlag, New-York (1986).

[15] G. Marsaglia, Proc. Nat. Acad. Sci. **61**, 25 (1968); G. Marsaglia, *Applications of Number Theory to Numerical Analysis*, Academic, New York (1972). G. Marsaglia, *Generating exponential random variables*, Ann. Math. Stat., **32**, 899 (1961); J. H. Ahren, and U. Dieter, *Extensions of Forsythe's method for random sampling from normal distribution*, Math. Comp., **27**, 927 (1973); G. E. Forsythe, *Von Neumann's comparison method for random sampling and from the normal and other distributions*, Math. of Comp. **26**, 817 (1972).

[16] A.M. Ferrenberg, D.P. Landau and Y. J. Wong, *Monte-Carlo simulations: Hidden errors from "Good" Random Number Generators*, Phys. Rev. Lett. **69**, 3382 (1992).

[17] G. L. Buffon, *Editor note concerning a lecture given in 1733 by Mr. Leclerc de Buffon to the Royal Academy of Sciences in Paris*. Histoire de l'Académie. Royale des Sciences, 43-45. G. L. Buffon, *Essai darithmétique morale*, Histoire naturelle, générale er particulière, Supplément 4, 46123. (1777).

[18] Knuth, D.E. 1981, Seminumerical Algorithms, 2nd ed., *The Art of Computer Programming* vol. 2 (Addison-Wesley).

[19] G. F. Arfken and H. J. Weber, *Mathematical methods for physicists*, Sixth Edition, Elsevier Academic Press (London, 2005).

[20] M. N. O. Sadiku, *Monte-Carlo methods for Electromagnetics*, CRC Press, Taylor & Francis Group Boca Raton, Florida (2009).

[21] A. George and J. Liu, *Computer solution of large sparse positive definite systems*. Prentice Hall, Englewood Cliffs, NJ, (1981).

[22] G. E. Forsythe and R. A. Leibler, *Matrix inversion by a Monte-Carlo method*, Mathematical Tables and other Aids to Computation, 4(31):127-127, (1950).

[23] G. Oktent, *Solving Linear equations by Monte-Carlo simulation*, SIAM J. Sci. Comput. **27**, 511 (2005).

[24] T. P. Plaks, G. M. Megson, and V. N. Alexandrov, *A family of efficient regular arrays for band and sparse matrix inversion using Monte-Carlo methods*, in Proceedings of the 16th IMACS World Congress 2000 on Scientific Computation, Applied Mathematics and Simulation, M. Deville and R. Owens, editors, Lausanne, Switzerland (2000).

[25] N. J. Giordano, *Computational physics*, Prentice-Hall, Upper Saddle River, New Jersey [1997].

[26] L. D. Landau and E. M. Lifshitz, *Statistical Physics*, Vol. 5, Pergamon, Oxford (1975).

[27] D. A. McQuarrie, *Statistical Mechanics*, Harper and Row, New York (1973).

[28] J. Stoer and R. Bulirsch, *Introduction to Numerical Analysis*, Second Edition, Springer-Verlag, (New York, 1992).

[29] D.E. Potter, *Computational physics*, John Wiley & Sons, London, New York [1973].