Welcome, guest | Sign In | My Account | Store | Cart
```# 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)
```

#### Diff to Previous Revision

```--- revision 1 2010-06-14 00:05:14
+++ revision 2 2010-06-16 19:31:17
@@ -29,8 +29,10 @@
x = xmin + (xmax - xmin) * random.random()
y = ymin + (ymax - ymin) * random.random()
if math.fabs(y) <= math.fabs(f(x)):
-        # area over x-axis is positive, and under is negative
-        ctr += math.copysign(1, y)
+    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)
```