Simpson方法自适应实现epsilon=1e8，只完成一次迭代

Question

我试图在逻辑中实现这个问题，用epsilon=1e-8计算Simpson的方法自适应实现。解释如下：

"""The approximate definite integral of function from a to b using Simpson's method. 
This function is vectorized, it uses numpy array operations to calculate the approximation.
This is an adaptive implementation, the method starts out with N=2 intervals, and try
successive sizes of N (by doubling the size), until the desired precision, is reached.
This adaptive solution uses our improved approach/equation for Simpson's method, to avoid 
unnecessary recalculations of the integrand function.

a, b - Scalar float values, the begin, and endpoints of the interval we are to
        integrate the function over.
f - A vectorized function, should accept a numpy array of x values, and compute the
      corresponding y values for all points according to some function.
epsilon - The desired precision to calculate the integral to.  Default is 8 decimal places
      of precision (1e-8)

returns - A tuple, (ival, error).  A scalar float value, the approximated integral of 
       the function over the given interval, and a scaler float value of the 
       approximation error on the integral"""

以下是我对辛普森方法的代码：

import pylab as pl
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
%matplotlib inline
def simpsons_adaptive_approximation(a, b, f, epsilon=1e-8):

    N_prev = 2 # the previous number of slices
    h_prev = (b - a) / N_prev # previous interval width
    x = np.arange(a+h_prev, b, h_prev) # x locations of the previous interval
    I_prev =  h_prev * (0.5 * f(a) + 0.5 * f(b) + np.sum(f(x)))

    # set up variables to adaptively iterate successively better approximations
    N_cur  = 4 # the current number of slices
    I_cur  = 0.0 # calculated in loop iteration
    error  = 1.0 # calculated in loop iteration
    itr    = 1 # keep track of the number of iterations we perform, for display/debug

    h = (b-a)/float(epsilon)
    I_cur = f(a) + f(b)

    for i in frange(1,epsilon,1):
        if(i%2 ==0):
            I_cur = I_cur + (2*(f(a + i*h)))
        else:
            I_cur = I_cur + (4*(f(a + i*h)))
    error = np.abs((1.0/3.0) * (I_cur - I_prev))
    print("At iteration %d (N=%d), val=%0.16f prev=%0.16f error=%e" % (itr, N_cur, I_cur, I_prev, error) )

    I_cur *= (h/3.0)

    I_prev = I_cur
    N_prev = N_cur
    N_cur *= 2
    itr += 1

return (I_cur, error)

下面是我的f2（x）函数：

def f2(x):
    return x**4 - 2*x + 1

a = 0.0
b = 2.0
eps = 1e-10
(val, err) = simpsons_adaptive_approximation(a, b, f2, eps)
print( "Calculated value: %0.16f  error: %e  for an epsilon of: %e" % (val, err, eps) ) 

**结果只有1次迭代**：

At iteration 1 (N=4), val=14.0000000000000000 prev=7.0000000000000000 error=2.333333e+00
Calculated value: 93333333333.3333435058593750  error: 2.333333e+00  for an epsilon of: 1.000000e-10

多次迭代的预期结果如下：

At iteration 1 (N=...), val=........ prev=...... error=1.552204e-10
At iteration n (N=...), val=........ prev=...... error=3.880511e-11
Calculated value: 0.0000000000  error: 3.880511e-11  for an epsilon of: 1.000000e-10

我想了解一些代码，因为这个问题已经被问到了其他值n，而不是数组