/usr/share/pyshared/Scientific/Functions/Romberg.py is in python-scientific 2.8-4.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 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 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 | # Romberg quadratures for numeric integration.
#
# Written by Scott M. Ransom <ransom@cfa.harvard.edu>
# last revision: 14 Nov 98
#
# Cosmetic changes by Konrad Hinsen <hinsen@cnrs-orleans.fr>
# last revision: 2006-11-23
#
"""
Numerical integration using the Romberg algorithm
"""
def trapezoid(function, interval, numtraps):
"""
Numerical X{integration} using the X{trapezoidal} rule
Example::
>>>from Scientific.Functions.Romberg import trapezoid
>>>from Scientific.N import pi, tan
>>>trapezoid(tan, (0.0, pi/3.0), 100)
yields 0.69317459482518262
@param function: a function of one variable
@type function: callable
@param interval: the lower and upper limit of the integration interval
@type interval: sequence of two C{float}s
@param numtraps: the number of trapezoids
@type numtraps: C{int}
@returns: the numerical integral of the function over the interval
@rtype: number
"""
lox, hix = interval
h = float(hix-lox)/numtraps
sum = 0.5*(function(lox)+function(hix))
for i in range(1, numtraps):
sum = sum + function(lox + i*h)
return h*sum
def _difftrap(function, interval, numtraps):
# Perform part of the trapezoidal rule to integrate a function.
# Assume that we had called _difftrap with all lower powers-of-2
# starting with 1. Calling _difftrap only returns the summation
# of the new ordinates. It does _not_ multiply by the width
# of the trapezoids. This must be performed by the caller.
# 'function' is the function to evaluate.
# 'interval' is a sequence with lower and upper limits
# of integration.
# 'numtraps' is the number of trapezoids to use (must be a
# power-of-2).
if numtraps<=0:
print "numtraps must be > 0 in _difftrap()."
return
elif numtraps==1:
return 0.5*(function(interval[0])+function(interval[1]))
else:
numtosum = numtraps/2
h = float(interval[1]-interval[0])/numtosum
lox = interval[0] + 0.5 * h;
sum = 0.0
for i in range(0, numtosum):
sum = sum + function(lox + i*h)
return sum
def _romberg_diff(b, c, k):
# Compute the differences for the Romberg quadrature corrections.
# See Forman Acton's "Real Computing Made Real," p 143.
tmp = 4.0**k
return (tmp * c - b)/(tmp - 1.0)
def romberg(function, interval, accuracy=1.0E-7):
"""
Numerical X{integration} using the X{Romberg} method
Example::
>>>from Scientific.Functions.Romberg import romberg
>>>from Scientific.N import pi, tan
>>>romberg(tan, (0.0, pi/3.0))
yields '0.693147180562'
@param function: a function of one variable
@type function: callable
@param interval: the lower and upper limit of the integration interval
@type interval: sequence of two C{float}s
@param accuracy: convergence criterion (absolute error)
@type accuracy: C{float}
@returns: the numerical integral of the function over the interval
@rtype: number
"""
i = n = 1
intrange = interval[1] - interval[0]
ordsum = _difftrap(function, interval, n)
result = intrange * ordsum
resmat = [[result]]
lastresult = result + accuracy * 2.0
while (abs(result - lastresult) > accuracy):
n = n * 2
ordsum = ordsum + _difftrap(function, interval, n)
resmat.append([])
resmat[i].append(intrange * ordsum / n)
for k in range(i):
resmat[i].append(_romberg_diff(resmat[i-1][k],
resmat[i][k], k+1))
result = resmat[i][i]
lastresult = resmat[i-1][i-1]
i = i + 1
return result
# Test code
if __name__ == '__main__':
from math import tan, pi
print ''
r = romberg(tan, (0.5, 1.0))
t = trapezoid(tan, (0.5, 1.0), 1000)
theo = 0.485042229942291
print ''
print 'Trapezoidal rule with 1000 function evaluations =', t
print 'Theoretical result (correct to all shown digits) = %15.14f' % theo
print ''
|