Numerical Integration using Monte Carlo method.
Python, 38 lines
Download
Copy to clipboard1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 | # Numerical Integration using Monte Carlo method
# FB - 201006137
import math
import random
# define any function here!
def f(x):
return math.sin(x)
# define any xmin-xmax interval here! (xmin < xmax)
xmin = 0.0
xmax = 2.0 * math.pi
# find ymin-ymax
numSteps = 1000000 # bigger the better but slower!
ymin = f(xmin)
ymax = ymin
for i in range(numSteps):
x = xmin + (xmax - xmin) * float(i) / numSteps
y = f(x)
if y < ymin: ymin = y
if y > ymax: ymax = y
# Monte Carlo
rectArea = (xmax - xmin) * (ymax - ymin)
numPoints = 1000000 # bigger the better but slower!
ctr = 0
for j in range(numPoints):
x = xmin + (xmax - xmin) * random.random()
y = ymin + (ymax - ymin) * random.random()
if math.fabs(y) <= math.fabs(f(x)):
if f(x) > 0 and y > 0 and y <= f(x):
ctr += 1 # area over x-axis is positive
if f(x) < 0 and y < 0 and y >= f(x):
ctr -= 1 # area under x-axis is negative
fnArea = rectArea * float(ctr) / numPoints
print "Numerical integration = " + str(fnArea)
|
I certainly did not write this code for any practical purpose. Monte Carlo is probably the most inefficient solution method for almost any problem. I only put it here thinking some people may find it interesting, that's all.
For higher-dimensional integrals, Monte Carlo is often the tool of choice. Yes, it's inefficient for single integrals, but it's a great thing for students to look at because a) it's simple to understand (no need of calculus) and b) it's easy to code.
If you're not using python 3, you should get in the habit of using xrange instead of range in your for loops. The range() built-in creates a large list of numbers, whereas xrange uses lazy evaluation. Otherwise, some programs may experience an out of memory condition (yes, it does happen when dealing with large arrays, even on a Linux box with 4 GB of memory).
thanks for this post!