Monte Carlo Method of Computing Pi – Numerical Explorations

Triple Integration of Three Functions Using a Monte Carlo Algorithm and an Adaptive Quadrature Method (c) January 15 – 17, January 2024 by James Pate Williams, Jr.

Back in 2015 I translated a multiple integration FORTRAN subroutine to C# using switch statements to emulate conditional and unconditional go to statements. On January 15 – January 16, 2024, I translated my C# source code to vanilla C. Below is the FORTRAN subroutine’s website:

Remarks on algorithm 006: An adaptive algorithm for numerical integration over an N-dimensional rectangular region – ScienceDirect

Here are some web pages with examples of triple integration:

3.4: Numerical Approximation of Multiple Integrals – Mathematics LibreTexts

3.3: Triple Integrals – Mathematics LibreTexts

https://tutorial.math.lamar.edu/Classes/CalcIII/IteratedIntegrals.aspx

https://tutorial.math.lamar.edu/Solutions/CalcIII/TripleIntegrals/Prob1.aspx

https://tutorial.math.lamar.edu/Problems/CalcIII/TripleIntegrals.aspx

I created a C test application to attempt to verify some of the results in online sources. For one example I also used my C# program.

== Menu ==

1 f(x, y, z) = 8 * x * y * z

2 f(x, y, z) = 4 * x * x * y - z * z * z

3 f(x, y, z) = x * y + z

5 Exit

Option: 1

f(x, y, z) = 8 * x * y * z

[0.000000, 1.000000] x [1.000000, 2.000000] x [2.000000, 3.000000]

N       Integral        -+Error         % Error

10      0.754706        0.250631        94.968627

100     1.050263        0.126084        92.998246

1000    1.032125        0.036751        93.119165

10000   1.015068        0.011748        93.232880

100000  1.004080        0.003706        93.306131

1000000 1.000650        0.001169        93.329000

nQuadrature Integral Value and Percent Error: 15.000000 0.000000

== Menu ==

1 f(x, y, z) = 8 * x * y * z

2 f(x, y, z) = 4 * x * x * y - z * z * z

3 f(x, y, z) = x * y + z

5 Exit

Option:

== Menu ==

1 f(x, y, z) = 8 * x * y * z

2 f(x, y, z) = 4 * x * x * y - z * z * z

3 f(x, y, z) = x * y + z

5 Exit

Option: 2

f(x, y, z) = 4 * x * x * y - z * z * z

[2.000000, 3.000000] x [-1.000000, 4.000000] x [1.000000, 0.000000]

N       Integral        -+Error         % Error

10      -12.087502      -3.929801       93.596026

100     -17.622371      -1.969516       90.663644

1000    -17.963837      -0.612528       90.482735

10000   -18.129559      -0.197880       90.394936

100000  -18.018231      -0.062690       90.453917

1000000 -17.939045      -0.019779       90.495870

nQuadrature Integral Value and Percent Error: 188.750000        200.000000

== Menu ==

1 f(x, y, z) = 8 * x * y * z

2 f(x, y, z) = 4 * x * x * y - z * z * z

3 f(x, y, z) = x * y + z

5 Exit

Option:

== Menu ==

1 f(x, y, z) = 8 * x * y * z

2 f(x, y, z) = 4 * x * x * y - z * z * z

3 f(x, y, z) = x * y + z

5 Exit

Option: 3

f(x, y, z) = x * y + z

[0.000000, 3.000000] x [0.000000, 2.000000] x [0.000000, 1.000000]

N       Integral        -+Error         % Error

10      9.243858        2.121562        22.967849

100     12.062489       0.808150        0.520745

1000    12.096541       0.255335        0.804509

10000   12.094473       0.081030        0.787275

100000  12.051033       0.025704        0.425271

1000000 12.009588       0.008121        0.079900

nQuadrature Integral Value and Percent Error: 12.000000 0.000000

== Menu ==

1 f(x, y, z) = 8 * x * y * z

2 f(x, y, z) = 4 * x * x * y - z * z * z

3 f(x, y, z) = x * y + z

5 Exit

Option:

desired relative error 0.001



f#	integral		epsilon		number	err code

1	+1.4346639496E+000	1.295116E-003	257	0

2	+5.7531639665E-001	5.745470E-004	97	0

3	+2.1527578485E+000	1.793429E-003	45	0

4	+1.5998921741E+001	1.579802E-002	151	0

5	+1.8390615688E-001	6.510967E-005	97	0

6	-4.0003324629E+000	3.747971E-003	17	0

7	+8.6330831791E-001	8.563683E-004	45	0

8	+1.5000000000E+001	8.850520E-011	45	0

9	+1.8875000000E+002	1.109797E-009	45	0

10	+1.2000000000E+001	7.081136E-011	45	0



f#	abs errror		percent error

1	+9.7938837210E-005	+6.8261387483E-003

2	+4.7748257754E-005	+8.2987892410E-003

3	+6.1501589893E-004	+2.8576909005E-002

4	+1.0782593204E-003	+6.7391207526E-003

5	+9.9602348685E-007	+5.4159040098E-004

6	+3.3246288154E-004	+8.3115720386E-003

7	+2.6210055455E-004	+3.0369237392E-002

8	+8.8506979523E-011	+5.9004653015E-010

9	+1.1098109098E-009	+5.8797928998E-010

10	+7.0812689046E-011	+5.9010574205E-010

Numbers 8-10 in the preceding data correspond to the three three-dimensional functions in the menu illustrated above. The website claims option 2 in the Menu integral is -755 / 4 = -188.75. My calculation is the same magnitude but a positive sign.

Double Integration of Six Functions Using Three Different Quadrature Algorithms (c) January 15, 2024, by James Pate Williams, Jr. All Applicable Rights Reserved

I utilized six functions that have two independent variables. The integration methods were a Monte Carlo probabilistic procedure, TRICUB from “A Numerical Library in C for Scientists and Engineers” by H. T. Lau, PhD, and an algorithm from the ACM publication. I translated the ACM method from FORTRAN to C. See the following online references:

https://math.libretexts.org/Bookshelves/Calculus/Vector_Calculus_(Corral)/03%3A_Multiple_Integrals/3.04%3A_Numerical_Approximation_of_Multiple_Integrals
https://tutorial.math.lamar.edu/Classes/CalcIII/IteratedIntegrals.aspx

== Menu ==

1 8 * x + 6 * y

2 6 * x * y * y

3 2 * x - 4 * y * y * y

4 x * x * y * y + cos(pi * x) + sin(pi * y)

5 pow(2 * x + 3 * y, -2)

6 x * exp(x * y)

7 Exit

Option: 1

f(x, y) = 8x + 6y

[0.000000, 1.000000] x [0.000000, 2.000000]

N       Integral        +-Error         % Error

10      18.398218       2.687822        8.008911

100     19.546981       0.838888        2.265096

1000    19.859442       0.266101        0.702789

10000   20.152819       0.082772        0.764093

100000  20.014598       0.026273        0.072988

1000000 20.013430       0.008323        0.067150

NUMAL Integral Value and Percent Error: 22.666667       13.333333

nQuadrature Integral Value and Percent Error: 19.999968 0.000160

== Menu ==

1 8 * x + 6 * y

2 6 * x * y * y

3 2 * x - 4 * y * y * y

4 x * x * y * y + cos(pi * x) + sin(pi * y)

5 pow(2 * x + 3 * y, -2)

6 x * exp(x * y)

7 Exit

Option:

== Menu ==

1 8 * x + 6 * y

2 6 * x * y * y

3 2 * x - 4 * y * y * y

4 x * x * y * y + cos(pi * x) + sin(pi * y)

5 pow(2 * x + 3 * y, -2)

6 x * exp(x * y)

7 Exit

Option: 2

f(x, y) = 6 * x * y * y

[1.000000, 2.000000] x [2.000000, 4.000000]

N       Integral        +-Error         % Error

10      6.400786        2.133333        92.380017

100     7.528110        0.951714        91.037965

1000    7.951653        0.294436        90.533746

10000   8.083377        0.094959        90.376933

100000  8.001852        0.029859        90.473986

1000000 8.010457        0.009474        90.463742

NUMAL Integral Value and Percent Error: 204.800000      143.809524

nQuadrature Integral Value and Percent Error: 167.999730        99.999679

== Menu ==

1 8 * x + 6 * y

2 6 * x * y * y

3 2 * x - 4 * y * y * y

4 x * x * y * y + cos(pi * x) + sin(pi * y)

5 pow(2 * x + 3 * y, -2)

6 x * exp(x * y)

7 Exit

Option:

== Menu ==

1 8 * x + 6 * y

2 6 * x * y * y

3 2 * x - 4 * y * y * y

4 x * x * y * y + cos(pi * x) + sin(pi * y)

5 pow(2 * x + 3 * y, -2)

6 x * exp(x * y)

7 Exit

Option: 3

f(x, y) = 2 * x - 4 * y * y * y

[-5.000000, 4.000000] x [0.000000, 3.000000]

N       Integral        +-Error         % Error

10      -645.214459     239.379239      14.654172

100     -436.170436     84.138767       42.305498

1000    -471.478924     25.861533       37.635063

10000   -486.066536     8.338088        35.705485

100000  -484.301805     2.644031        35.938915

1000000 -486.527596     0.839338        35.644498

NUMAL Integral Value and Percent Error: -765.000000     1.190476

nQuadrature Integral Value and Percent Error: -755.998734       0.000167

== Menu ==

1 8 * x + 6 * y

2 6 * x * y * y

3 2 * x - 4 * y * y * y

4 x * x * y * y + cos(pi * x) + sin(pi * y)

5 pow(2 * x + 3 * y, -2)

6 x * exp(x * y)

7 Exit

Option:

== Menu ==

1 8 * x + 6 * y

2 6 * x * y * y

3 2 * x - 4 * y * y * y

4 x * x * y * y + cos(pi * x) + sin(pi * y)

5 pow(2 * x + 3 * y, -2)

6 x * exp(x * y)

7 Exit

Option: 4

f(x, y) = x * x * y * y + cos(pi * x) + sin(pi * y)

[-2.000000, -1.000000] x [0.000000, 1.000000]

N       Integral        +-Error         % Error

10      1.103184        0.173977        22.003269

100     0.708005        0.063869        49.942995

1000    0.757961        0.021546        46.411062

10000   0.744889        0.006866        47.335224

100000  0.748463        0.002171        47.082576

1000000 0.746892        0.000686        47.193609

NUMAL Integral Value and Percent Error: 1.449491        2.481177

nQuadrature Integral Value and Percent Error: 1.414395  0.000160

== Menu ==

1 8 * x + 6 * y

2 6 * x * y * y

3 2 * x - 4 * y * y * y

4 x * x * y * y + cos(pi * x) + sin(pi * y)

5 pow(2 * x + 3 * y, -2)

6 x * exp(x * y)

7 Exit

Option:

== Menu ==

1 8 * x + 6 * y

2 6 * x * y * y

3 2 * x - 4 * y * y * y

4 x * x * y * y + cos(pi * x) + sin(pi * y)

5 pow(2 * x + 3 * y, -2)

6 x * exp(x * y)

7 Exit

Option: 5

f(x, y) = pow(2 * x + 3 * y, -2)

[0.000000, 1.000000] x [1.000000, 2.000000]

N       Integral        +-Error         % Error

10      0.555876        0.235032        1394.669315

100     0.752953        0.323309        1924.579264

1000    0.826769        0.275474        2123.058629

10000   0.927229        0.177631        2393.180482

100000  0.995193        0.097358        2575.927399

1000000 1.291253        0.127325        3371.987619

NUMAL Integral Value and Percent Error: 0.034125        8.243058

nQuadrature Integral Value and Percent Error: 0.037190  0.000639

== Menu ==

1 8 * x + 6 * y

2 6 * x * y * y

3 2 * x - 4 * y * y * y

4 x * x * y * y + cos(pi * x) + sin(pi * y)

5 pow(2 * x + 3 * y, -2)

6 x * exp(x * y)

7 Exit

Option:

== Menu ==

1 8 * x + 6 * y

2 6 * x * y * y

3 2 * x - 4 * y * y * y

4 x * x * y * y + cos(pi * x) + sin(pi * y)

5 pow(2 * x + 3 * y, -2)

6 x * exp(x * y)

7 Exit

Option: 6

f(x, y) = x * exp(x * y)

[0.000000, 1.000000] x [-1.000000, 2.000000]

N       Integral        +-Error         % Error

10      3.217934        1.437426        inf

100     5.562365        0.873166        inf

1000    5.325378        0.225123        inf

10000   5.497069        0.075770        inf

100000  5.362919        0.023344        inf

1000000 5.371802        0.007410        inf

NUMAL Integral Value and Percent Error: 6.162907        inf

nQuadrature Integral Value and Percent Error: 2.562339  inf

Unknown option, please try again.== Menu ==

1 8 * x + 6 * y

2 6 * x * y * y

3 2 * x - 4 * y * y * y

4 x * x * y * y + cos(pi * x) + sin(pi * y)

5 pow(2 * x + 3 * y, -2)

6 x * exp(x * y)

7 Exit

Option:

Three Ways of Computing a Few Digits of Pi by James Pate Williams, Jr.

I wrote a short C++ program to calculate a few digits of pi, a famous transcendental number. The algorithms are as follows:

Monte Carlo Method
Leibniz’s Infinite Series
Nilakantha’s Infinite Series

First the results then the C++ source code listing.

Calculating a Few Digits of the Transcendental Number Pi by Throwing Darts by James Pate Williams, BA, BS, MSwE, PhD

Suppose you have a unit square with a circle of unit diameter inscribed . You can compute a few digits of the transcendental number, pi, 3.1415926535897932384626433832795…, by using the algorithm described as follows. Let n be the number of darts to throw and h be the number of darts that land within the inscribed circle.

h = 0

for i = 1 to n do

Choose two random numbers x and y such that x and y are contained in the interval 0 to 1 inclusive that is x and y contained in [0, 1]

Let u = x – 1 / 2 and v = y – 1 / 2

if u * u + v * v <= 0.25 = 1 / 4 then h = h + 1

next i

pi = 4 * h / n

Below are the results of a C# Microsoft Visual Studio simulation project. In the first case we throw 100,000 darts and get two significant digits of pi and then we throw a 1,000,000 darts and five significant digits of pi are computed. Of course, in a previous entry by this author we can calculate hundreds or thousands of digits of pi in relatively little time:

https://jamespatewilliamsjr.wordpress.com/2018/07/01/the-bailey-borwein-plouffe-formula-for-calculating-the-first-n-digits-of-pi/

MainForm 10_16_2018 3_04_04 AM MainForm 10_16_2018 3_05_45 AM MainForm 10_16_2018 3_09_36 AM MainForm 10_16_2018 3_09_54 AM