Unformatted text preview: Numerical Integration Background Trapezoid Simpson’s 1/3 rule Other methods: Monte Carlo HW#4 due today. You should be reading in Appendix A.3 Integration is the inverse operation to differentiation. Integral of a function is a measure of the “area” bounded by the function and the horizontal axis. Many analytic functions can integrated analytically to get a closed form solution. Area above axis is (+), area below axis is (). b f (x)dx = area a 1 1 dx = ln(x) xm dx = xm+1 x m+1 sin(x)dx = −cos(x) π sin(x)dx = 2.00
0 b a 1 m+1 x dx = x m+1
m b a On a computer, integration becomes summation. Numerical integration is also called quadrature. Break integral range into multiple “panels” with a simple shape. Add areas of all panels to get total area. NewtonCotes Method: top of each panel defined by a polynomial of order n. n=1: line  Trapeziod method n=2: parabola – Simpson’s 1/3 n=3: cubic – Simpson’s 3/8 For many applications, Trapezoid method is good enough. Geometry of panels is trapezoidal. Area of a trapezoid: Total area for N panels: h [f (a + h) + f (a)] 2 I = I1 + I2 + I3 + ... + IN N −1 h [f (a + ih) + f (a + (i + 1)h)] = 2 i=1 def trap(f,xmin,xmax,numPanels): Isum=0. h=(xmaxxmin)/float(numPanels) for i in range(numPanels): x=xmin+i*h Isum+=(f(x+h)+f(x))*h/2. return Isum
Trapezoid error is of order ~h2 Trapezoid method is can be coupled with Richardson extrapolation to improve error further known as the Romberg method Panel cap is a quadratic curve rather than a straight line. Requires evaluation at 3 points rather than 2. Error is of order ~h4. b a a+b h f (x)dx f (a) + 4f ( ) + f (b) 2 3 def simp13(f,xmin,xmax,numPanels): Isum=0. h=(xmaxxmin)/float(numPanels) for i in range(0,numPanels,2): x=xmin+i*h Isum+=(f(x)+4*f(x+h)+f(x+2*h))*h/3 return Isum Use the trapezoid method to estimate the integral: Compute and plot the % error vs number of panels. 1 e −x 0 dx = −[e −1 − e ] = 0.632120558829
0 ...
View
Full Document
 Spring '09
 Gladden
 trapezoid method, Isum, Integration Background

Click to edit the document details