Andy R Terrel - CMSC 16100 Lab 3 - Numerical Integration

Simple Numerical Integration rules

Numerical integration routines are as old as Archimedes and are a great example of treating functions as higher order objects. The difference between these different rules come in the number of function evaluations and the order of error dependent on the step size. For this lab, you will be implementing several simple integration rules, use a routine to see the difference in error, and then use these rules with a heirarchical algorithm.

To compute the integral,

I = ∫_a^b f(x) dx,

We will approximate it with a series of sums, based upon function evaluation. For example:

Midpoint rule: I_h = h Σ f(x_i+1/2) , for i=0,1,2,...,n-1
Trapezoid rule: I_h = h/2 Σ f(x_i) + f(x_i+1), for i=0,1,2,...,n-1
Simpson 1/3 rule: I_h = h/3 Σ f(x_i-1) + 4 f(x_i) + f(x_i+1), for i=1,3,5,...,n-1

For each rule make an integration function that takes the following inputs and returns I_h.

function (f(x)),
start (a),
stop (b),
k (where h=(b-a)/2^k),

Checking Convergence

To numerically see the trade off between order of accuracy and number function evaluations, create a function to measure the error different rules give on with a range of values for k. The function should accept an integration rule that it can input the testing f(x). For an example, you can integrate x⁷ from 0 to 1 which has the exact solution .125. For each integration rule, add to your Test Cases section an evalution with k = 1,2,3, and separate the different integration rules with a display statement, for example:

(display "Midpoint\n")
(chk_error ...)
(chk_error ...)
(chk_error ...)
(chk_error ...)
(display "Trapezoid\n")
(chk_error ...)
...

Romberg's Rule

Romberg's rule is a heirarchical method that for each level corrects with the solution of the previous level. It repeats this process m times, and for m=0 applies the provided integration rule with 2^k points. Thus the algorithm has two cases:

R(I_h,f,0,k,a,b) = I_h(f, k, a,b)
R(I_h,f,m,k,a,b) = D (C * R(I_h,f,m-1,k,a,b) - R(I_h,f,m-1,k-1,a,b))
where C = 4^m and D = 1/(C-1)

Please note k must be larger than m, have your code display an error (search the helpdesk) if this inequality is not satified. Implement Romberg's rule, and add calls to your chk_error function for (m,k) = (1,2),(2,3),(2,5),(3,4) for each simple rule (again separated by the display)

Extra-Credit: Integrating Mulitple Dimensions

Doing multiple dimensions can be done applying the rule to one direction at a time. This can be done in our setup by accumulating the error for a value along one direction (by applying the appropriate I_h), then combining the results these results for the other direction is the same way as done for one direction.

Make a two dimensional Midpoint and Trapezoid rule, add an error checker, and add appropriate Test Cases.

Turning Everything in

Add the appropriate test cases to the end of your definitions window so I can see your examples:

  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  ;;                            Test Cases                                     ;;
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

Save the definitions window as (your_cnet_id)-Lab03.scm and put it in the Lab 3 assignment on chalk. So for example I would save my file as aterrel-Lab03.scm. Also do yourself a favor and save your file somewhere you can look back at it (email yourself or store it on your own media).