使用梯形法则计算积分（积分的近似值和迭代次数）

Question

程序需要计算定义积分与预定 梯形规则的准确性（eps），我的函数需要返回：

1.积分的近似值。

2.迭代次数。

我的代码：

    from math import *
    def f1(x):
       return (x ** 2 - 1)**(-0.5)

    def f2(x):
       return (cos(x)/(x + 1))

    def integral(f,a,b,eps):
       n = 2
       x = a
       h = (b - a) / n
       sum = 0.5 * (f(a) + f(b))
       for i in range(n):
          sum = sum + f(a + i * h)
       sum_2 = h * sum
       k = 0
       flag = 1

       while flag == 1:
          n = n * 2
          sum = 0
          k = k + 1
          x = a
          h = (b - a) / n
          sum = 0.5 * (f(a) + f(b))
          for i in range(n):
             sum = sum + f(a + i * h)
          sum_new = h * sum
          if eps > abs(sum_new - sum_2):
             t1 = sum_new
             t2 = k
             return t1, t2
          else:
             sum_2 = sum_new

   x1 = float(input("First-begin: "))
   x2 = float(input("First-end: "))
   y1 = float(input("Second-begin: "))
   y2 = float(input("Second-end: "))
   int_1 = integral(f1,x1,y1,1e-6)
   int_2 = integral(f2,x2,y2,1e-6)
   print(int_1)
   print(int_2) 

它不能正常工作。请帮忙！

Answer 1

您的数学运算错误。错误在行中

for i in range(n):
    sum = sum + f(a + i * h)

range(n) 总是从 0 开始，因此在您的第一次迭代中，您只需再次添加 f(a) 术语。 如果您将其替换为

for i in range(1, n):
      sum = sum + f(a + i * h)

它有效。

此外，您还有大量冗余代码；您基本上对集成算法的核心进行了两次编码。尝试关注DRY-principle。

Answer 2

积分的梯形规则简单地说，积分 $\int_a^b f(x) dx$ 的近似值是(b-a) (f(a)+f(b))/2。误差与(b-a)^2 成正比，因此可以使用复合 规则进行更好的估计，即将初始区间细分为多个较短的区间。

是否可以使用更短的间隔并仍然重用先前计算的函数值，从而最大限度地减少函数评估的总数？

是的，如果我们将每个区间分成两等份是可能的，因此在阶段0 我们使用1 个区间，在阶段1 2 等区间，一般来说，在阶段n，我们使用2ⁿ 个等间隔。

让我们从一个简单的问题开始，看看是否可以概括这个过程……

a, b = 0, 32
L = b-a = 32

根据梯形规则，初始近似值是I0，由下式给出

I0 = L * (f0+f1)/2
   = L * S0

S0 = (f0+f1)/2;下面是实轴、区间极值坐标和评估函数的图形表示

x0                             x1
01234567890123456789012345679012
f0                             f1

接下来，我们将原来的区间一分为二，

L = L/2

x0              x2             x1
01234567890123456789012345679012
f0              f2             f1

新的近似值，阶段n=1，是使用梯形规则的两倍并应用一点代数得到的

I1 = L * (f0+f2)/2 + L * (f2+f1)/2
   = L * [(f0+f1)/2 + f2]
   = L * [S0 + S1]

S1 = f2

另一个细分，舞台n=2、L = L/2和

x0      x3      x2      x4      x1  
012345678901234567890123456789012
f0      f3      f2      f4      f1

I2 = L * [(f0+f3) + (f3+f2) + (f2+f4) + (f4+f1)] / 2
   = L * [(f0+f1)/2 + f2 + (f3+f4)]
   = L * [S0+S1+S2]

S2 = f3 + f4.

不难，鉴于这张图，

x0  x5  x3  x6  x2  x7  x4  x8  x1
012345678901234567890123456789012
f0  f5  f3  f6  f2  f7  f4  f8  f1

了解我们的下一个近似值可以如下计算

L = L/2
S3 = f5+f6+f7+f8
I3 = L*[S0+S1+S2+S3]

现在，我们必须了解如何计算Sn, n = 1, … 的泛化——对我们来说，伪代码是

L_n = (b-a)/2**n
list_x_n = list(a + L_n + 2*Ln*j for j=0, …, 2**n-1) 
Sn = Sum(f(xj) for each xj in list_x_n)

对于n = 3，L = (b-a)/8 = 4，我们从上面的公式得到list_x_n = [4, 12, 20, 28]，请看图...

现在我们准备好用 Python 编写算法了

def trapaezia(f, a, b, tol):
     "returns integ(f, (a,b)), estimated error and number of evaluations"

     from math import fsum # controls accumulation of rounding errors in sums

     L = b - a
     S = (f(a)+f(b))/2
     I = L*S
     n = 1
     while True:
         L = L/2
         new_points = (a+L+j*L for j in range(0, n+n, 2))
         delta_S = fsum(f(x) for x in new_points)
         new_S = S + delta_S
         new_I = L*new_S
         # error is estimated using Richardson extrapolation (REP)
         err = (new_I - I) * 4/3
         if abs(err) > tol:
             n = n+n
             S, I = new_S, new_I
         else:
             # we return a better estimate using again REP
             return (4*new_I-I)/3, err, n+n+1

如果您对 Richardson 外推法感到好奇，我推荐 this document，它正好处理 REP 在梯形规则正交算法中的应用。

如果您对math.fsum、the docs 感到好奇，请不要说太多，但the link 对原始实现也包括对所有相关问题的扩展解释。