递归地与 Python 中的辛普森规则集成

Question

我想使用 Python 来实现 Simpson 集成。如果我不需要它自动收敛（这里我需要abs(result - my_expect) < 0.001）并不难。 但我想要一个用 Python 编写的自动融合程序。所以我尝试使用我从SICP中学到的方法——递归实现。

# imports
from __future__ import division
import numpy as np

# constants
MU = 50
SIGMA = 15
X0 = 0
XN = 100


# interface functions
def f(x):
    """Integrand.

    It will be used for Lagrangian interpolation and
    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)
                     for 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()

当我运行它时，出现了很多错误：

  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) for 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 a Python object
Process finished with exit code 1

那我该怎么办？有什么办法可以改正吗？

编辑：

好吧，我对函数integrate_iter 进行了更改，并添加了一个函数is_rough。我的最终代码是这样的：

...
    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])

我对递归的次数不是很满意，不知道是不是应该收敛这么慢（根据我的经验）。

Answer 1

integrate_iter 内部需要两个修复。

首先，当您增加i 时，您不仅需要将1 加到i，还需要将i 重新定义为i + 1。

其次，在您发布的代码中，integrate_iter 在执行else 块时返回None。这可以防止在integrate 中获得改进的结果，因为它是需要的，以便可以在print 语句中返回。

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)

虽然上述方法可行，但我认为这样写更清楚：

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 is good enough; return.
    return summation

Answer 2

你错过了堆栈跟踪中的主线：

RuntimeError: maximum recursion depth exceeded while calling a Python object

您的递归 integrate_iter 不会终止。

Answer 3

除非这是一个学习练习，否则您应该改用scipy.integrate.quad --- 它以更高级的方式进行自适应集成。