Project 5: Different Methods of Numerical Integration

Due Dates

Wednesday, April 9, at 4:00 PM: Pseudocode
Monday, April 14, at 9:00 AM: The finished MATLAB code and analysis

Your job in this project: write programs to

integrate the Planck function over a given interval, using three different methods
integrate a given function of one variable over a given interval, using three different methods

You may not use MATLAB's built-in functions for finding maxima or performing numerical integration -- instead, please implement these algorithms yourself.

You'll need to modify one of your existing MATLAB source code files this week, plus write three more. Don't worry -- they are all very much alike, and you'll be able to cut and paste.

planck_integ.m (modify this one from last week)
- (you don't need to hand this one in, just use it to calculate some values as described below)
rect_integ.m
trap_integ.m
simpson_integ.m

Once again, you may use my Planck function in this assignment:

planck.m

Integrating the Planck function with different methods

For last week's assignment, you wrote a function which integrated power per unit area radiated by a blackbody at a given temperature:

 function energy = planck_integ(start_lambda, end_lambda, temperature, epsilon)

This week, I want you to consider a single case (the Sun, over the visible range of wavelengths):

temp 5,600 K, start_lambda 400 nm, end_lambda 800 nm

Set epsilon = 1e-6, as before.

Questions:

Perform this integration using three different methods:
1. rectangle method (same as last week)
2. trapezoid method
3. Simpson's method
For each of the three methods, you must divide the interval into 2, 4, 8, 16, and 32 little pieces. Make a table which shows for each method
```
     # of slices                      fractional change      
                                    since previous iteration
   ----------------------------------------------------------
        2
        4
        8
       16
       32
```
Using your tables as evidence, answer the following question:
- Which method converges most rapidly?
- How fast does each method converge? That is, if you double the number of slices, does the fractional change go down by a factor of 2? Or a factor of 4? Or what?
Finally, for each method, find out how many pieces into which you must break the interval in order to reach a fractional change of 1e-6. Make a table which shows the number the pieces, and the time it took your program to reach this point and evaluate the integral, for each method. That is, include the time it takes your program to iterate until it reaches the number of pieces required to satisfy the criterion. Which is fastest?

Integrate a given function of one variable over a given interval, using three different methods

Now, with the experience you've gained so far, it should be easy to write a program which integrates an arbitrary function of one variable, instead of the Planck function. The user will pass the name of the function to your program, which will then evaluate it over a given interval.

In fact, write three programs:

 function integ = rect_integ(start_x, end_x, func, epsilon)

 function integ = trap_integ(start_x, end_x, func, epsilon)

 function integ = simpson_integ(start_x, end_x, func, epsilon)

in which

          start_x             is the starting point of the interval over
	                             which to integrate

          end_x               is the ending point of the interval over
	                             which to integrate

          func                is the name of a function to integrate

          epsilon             denotes the termination criterion: keep going
                                until the fractional change in the
                                integral drops to this level

          integ               (the output argument) is the result of 
	                             the integration

For example, I might integrate sine(x) over the interval x=0 to x=10 by calling

          a = rect_integ(0, 10, 'sin', 1e-6);

or I might integrate a simple little function stored in file simple.m like this:

          function y = simple(x)
          y = sin(100*x)*exp(-x);

by calling

          a = rect_integ(0, 10, 'simple', 1e-6);

I suggest that you test your functions by using them to integrate simple functions, to which you know the correct answer.

Use each of your functions to integrate tan(x) over the interval x=0.1 to x=1.4, where "x" is in radians (this is how MATLAB does all its trigonometry). Make a table which shows

the number of pieces required to reach fractional changes of 1e-2, 1e-4, 1e-6
the time it takes to perform the integral in each case (again, include all the time it takes your program to work its way up to the critical number of pieces)

Make a graph which shows the fractional change on the x-axis, and the time required to find the integral to that precision on the y-axis. Use logarithmic axes in both directions. Plot the points for each method together on the same graph, so that you can compare the methods directly. What does the graph tell you about the efficiency of the different methods?

Bells and Whistles

Make a graph with MATLAB which shows the Planck function between 400 nm and 800 nm. Show the approximation to the area under the curve using the rectangular method, for N = 16 slices.
Find a simple formula which yields the approximate number of pieces, N, the rectangular method needs to integrate tan(x) over the interval x=0.1 to x=1.4 to a fractional error "epsilon". That is, some function which approximates N = f(epsilon). Then do the same for Simpson's method. Comment on the two functions.
Integrate the function
```
         y = sin(100*x) * exp(-x)
```
(where the argument of sin is in radians, as usual) over the range from 0 to 10. Calculate the answer analytically and compare to your numerical method. How many pieces do you need to get an accurate result? Why is this function so hard to integrate?

This page maintained by Michael Richmond. Last modified Apr 16, 2003.