Error for trapezoidal method increases with N if integral equals 0

Question

I've written a simple implementation of the trapezoidal method to find the integral of sine:

def trapezoidal_method(a: float, b: float, n: int) -> float:
    length = (b - a)/n
    integral = 0
    start = a
    integral += math.sin(a)/2
    for _ in range(1, n):
        integral += math.sin(start + length)
        start += length
    integral += math.sin(b)/2
    return integral * length

It converges as expected for most situations... but it goes crazy when the result should be 0 (like integrating from -1 to 1):

How do I fix this? Tried explicitly casting n to float and tried using the decimal library, didn't change anything

but if you look at the scale factor of the y-axis of the 2nd plot is almost 0 (order of e-15) — cards, Dec 30 '22 at 18:13
Because of numerical instability. It is simply not possible to expect floating-point types to work well in this sort of environment, because of what they fundamentally **are**. if the linked duplicate doesn't adequately answer the question for you, it would be better rephrased as a question for [math.se]. — Karl Knechtel, Dec 30 '22 at 23:08
Probably the error increases with N because it is accumulated with the iterations. More about Python's [floating](https://stackoverflow.com/a/455634/16462878) numbers — cards, Dec 30 '22 at 23:20

cards · Answer 1 · 2022-12-30T23:35:44.837

That's (mainly) due to floating number precision. In the case of functions as sine, which are transcendental by nature, only an approximation of it is possible (or depends on its implementation).

By adding a (very quick!) a print of the x and y values at each term of the sequence you will see the odd symmetry is not respected... or at least with a small error. So, f(x) = -f(-x) + error.

import math


def trapezoidal_method(func, a: float, b: float, n: int) -> float:
    length = (b - a)/n
    integral = 0
    start = a
    integral += func(a)/2
    for _ in range(1, n):
        integral += func(start + length)
        start += length
    integral += func(b)/2

    # evaluation for each term of the sequence
    print(0, f'x={a}  y={func(a)}') 
    x = a
    for i in range(1, n):
        x += length
        y = func(x)
        print(i, f'{x=}  {y=}')     

    print(n, f'x={b}  y={func(b)}')
    return integral * length

A test with sin

a, b, n = -1, 1, 7
func = math.sin
tm = trapezoidal_method(func, a, b, n)

Output

0 x=-1                  y=-0.8414709848078965
1 x=-0.7142857142857143 y=-0.6550778971785186
2 x=-0.4285714285714286 y=-0.41557185499305205 # <-
3 x=-0.1428571428571429 y=-0.1423717297922637
4 x= 0.1428571428571428 y=0.1423717297922636
5 x= 0.4285714285714285 y=0.41557185499305194  # <-
6 x= 0.7142857142857142 y=0.6550778971785185
8 x= 1                   y=0.8414709848078965
-1.586032892321652e-16

Notice that also the x-values maybe slightly different!

By choosing suitable limits of integration and amount of strips then the in some cases the integral maybe 0. Here an example with the diagonal function

a, b, n = -1, 1, 6
func = lambda x: x
tm = trapezoidal_method(func, a, b, n)

Output

0 x=-3.0 y=-3.0
1 x=-2.0 y=-2.0
2 x=-1.0 y=-1.0
3 x=0.0 y=0.0
4 x=1.0 y=1.0
5 x=2.0 y=2.0
7 x=3.0 y=3.0
0.0

Error for trapezoidal method increases with N if integral equals 0

1 Answers1