The foundamental idea of numerical integration is to estimate the area of the region in the *xy*-plane bounded by the graph of function *f(x)*. The integral was esimated by divide *x* to small intervals, then add all the small approximations to give a total approximation.

Numerical integration can be done by trapezoidal rule, simpson's rule and quadrature rules. R has a built-in function integrate, which performs adaptive quadrature.

Trapezoidal rule works by approximating the region under the graph *f(x)* as a trapezoid and calculating its area.

^{?}View Code RSPLUS

1 2 3 4 5 6 7 8 9 10 | trapezoid < - function(fun, a, b, n=100) { # numerical integral of fun from a to b # using the trapezoid rule with n subdivisions # assume a < b and n is a positive integer h <- (b-a)/n x <- seq(a, b, by=h) y <- fun(x) s <- h * (y[1]/2 + sum(y[2:n]) + y[n+1]/2) return(s) } |

Simpson's rule subdivides the interval [a,b] into *n* subintervals, where *n* is even, then on each consecutive pairs of subintervals, it approximates the behaviour of *f(x)* by a parabola (polynomial of degree 2) rather than by the straight lines used in the trapezoidal rule.

^{?}View Code RSPLUS

1 2 3 4 5 6 7 8 9 10 11 12 13 | simpson < - function(fun, a, b, n=100) { # numerical integral using Simpson's rule # assume a < b and n is an even positive integer h <- (b-a)/n x <- seq(a, b, by=h) if (n == 2) { s <- fun(x[1]) + 4*fun(x[2]) +fun(x[3]) } else { s <- fun(x[1]) + fun(x[n+1]) + 2*sum(fun(x[seq(2,n,by=2)])) + 4 *sum(fun(x[seq(3,n-1, by=2)])) } s <- s*h/3 return(s) } |

## Recent Comments