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 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 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 0.746824132812427 |
|
|