# "Numerical Recipes in C", 2nd ed, trapzd() on page 137
# We need to pass a function with one argument to it.
# We define a sample square function first.
import math
# The square function
def FUNC(x):
return (x*x)
# trapzd takes a function as input, with a limit (a,b), at n-th iteration.
# C uses static variable s. But in Python, I better pass it as argument.
# the last argument s is the value calculated in previous iteration.
def trapzd(FUNC, a, b, n, s):
if n==1:
s = 0.5*(b-a)*(FUNC(a)+FUNC(b))
return s
else:
it = 2**(n-2)
d = (b-a)/it
x = a + 0.5*d
sum = 0.0
for j in range(0,it):
sum += FUNC(x)
x += d
s = 0.5*(s + (b-a)*sum/it)
return s
# end if
# end def trapzd
# Call trapzd n time to get to the desired accuracy
# qtrap() returns the integral of the function FUNC from a to b.
# The parameter EPS can be set to the desired relative accuracy
# JMAX is the maximum allowed number of steps, refining the interval
# to (b-a)/2**(JMAX-1). Integration is performed by the trapzoidal rule.
def qtrap(FUNC, a, b):
EPS = 1.0e-8
JMAX = 20
olds = -1.0e300
s = 0.0
for j in range(1, JMAX):
s = trapzd(FUNC,a,b,j,s)
if (math.fabs(s-olds) < EPS*math.fabs(olds)):
return s
if(s == 0.0 and olds == 0.0 and j > 6):
return s
olds = s
# end def qtrap
# testing int_0^1 x^2 dx:
res = qtrap(FUNC,0.0,1.0)
print(res)