Approximating definite integrals using midpoint Riemann sums and Simpson's  Rule.

I. Midpoint Riemann sum approximations.

f[x_] = 1/(1 + x^2) ;

a = 0 ;

b = 1 ;

MidpointSum[n_] := Sum[f[a + (2i - 1) * (b - a)/(2 * n)], {i, 1, n}] * (b - a)/n ;

Kmid = Maximize[Abs[D[f[x], {x, 2}]], x≥a&&x≤b, x][[1]] ;

MidpointError[n_] := Kmid * (b - a)^3/(24 * n^2) ;

Note: D[f[x],{x,2}] gives the second derivative of f with respect to x, same as f''[x].

N[MidpointSum[100], 20]

N[MidpointError[100], 20]

0.78540024673078116242

8.3333333333333333333*10^-6

Notice that the above shows that Pi/4 equals .7854002... with an
accuracy guaranteed to 5 decimal places. Let's check the acutal
difference:

N[Pi/4 - MidpointSum[100], 20]

-2.0833333328528025817*10^-6

The following took a couple of minutes of runtime on my PC. If you don't want to wait on a calculation like this--from "Kernel" on
the main MATHEMATICA menu choose "Abort Evaluation".

N[MidpointSum[20000], 20]

N[MidpointError[20000], 20]

N[Pi/4 - %%, 20]

0.78539816344953164295

2.0833333333333333333*10^-10

-5.208333333*10^-11

Analyzing approximations for different values of n in tabular form is sometimes useful for getting an idea of the closeness of the approximations:

tt = Table[Print[N[MidpointSum[100 * k], 12], "  ", N[MidpointError[100 * k]]], {k, 1, 10}] ;

0.785400246731  8.33333*10^-6

0.785398684231  2.08333*10^-6

0.785398394879  9.25926*10^-7

0.785398293606  5.20833*10^-7

0.785398246731  3.33333*10^-7

0.785398221268  2.31481*10^-7

0.785398205914  1.70068*10^-7

0.785398195950  1.30208*10^-7

0.785398189118  1.02881*10^-7

0.785398184231  8.33333*10^-8

II. Approximation using Simpson's Rule.

f[x_] = 1/(1 + x^2) ;

a = 0 ;

b = 1 ;

SimpsonSum[n_] := (4 * Sum[f[a + (2i - 1) * (b - a)/(2 * n)], {i, 1, n}] + f[a] + f[b] + 2 * Sum[f[a + (2i) * (b - a)/(2 * n)], {i, 1, n - 1}]) * (b - a)/(6 * n) ;

Ksimp = Maximize[Abs[f''''[x]], x≥a&&x≤b, x][[1]] ;

SimpsonError[n_] := Ksimp * (b - a)^5/(180 * n^4) ;

N[SimpsonSum[100], 20]

N[SimpsonError[100], 20]

0.78539816339744815461

1.3333333333333333333*10^-9

In this example Ksimp equals 24. Here's a picture showing it:

Ksimp

Simplify[D[1/(1 + x^2), {x, 4}]]

Plot[%, {x, -2, 2}, ImageSize→500] ;

Ksimp

(24 (1 - 10 x^2 + 5 x^4))/(1 + x^2)^5

[Graphics:HTMLFiles/ComputerLab3_45.gif]

Table[Print[N[SimpsonSum[100 * k], 10], "  ", N[SimpsonError[10 * k]]], {k, 1, 15}] ;

0.7853981634  0.0000133333

0.7853981634  8.33333*10^-7

0.7853981634  1.64609*10^-7

0.7853981634  5.20833*10^-8

0.7853981634  2.13333*10^-8

0.7853981634  1.02881*10^-8

0.7853981634  5.55324*10^-9

0.7853981634  3.25521*10^-9

0.7853981634  2.03221*10^-9

0.7853981634  1.33333*10^-9

0.7853981634  9.10685*10^-10

0.7853981634  6.43004*10^-10

0.7853981634  4.66837*10^-10

0.7853981634  3.47078*10^-10

0.7853981634  2.63374*10^-10

Created by Mathematica  (December 5, 2006) Valid XHTML 1.1!