**Numerical Integration**

**Introduction
**

- Given a function y=f(x)

- Want area under curve between x=a and x=b

- Assume we cannot obtain an exact solution, or perhaps f(x) is unknown, and we only have data points (xi,yi).

- Main idea: Break area under curve into n
panels of width h: A = Ai , h = (b-a)/n

- How to calculate Ai ?

- Many methods of approximation

- All estimates improve with larger n, as
h0 .

**Midpoint Rule
**

- Use rectangular panels (i.e. assume f(x) is constant and equal to y* between xi-1 and xi )

- Three possibilities:

- Left endpoint: y* = yi-1 = f(xi-1)
- Right endpoint: y* = yi = f(xi) = f(xi-1 + h)
- Midpoint: y* = yi - ½ = f( (xi-1 + xi)/2 ) = f(xi-1 + h/2)
- On average, midpoint looks like best guess

**Midpoint Rule Algorithm
**We wish to integrate y=f(x)
from x=a to x=b using n rectangular panels of equal width.

h = (b-a) /n

x0 = a

xi = xi-1 + h

Ai = h* f( (xi-1 + xi)/2 )

A = Sum of Ai (i=1 -> i=n)

A = h* [Sum of f( (xi-1 + xi) /2 )] (i=1 -> i=n)

REAL FUNCTION MidPt(x0, x1, N) ! Estimate the area under F between x0 and x1 ! using the Midpoint Rule with N panels. REAL x0, x1, x, h, Sum INTEGER N h = (x1 - x0) / N Sum = 0.0 ! Add area of each panel, assuming unit width ! and multiplying by the true width h later. x = x0+h/2.0 ! Initial start is middle of first panel L1: DO k= 1,N Sum = Sum + F(x) x = x + h END DO L1 MidPt = h*Sum RETURN

C Version

double midpoint(double x0, double x1, int n) { double h, sum; /* Calculate panel width h; initialize Sum. */ h= (x1 - x0) / n; x0= x0 + h/2; sum= 0; /* Add area of each panel, multiplying by the true width h only once later. */ do { sum= sum + f(x0); x0= x0 + h; }while(x0 x1); return(sum * h); }

**Trapezoidal Rule
**

Try taking average of the expressions for rectangular
Ai, with f evaluated at left and right endpoints

Ai = [ h f(xi-1) + h f(xi) ]

Ai = h/2 [ f(xi-1) + f(xi) ]

Approximate Ai, as the area of trapezoid abcd.

Trapezoidal method appears to take 2 function evaluations per
panel, but

Ai + Ai+1 = h/2 [ f(xi-1) + f(xi) ] + h/2 [
f(xi) + f(xi+1) ]

Ai + Ai+1 = h/2 [ f(xi-1) + 2f(xi) + f(xi+1)
]

For n panels, need only n+1 evaluations of f.

**Trapezoidal Rule Algorithm
**

We wish to integrate y=f(x) from x=a to x=b
using n trapezoidal panels of equal width.

h= (b-a)/n

x0 = a

xn = b

xi = xi-1 + h

Ai = h/2 [ (f(xi-1) + f(xi) ]

A Ai

n - 1

A h [ ( f(x0) + f(xn) ) + f(xi) ]

i = 1

**Simpson's Rule
**

In Trapezoidal rule, successive points (xi, yi) were joined with straight lines

Now, fit a parabola to 3 points (xi-2, yi-2), (xi-1, yi-1), and (xi, yi) for better accuracy

One panel (width 2h) now includes 3 points

h = (b-a)/n

b-a = nh

n must therefore be even

From geometry or Taylor series expansion, we find

Ai = h/3 [ f(xi-2) + 4f(xi-1) + f(xi) ] , for i = 2,4,6,…,n

Once more, combine expressions for adjacent panels and simplify.

Ai + Ai+2 = h/3 [ f(xi-2) + 4f(xi-1) + 2f(xi)
+ 4f(xi+1) + f(xi+2) ]

For n panels, again need only n+1 evaluations of f.

**Simpson's Rule Algorithm
**

We wish to integrate y=f(x) from x=a to x=b
using n/2 double parabolic panels of equal width.

h = (b-a)/n

x0= a

xn= b

xi= xi-1 + h

Ai = h/3 [ f(xi-2) + 4f(xi-1) + f(xi) ]

A Ai

n - 1 n - 2

A h/3 [ f(x0) + f(xn) + 4 f(xi) + 2 f(xi) ]

i = 1 i = 2

(i odd) (i even)

**Adaptive Integration Procedure
**

- Arbitrary precision

**Advanced, recursive adaptive integration
technique**

/* more sophisticated adaptive integration method

using recursion */

Go back to lecture menu

Go back to main page