Computing integral with the Trapezoidal Rule (the approximate value of the integral, and the number of iterations)

Question

The program needs to compute define integral with a predetermined accuracy (eps) with the Trapezoidal Rule and my function needs to return:

1.the approximate value of the integral.

2.the number of iterations.

My code:

    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) 

It doesn't work correct. Help, please!

Answer 1

You implemented the math wrong. The error is in the lines

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

range(n) always starts at 0, so in your first iteration you just add the f(a) term again. If you replace it with

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

it works.

Also, you have a ton of redundant code; you basically coded the core of the integration algorithm twice. Try to follow the DRY-principle.

Answer 2

The trapezoidal rule of integration simply says that an approximation to the integral $\int_a^b f(x) dx$ is (b-a) (f(a)+f(b))/2. The error is proportional to (b-a)^2, so that it is possible to have a better estimate using the composite rule, i.e., subdividing the initial interval in a number of shorter intervals.

Is it possible to use shorter intervals and still reuse the function values previously computed, so minimizing the total number of function evaluation?

Yes, it is possible if we divide each interval in two equal parts, so that at stage 0 we use 1 intervals, at stage 1 2 equal intervals and in general, at stage n, we use 2ⁿ equal intervals.

Let's start with a simple problem and see if it possible to generalize the procedure…

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

by the trapezoidal rule the initial approximation say I0, is given by

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

with S0 = (f0+f1)/2; a pictorial representation of the real axis, the coordinates of the interval extremes and the evaluated functions follows

x0                             x1
01234567890123456789012345679012
f0                             f1

Next, we divide the original interval in two,

L = L/2

x0              x2             x1
01234567890123456789012345679012
f0              f2             f1

and the new approximation, stage n=1, is obtained using two times the trapezoidal rule and applying a bit of algebra

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

with S1 = f2

Another subdivision, stage n=2, L = L/2 and

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]

with S2 = f3 + f4.

It is not difficult, given this picture,

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

to understand that our next approximation can be computed as follows

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

Now, we have to understand how to compute a generalization of Sn, n = 1, … — for us, the pseudocode is

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)

For n = 3, L = (b-a)/8 = 4, we have from the formula above list_x_n = [4, 12, 20, 28], please check with the picture...

Now we are ready to code our algorithm in 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

If you are curious about Richardson extrapolation, I recommend this document that deals exactly with the application of REP to the trapezoidal rule quadrature algorithm.

If you are curious about math.fsum, the docs don't say too much but the link to the original implementation that also includes an extended explanation of all the issues involved.