4005.sagews

The trapezoid approximation: $f$ is the function, $[a,b]$ the interval, $k$ the number of subintervals.

def trapezoid(f, a, b, k):
    sum = 0;
    for i in range(k+1):
        x = a + (b - a)*i/k;
        y = f(x);
        if (i == 0 or i == k):
            y = y/2;
        sum += y;
    sum = sum * (b-a) / k
    return sum

The trapezoid error: $M$ is an upper bound on the second derivative of $f$.

def trapezoid_error(M,a,b,k):
    return M*(b-a)^3/12/k^2

An example: Let $f=e^{-x^2}$, on the interval $[0,1]$.

f(x)=exp(-x^2); x0=0; x1=1;

Compute and plot the second derivative to find the maximum value of $|f^{(2)}|$.

fpp=diff(f,x,2); plot(fpp,x0,x1)

So the maximum appears to be 2. Now compute the approximation until the error is small; try increasing $j$ until the two values $A-E$ and $A+E$ round to the same value to two decimal places.

j=6; A=n(trapezoid(f,x0,x1,j)); E=n(trapezoid_error(2,x0,x1,j)); A; E; A-E; A+E;

0.745119412436179 0.00462962962962963 0.740489782806550 0.749749042065809

The Simpson's rule function and error computation are similar. You must put in an even number for $k$; the function doesn't check for this.

def simpson(f, a, b, k):
    sum = 0;
    for i in range(k+1):
        x = a + (b - a)*i/k;
        y = f(x);
        if (i > 0 and i < k):
            y = 2*(i%2+1)*y
        sum += y;
    sum = sum * (b-a) / k /3
    return sum

def simpson_error(M,a,b,k):
    return M*(b-a)^5/180/k^4

f4p=diff(f,x,4); plot(f4p,x0,x1)

The maximum value of $|f^{(4)}|$ appears to be 12.

j=4; A=n(simpson(f,x0,x1,j)); E=n(simpson_error(12,x0,x1,j)); A; E; A-E; A+E;

0.746855379790987 0.000260416666666667 0.746594963124320 0.747115796457654

Let's see what Sage thinks the value is.

n(integral(f,x,x0,x1))

0.746824132812427