> Here's my implementation of Simpson's, except it divides the interval into
> 2n segments.
>
> >>> def simpson(f,a,b,n):
> h = float(b-a)/(2*n)
> sum1 = sum(f(a + 2*k *h) for k in range(1,n))
> sum2 = sum(f(a + (2*k-1)*h) for k in range(1,n+1))
> return (h/3)*(f(a)+f(b)) + (2*h/3)*sum1 + (4*h/3)*sum2
>
> >>> def g(x): return x*x
>
> >>> simpson(g,0,3,10000)
> 8.9999999999999893
> >>> simpson(g,0,3,5000)
> 8.9999999999999787
>
> Kirby
And with a decorator:
def simpson(f):
def defint(a,b,n=1000):
h = float(b-a)/(2*n)
sum1 = sum(f(a + 2*k *h) for k in range(1,n))
sum2 = sum(f(a + (2*k-1)*h) for k in range(1,n+1))
return (h/3)*(f(a)+f(b)) + (2*h/3)*sum1 + (4*h/3)*sum2
return defint
@simpson
def g(x): return x*x
>>> g(0, 3)
9.0000000000000036
>>> g(0, 4)
21.333333333333329
>>> g(1, 7)
113.99999999999983
>>> g(1, 7, n=10000)
114.0
Kirby
_______________________________________________
Edu-sig mailing list
Edu-sig@python.org
http://mail.python.org/mailman/listinfo/edu-sig
