Transcript Document

Revisiting Numerical Integration: Getting More from Fewer Points

Bruce Cohen Lowell High School, SFUSD [email protected]

http://www.cgl.ucsf.edu/home/bic Asilomar - December 2007 David Sklar San Francisco State University [email protected]

Ver. 3.00

Romberg Integration

The method arises from a technique called Richardson extrapolation which can be used whenever the error

E

(

h

) can be expanded in a series of the form 

c h

1 

c h

2 2 

c h

3 3 

c h

4 4  We’ll illustrate Richardson’s technique by applying it to the trapezoidal rule.

I

 

b a

  , with    In fact we can show that if

f

subinterval then

E

can be expanded in a Taylor series on each can be expanded in a series of the form 

c h

2 2 

c h

4 4 

c h

6 6 

c h

8 8 

Romberg Integration

To implement the method we don’t need to know the coefficients, we need only

I

know that they exist. Assume that we have computed trapezoidal estimates for for

h, h/2, h/4, h/8,

       

looking at the first two we have

I

   

c h

2 2 

c h

4 4 

c h

6 6 

c h

8 8  and

I

T

E

T

c

2

h

2 2 2 

c

4

h

4 2 4 

c

6

h

6 2 6 

c

8

h

8 2 8  We can eliminate the

h

2 term by multiplying the second by 2 2 and subtracting the first giving  2 2  1 

I

 2 2

T

 

c

4 

h

4 2 4 

c

6 

h

6 2 6 

c

8 

h

8 2 8  Which we can rewrite as  2 2  1 

I

  2 2  1 

T

T

 

c

4 

h

4 2 4 

c

6 

h

6 2 6 

c

8 

h

8 2 8 

Romberg Integration

The last equation is equivalent to

I

T

T h

2  2 2  1  Which is in the form

I

 estimate  error 

k

4

h

4 2 4 

k

6

h

6 2 6 

k

8

h

8 2 8   2   

k

4

h

4 2 4 

k

6

h

6 2 6 

k

8

h

8 2 8 

T

2 (

h

) is an estimate for

I

with an error that is

O

(

h

4 )

Romberg Integration

We can repeat the process beginning with

I

 2   

k

4

h

4 2 4 

k

6

h

6 2 6 

k

8

h

8 2 8 

I

 3

 

k

 6

h

2 6 6

 

T

2  

k

8 

h

8 2 8 

T

2

h

2  2 4  1 

 

Continuing in this way we get a sequence of estimates for the integral

I

1

 

 ,

T h j

 

T j

 1 Where for each estimate we have

I

 

T j

 1 

h

2 2 2 

j

 1 

T j

 1  1 

j

  

E j

, with

E j

Romberg Integration -- An example

To see how Romberg integration works we organize our estimates into an array as follows

T

1 1    

T j

 1 

T j

 1  2

h

2 2 

j

 1 

T j

 1  1 

T

1

T

2 3  

T

1

T

2

T

3 4   The accuracy increases as we move down and as we move to the right.

Next we perform Romberg integration for a familiar example.  1 3 We place our trapezoidal estimates in the first column. The first estimate corresponds to a partition of [1, 3] into 4 subintervals, an

h

value of ½. Accuracy improves as we move down the column  1 3  2 5

x

5 2 1 3  2 5  3 5 2  1 5 2   5.835382907

5.858234

5.841100

5.835388

5.836813

5.835384

5.835384

5.835740

5.835382

5.835382 5.835382

Romberg Integration -- An example

We can improve our Romberg result if we start with more precise (though quite inaccurate) trapezoidal estimates in the first column. The averaging method used to get the more accurate estimates involves the subtraction of nearly equal numbers which can lead to precision loss as the computation progresses. The distinction between precision and accuracy in computation is important in practical numerical analysis.

5.858234

5.841100

5.835388

 1 3 3 2

x dx

 5.835382907

5.85823385920

5.836813

5.835384

5.835384

5.84110042206

5.83538927635

5.835740

5.835382

5.835382 5.835382

5.83681259147

5.83538331460

5.83538291715

5.83574034752

5.83538293287

5.835382907425.83538290726

Note : The correct value to 12 significant digits is 5.83538290725

Gaussian Quadrature

Gaussian Quadrature

In Gaussian quadrature the approximation to the integral is also a weighted sum of function values. The evaluation points (which are not equally spaced ) and weights are chosen to provide a very efficient estimate.

One way to implement Gaussian quadrature is to first transform the integral to an integral from –1 to 1, and then, for appropriate numbers of points

n

, get the gauss points and weights from the table. For our example we have

I

  1 3 3 2

x dx

  1  1 

u

 2  3 2

du

, using the linear substitution

u

2 So using

n

 5 points

I

   1 1   1  1 

u

 2  3 2

du

i

5   1

i

 

i

 .23692885

f

  .906179846

.5688888889

f

   .478628671

.23692885

f f

  .538469310

 .906179846

   .478628671

f

 .538469310

  5.8353829 which is correct to 7 decimal places.

Bibliography

1. M. Abramowitz, I. Stegun,

Handbook of Mathematical Functions

, Dover, New York, 1965 2. R. Burden, J. Faires,

Numerical Analysis

, 7 th Edition, Brooks/Cole, 2001 3. F. B. Hildebrand,

Introduction to Numerical Analysis

, McGraw-Hill, New York, 1956 4. http://numericalmethods.eng.usf.edu/anecdotes/romberg.doc

5. http://www-history.mcs.st-andrews.ac.uk/Biographies/ 6.

R Development Core Team,

R: A language and environment for statistical computing

, R Foundation for Statistical Computing, Vienna, Austria, 2006, 7.

8.

L. F. Richardson and J. L. Gaunt,

The deferred approach to the limit,

Philosophical Transactions of the Royal Society of London 226A, 1927 W. Romberg

, W. (1955), "Vereinfachte numerische Integration",

Norske Videnskabers Selskab Forhandlinger (Trondheim)

28

(7), 1955