MANE 3351
Lecture 15
Classroom Management
Agenda
- Romberg Integration
- Lab 7, if not completed
- Homework 4 due Monday
Calendar
| Date | Lecture | Lab |
|---|---|---|
| 10/14 | Simpson's Rule, Homework 4 | Lab 7 |
| 10/16 | Romberg Integration | Lab 7, continued |
| 10/21 | No Class, Dr. Timmer on ABET Visit | No Class, Dr. Timmer on ABET Visit |
| 10/23 | Gaussian Quadrature, Homework 5 | Lab 8 |
| 10/28 | Numerical Differentiation (not on Test 2) | Lab 8, continued |
| 11/4 | Numerical Integration (not on Test 2) | Lab 9 |
| 11/6 | Linear Algebra, part 1 (not on Test 2) | Lab 9, continued |
| 11/11 | Test 2 (Root Finding and Numerical Integration) | No Lab |
Resources
Handouts
Lecture 15
Today's topic is Romberg Integration
- Clever combination of trapezoid rule and Richardson's Extrapolation
- Highly accurate
- Cheney and Kincaid (2004)1 show example output in the form of a lower triangle from Romberg integration
Step 1
The first step is to calculate \(R(0,0)\)
- \(R(0,0)\) is the result of applying the Trapezoid rule with 1 interval
- \(R(0,0)=\frac{1}{2}(b-a)\left[f\left(a\right)+f\left(b\right)\right]\)
- For our example of the standard normal pdf with \(a=-5\), and \(b=0\), we observe
- \(R(0,0)=\frac{1}{2}(0-(-5.0))\left[f(-5.0)+f(0.0)\right]=\frac{1}{2}(5.0)\left[0.0+0.398942\right]=0.997355\)
- This is a very poor approximation to the true value of 0.5
Step 2
Start a second row and calculate \(R(1,0)\). For each new row, double the number of intervals used in the trapezoid rule
- \(R(1,0)\) is the trapezoid with two intervals
- The general formula for \(R(n,0)\) is
\[
R(n,0)=\frac{1}{2}R\left(n-1,0\right)+h\sum_{k=1}^{2^{n-1}}f\left[a+\left(2k-1\right)h\right]
\]
where \(h=\left(b-a\right)/2^n\) and \(n\geq 1\)
Step 3
Complete the second row and calculate \(R(1,1)\)
- The calculation of \(R(1,1)\) utilizes Richardson's extrapolation, \(R(1,1)=f\left[R(1,0),R(0,0)\right]\)
- The general formula for \(R(n,m)\) is
\[
R(n,m)=R(n,m-1)+\frac{1}{4^m-1}\left[R(n,m-1)-R(n-1,m-1)\right]
\]
Error Analysis
- Cheney and Kincaid (2004)1 reports the following errors
- The error for the first column is \(\mathcal{O}(h^2)\)
- The error for the second column is \(\mathcal{O}(h^4)\)
- The error for the third column is \(\mathcal{O}(h^8)\)
- and so on
Pseudo-code
Cheney and Kincaid (2004)1 provide the following pseudo-code
Python Code for Romberg Integration
import math
import numpy as np
def f(z):
return (math.exp(-0.5*z**2)/((2.0*math.pi)**0.5))
#
a=float(input("Enter the lower limit of the integral: "))
b=float(input("Enter the upper limit of the integral: "))
n=int(input("enter the number of interations (n): "))
#
# initialize matrix r
r=np.zeros(shape=(n+1,n+1))
h=b-a
#find R(0,0)
r[0][0]=(h/2.0)*(f(a)+f(b))
for i in range(1,n+1):
h=h/2.0
sum=0.0
for k in range(1,2**i,2):
sum=sum+f(a+k*h)
r[i][0]=0.5*r[i-1][0]+sum*h
for j in range(1,i+1):
r[i][j]=r[i][j-1]+(r[i][j-1]-r[i-1][j-1])/(4**j-1)
print(r)