[Edu-sig] Re: Integration correction

[André Roberge]

Kirby Urner wrote:
>>From: Kirby Urner <urnerk at qwest.net>
>[snip   ... Kirby, Art and many others including myself discuss
>  the possible misuse of decorators in the context of calculating
 > derivatives and integrals numerically .... end snip]
> 
> Now, if g(x) really *did* go on for 30-40 lines, OK, then maybe a decorator
> adds to readability.
> 
> Something to think about.
> 
> Kirby

Well, I don't think that decorators would add to readability once you've 
defined your functions.  I've played around with using decorators (and 
not) for a while, and here's what I would like to submit to this list
as a possible "summary".  This non-decorator approach is very clean imho.


def derivative(f):
     """
     Computes the numerical derivative of a function.
     Intended to be used as a decorator.
     """
     def df(x, h=0.1e-5):
         return ( f(x+h/2) - f(x-h/2) )/h
     return df

def simpson_integral(f, a, b, n=1000):
     """
     Computes the numerical integral of a function
     using the extended Simpson's rule.  This gives exact results
     for polynomials of O(3).
     """
     h = float(b-a)/n
     sum1 = sum(f(a+k*h) for k in range(1, n, 2))
     sum2 = sum(f(a+k*h) for k in range(2, n, 2))
     return (h/3)*( f(a)+f(b) + 4*sum1 + 2*sum2 )

# test functions
def g2(x): return x*x
def g3(x): return x*x*x

# first derivatives
dg2 = derivative(g2)
dg3 = derivative(g3)
# second derivatives
d2g2 = derivative(derivative(g2))
d2g3 = derivative(derivative(g3))
# integrals
intg2 = simpson_integral(g2, 0, 1)
intg3 = simpson_integral(g3, 0, 1)


print '         functions: ', g2(2.), g3(2.)
print ' first derivatives: ', dg2(2.), dg3(2.)
print 'second derivatives: ', d2g2(2.), d2g3(2.)
print 'numerical integral: ', intg2, intg3

#         functions:  4.0 8.0
# first derivatives:  4.00000000012 12.0000000008
#second derivatives:  2.00017780116 12.001066807
#numerical integral:  0.333333333333 0.25



#== WARNING ==
#-- derivatives computed numerically do not
#-- necessarily converge for small h --

h = .1
for i in range(1, 7):
     h /= 100
     print h, 'first = ', dg2(2, h), 'second =', dg3(2, h)

#0.001 first =  4.0 second = 12.00000025
#1e-005 first =  3.99999999989 second = 11.9999999997
#1e-007 first =  4.00000000234 second = 12.0000000159
#1e-009 first =  4.00000033096 second = 12.0000009929
#1e-011 first =  4.00000033096 second = 12.0000009929
#1e-013 first =  4.00568467285 second = 12.0170540185

#####################

André