numpy - Integrate with Simpson's Rule in Python recursively -

i want use python implement simpson integration. it's not hard if don't need converge (here need abs(result - my_expect) < 0.001) automatically. want automatic-converged program written in python. tried use method learned sicp -- make recursively.

# imports __future__ import division import numpy np  # constants mu = 50 sigma = 15 x0 = 0 xn = 100   # interface functions def f(x):     """integrand.      used lagrangian interpolation ,     calculating corresponding function value in function     `improve_precision`.      """     y = (1 / (sigma * np.sqrt(2 * np.pi))          * np.exp(-1/2 * (x-mu)**2 / sigma**2))     return y   # classes class integration:     def __init__(self, start, end):         self.start = start         self.end = end         self.a = (self.end - self.start) / 2         self.h = 2 * self.a         self.rp = 0         self.s = 0         self.i = 1      def integrate(self):         return self.integrate_iter(1, 0, 0.001)      def integrate_iter(self, expectation, summation, accuracy):         if accuracy > abs(summation - expectation):             return summation         else:             self.integrate_iter(expectation,                                 self.improve_precision(self.i+1),                                  accuracy)      def improve_precision(self, ii):         rc = np.sum([f(-self.a - self.h/2 + k*self.h)                      k in range(1, 2**ii + 1)])         self.s = self.h/6 * (self.rp + 4*rc)         self.h /= 2         self.rp = self.rp + 2*rc         return self.s  calc = integration(x0, xn) print calc.integrate() 

when ran it, lot of errors raised:

  file "/users/*/integrate.py", line 72, in <module>     print calc.integrate()   file "/users/*/integrate.py", line 56, in integrate     return self.integrate_iter(1, 0, 0.001)   file "/users/*/integrate.py", line 62, in integrate_iter     self.integrate_iter(expectation, self.improve_precision(self.i+1), accuracy)   file "/users/*/integrate.py", line 62, in integrate_iter     self.integrate_iter(expectation, self.improve_precision(self.i+1), accuracy)   ...   file "/users/*/integrate.py", line 65, in improve_precision     rc = np.sum([f(-self.a - self.h/2 + k*self.h) k in range(1, 2**ii + 1)])   file "/usr/local/lib/python2.7/site-packages/numpy/core/fromnumeric.py", line 1708, in sum     if isinstance(a, _gentype): runtimeerror: maximum recursion depth exceeded while calling python object process finished exit code 1 

so should do? there way make correct?

edit:

well, made change on function integrate_iter , added function is_rough. final code goes this:

...     def integrate(self):         """a function encapsulate calculating process."""         integrate_value = self.integrate_iter(0)         interval_count = self.i         return integrate_value, interval_count      def integrate_iter(self, summation):         """main part of recursive process."""         last_summation = summation         summation = self.improve_precision(self.i)         return self.is_good(summation, last_summation)      def is_good(self, summ, last_summ):         # if precision < 0.001, over; else, call `integrate_iter` again.         if abs(summ - last_summ) <= self.accuracy:             return summ         else:             self.i += 1             return self.integrate_iter(summ) ... calc = integration(x0, xn, 0.001) val = calc.integrate() scipy_result = scipy.integrate.quad(f, x0, xn, epsabs=0.001)[0] print 'the integration is: {0:.8f}'.format(val[0]) print 'the error is: {0:.8f}'.format(abs(val[0] - scipy_result)) print 'the number of recursions is: {0:d}'.format(val[1]) print 'the number of intervals is: {0:d}'.format(2**val[1]) 

i'm not quite satisfied number of recursion, don't know whether should converge slow(according experience).

two fixes needed inside integrate_iter.

first, when increment i, need not add 1 i, redefine i i + 1.

second, in posted code, integrate_iter returns none when else block executed. prevents improved result being obtained within integrate, needed can returned in print statement.

def integrate_iter(self, expectation, summation, accuracy):     if accuracy > abs(summation - expectation):         return summation     else:         self.i += 1         return self.integrate_iter(expectation,                                    self.improve_precision(self.i),                                     accuracy) 

although above works, think it's clearer write way:

def integrate_iter(self, expectation, summation, accuracy):     if accuracy <= abs(summation - expectation):         self.i += 1         # improve summation.         summation = self.integrate_iter(expectation,                                         self.improve_precision(self.i),                                          accuracy)     # summation enough; return.     return summation