Numerical differentiation and integration · Lecture 20

Numerical Integration

Areas, distances and probabilities from tables and from functions with no antiderivative: the trapezoidal rule, Simpson's 1/3 and 3/8 rules, their errors, Romberg integration and Gauss quadrature.

Md Mahbub E Noor · Assistant Professor
CSE-3107 Numerical Methods · B.Sc. 3rd Year 1st Semester
Department of CSE, University of Barishal

Why this lecture

The problem: distance from speedometer readings

A car accelerates from rest. Its speed is read every 10 seconds for one minute. How far did it travel?

  • Distance =∫060v(t) dt= \int_0^{60} v(t)\,dt, the area under the speed curve
  • There is no formula for v(t)v(t): only the readings
  • The same need appears for ∫01e−x2dx\int_0^1 e^{-x^2}dx (normal probabilities): a formula exists but no antiderivative
  • We need rules that turn a few function values into an accurate area, with an error estimate
tt (s)vv (m/s)
00
1012
2022
3030
4036
5040
6042
Md Mahbub E Noor
Lecture 20: Numerical Integration
University of Barishal
2 / 22

Today's goals

  • Derive the trapezoidal, Simpson 1/3 and 3/8 rules from Newton's forward formula
  • Apply the composite rules to tables and functions, and present the working clearly
  • State and use their error terms: O(h2)O(h^2) and O(h4)O(h^4)
  • Improve accuracy with Romberg integration; understand Gauss quadrature
  • Choose a rule for a given problem, and recognise when rules fail
Last time

Differentiation from tables: sensitive to noise and round-off.

Today

Integration from tables: stable, accurate, and the basis of the ODE methods that follow.

Numerical differentiation and integration · you are here
L19 · Numerical Differentiation→L20 · Numerical Integration
Builds on: L16, L19
Md Mahbub E Noor
Lecture 20: Numerical Integration
University of Barishal
3 / 22
Start here

The idea in one picture

xf(x)abarea under the curve = ?4 strips with straight tops (trapezoids)gapgapcurved tops (parabolas): Simpson. Gaps almost gone
An integral is the area under a curve.Cut the area into strips. One strip with a straight top is a trapezoid: easy area.Add the strips: that is the trapezoidal rule.Small gaps remain where the curve bends. More strips make them smaller.Simpson's rule uses curved tops (parabolas over two strips): far more accurate for the same work.
In plain words
  • An integral is an area
  • Cut it into thin strips whose areas are easy, and add them
  • Better-shaped strips (curved tops) need far fewer of them
New words today
quadrature
computing an integral numerically
trapezoidal rule
strips with straight tops
Simpson's rule
strips with parabolic tops
composite rule
the same rule repeated on many strips
Md Mahbub E Noor
Lecture 20: Numerical Integration
University of Barishal
4 / 22

The idea: integrate the interpolating polynomial

Integrate Newton's forward formula over n steps (x = x₀ + ph, dx = h dp)∫x0x0+nhy dx=h∫0n(y0+p Δy0+p(p−1)2Δ2y0+⋯)dp\int_{x_0}^{x_0+nh} y\,dx = h\int_0^n \Big(y_0 + p\,\Delta y_0 + \frac{p(p-1)}{2}\Delta^2 y_0 + \cdots\Big)dp
General quadrature formula (Newton–Cotes)=nh[y0+n2Δy0+n(2n−3)12Δ2y0+n(n−2)224Δ3y0+⋯]= nh\Big[y_0 + \frac n2\Delta y_0 + \frac{n(2n-3)}{12}\Delta^2y_0 + \frac{n(n-2)^2}{24}\Delta^3 y_0 + \cdots\Big]
  • n=1n = 1 (straight line): h2(y0+y1)\frac h2(y_0 + y_1), the trapezoidal rule
  • n=2n = 2 (parabola): h3(y0+4y1+y2)\frac h3(y_0 + 4y_1 + y_2), Simpson's 1/3 rule
  • n=3n = 3 (cubic): 3h8(y0+3y1+3y2+y3)\frac{3h}{8}(y_0 + 3y_1 + 3y_2 + y_3), Simpson's 3/8 rule
  • Then repeat the small rule along the table: the composite rules
Md Mahbub E Noor
Lecture 20: Numerical Integration
University of Barishal
5 / 22

The composite rules and their errors

Trapezoidal (any n): ends once, inside points twice, times h/2T=h2[y0+2(y1+⋯+yn−1)+yn],ET=−(b−a)h212f′′(ξ)T = \frac h2\big[y_0 + 2(y_1 + \cdots + y_{n-1}) + y_n\big], \qquad E_T = -\frac{(b-a)h^2}{12}f''(\xi)
Simpson 1/3 (n even): weights 1, 4, 2, 4, …, 2, 4, 1 times h/3S=h3[y0+4(y1+y3+⋯ )+2(y2+y4+⋯ )+yn]ES=−(b−a)h4180f(4)(ξ)\begin{aligned} S &= \frac h3\big[y_0 + 4(y_1 + y_3 + \cdots) + 2(y_2 + y_4 + \cdots) + y_n\big] \\ E_S &= -\frac{(b-a)h^4}{180}f^{(4)}(\xi) \end{aligned}
Simpson 3/8 (n a multiple of 3): weights 1, 3, 3, 2, 3, 3, 2, …, 1 times 3h/8S3/8=3h8[y0+3(y1+y2+y4+y5+⋯ )+2(y3+y6+⋯ )+yn]S_{3/8} = \frac{3h}{8}\big[y_0 + 3(y_1 + y_2 + y_4 + y_5 + \cdots) + 2(y_3 + y_6 + \cdots) + y_n\big]
  • Trapezoidal: halve hh, error ÷4\div 4; exact for straight lines
  • Simpson: halve hh, error ÷16\div 16; exact for cubics
  • Same number of function values, much better accuracy: Simpson is the default for smooth data
Md Mahbub E Noor
Lecture 20: Numerical Integration
University of Barishal
6 / 22
The method, step by step

Trapezoidal and Simpson's rules, step by step

  1. Choose nn (even for Simpson 1/3) and h=(b−a)/nh = (b - a)/n
  2. Tabulate xi=a+ihx_i = a + ih and yi=f(xi)y_i = f(x_i) for i=0,…,ni = 0, \dots, n
  3. Trapezoidal: sum = ends + 2 × inside; multiply by h/2h/2
  4. Simpson: sum = ends + 4 × odd + 2 × even (inside); multiply by h/3h/3
  5. Check: repeat with nn doubled; the change estimates the error
Pseudocode
input f, a, b, n   # n even
h ← (b − a) / n
T ← (f(a) + f(b)) / 2
S ← f(a) + f(b)
for i = 1 … n−1:
    y ← f(a + i·h)
    T ← T + y
    S ← S + (4 if i odd else 2)·y
return h·T,  h·S / 3

One loop, n+1n + 1 function values, both rules at once.

Md Mahbub E Noor
Lecture 20: Numerical Integration
University of Barishal
7 / 22

Explore: panels under the curve

  • Shaded shapes: what each rule actually integrates (trapezoids or parabolas)
  • Slide nn: watch the trapezoidal error fall by 4 per doubling, Simpson by 16
  • Simpson 3/8 needs nn divisible by 3
  • The table compares all rules with the same points
Md Mahbub E Noor
Lecture 20: Numerical Integration
University of Barishal
8 / 22

Problem 1: ∫06dx1+x2\int_0^6 \frac{dx}{1+x^2} with h=1h = 1

Evaluate ∫06dx1+x2\displaystyle\int_0^6 \frac{dx}{1+x^2} by the trapezoidal rule and both Simpson rules with h=1h = 1. Compare with tan⁡−16=1.405648\tan^{-1}6 = 1.405648.
1
Table: y=1, 0.5, 0.2, 0.1, 0.058824, 0.038462, 0.027027y = 1,\ 0.5,\ 0.2,\ 0.1,\ 0.058824,\ 0.038462,\ 0.027027
2
Trapezoidal: 12[1+0.027027+2(0.5+0.2+0.1+0.058824+0.038462)]=1.410799\tfrac12\big[1 + 0.027027 + 2(0.5 + 0.2 + 0.1 + 0.058824 + 0.038462)\big] = 1.410799
3
Simpson 1/3: 13[1.027027+4(0.5+0.1+0.038462)+2(0.2+0.058824)]=1.366173\tfrac13\big[1.027027 + 4(0.5 + 0.1 + 0.038462) + 2(0.2 + 0.058824)\big] = 1.366173
4
Simpson 3/8: 38[1.027027+3(0.5+0.2+0.058824+0.038462)+2(0.1)]=1.357081\tfrac38\big[1.027027 + 3(0.5 + 0.2 + 0.058824 + 0.038462) + 2(0.1)\big] = 1.357081
Errors: trapezoidal 0.005, Simpson 1/3 0.039, Simpson 3/8 0.049. With h=1h = 1 the curve bends sharply near 0, so the parabolas fit badly: Simpson is not automatically better when hh is too large.
Md Mahbub E Noor
Lecture 20: Numerical Integration
University of Barishal
9 / 22

Problem 2: the error falls as predicted

Compute ∫01e−x2dx=0.746824133\int_0^1 e^{-x^2}dx = 0.746824133 with n=2,4,8,16n = 2, 4, 8, 16. How does the error change when nn doubles?
nnTrapezoidalerrorratioSimpsonerrorratio
20.7313701.55×10−21.55\times10^{-2}—0.747180433.56×10−43.56\times10^{-4}—
40.7429843.84×10−33.84\times10^{-3}4.00.746855383.12×10−53.12\times10^{-5}11.4
80.7458669.59×10−49.59\times10^{-4}4.00.746826121.99×10−61.99\times10^{-6}15.7
160.7465852.40×10−42.40\times10^{-4}4.00.746824261.25×10−71.25\times10^{-7}15.9

Ratios of 4 and 16 confirm the orders h2h^2 and h4h^4. With 17 function values Simpson already gives 7 correct digits.

Md Mahbub E Noor
Lecture 20: Numerical Integration
University of Barishal
10 / 22

Romberg integration: Richardson applied to the trapezoidal rule

Halve h repeatedly, then remove the h², h⁴, … error termsRk,0=T(h/2k)Rk,j=Rk,j−1+Rk,j−1−Rk−1,j−14j−1\begin{aligned} R_{k,0} &= T(h/2^k) \\ R_{k,j} &= R_{k,j-1} + \frac{R_{k,j-1} - R_{k-1,j-1}}{4^j - 1} \end{aligned}
  • Column 1 is exactly Simpson's rule
  • R3,3R_{3,3} has 7 correct digits from only 9 function values
  • The same trick as Lecture 19, now in the stable direction
kkRk,0R_{k,0} (trap)Rk,1R_{k,1}Rk,2R_{k,2}Rk,3R_{k,3}
00.683940
10.7313700.747180
20.7429840.7468550.746834
30.7458660.7468260.7468240.7468240
Md Mahbub E Noor
Lecture 20: Numerical Integration
University of Barishal
11 / 22

Gauss quadrature: choose the points as well as the weights

Two-point Gauss–Legendre on [−1, 1]: exact for every cubic∫−11f(t) dt≈f(−13)+f(13)\int_{-1}^{1} f(t)\,dt \approx f\Big(-\frac{1}{\sqrt3}\Big) + f\Big(\frac{1}{\sqrt3}\Big)
Map [a, b] to [−1, 1]∫abf(x) dx=b−a2∫−11f(a+b2+b−a2t)dt\int_a^b f(x)\,dx = \frac{b-a}{2}\int_{-1}^{1} f\Big(\frac{a+b}{2} + \frac{b-a}{2}t\Big)dt
  • Newton–Cotes fixes equally spaced points; Gauss also chooses where to sample
  • nn points are exact for polynomials of degree 2n−12n - 1 (Simpson: 3 points, degree 3; Gauss: 2 points, degree 3)
  • ∫01e−x2dx\int_0^1 e^{-x^2}dx: 2-point Gauss 0.746595 (error 2.3×10−42.3\times10^{-4}); 3-point (weights 59,89,59\tfrac59, \tfrac89, \tfrac59) error 9.5×10−69.5\times10^{-6}
  • Needs ff at irrational points: fine for formulas, impossible for tables
Md Mahbub E Noor
Lecture 20: Numerical Integration
University of Barishal
12 / 22
Classical vs numerical

Same area, two approaches

Calculus

Classical (antiderivative)

  • ∫01xexdx\int_0^1 xe^x dx: integrate by parts
  • =[xex−ex]01=1= [xe^x - e^x]_0^1 = 1
  • Exact, if an antiderivative can be found
vs
Algorithm

Numerical (Simpson, n=4n = 4)

  • Five values of xexxe^x
  • S=1.000169S = 1.000169 (error 1.7×10−41.7\times10^{-4})
  • Same code for any integrand or table
When the classical approach cannot work

∫01e−x2dx\int_0^1 e^{-x^2}dx and the speedometer table

e−x2e^{-x^2} has no elementary antiderivative (Liouville's theorem): every normal-distribution table and every erf routine is computed numerically. The speed readings have no formula at all. Simpson: 0.7468554 with 5 values; distance 1620 m from 7 readings.

Md Mahbub E Noor
Lecture 20: Numerical Integration
University of Barishal
13 / 22

Trapezoidal, Simpson and Romberg in Python and C

  • One loop computes both composite rules
  • Romberg table built row by row
  • Applied to e−x2e^{-x^2} on [0,1][0, 1]
  • Lab: integrate the speedometer table
Md Mahbub E Noor
Lecture 20: Numerical Integration
University of Barishal
14 / 22
When quadrature goes wrong

Limitations shown failing

  • Aliasing: sin⁡2(20x)\sin^2(20x) on [0,π][0, \pi] with n=10n = 10 or 2020 samples only its zeros: every rule says 0 (true π/2\pi/2)
  • Singular derivatives: x\sqrt{x} near 0; Simpson's error falls only like h1.5h^{1.5}, not h4h^4
  • Simpson needs even nn; 3/8 needs a multiple of 3
  • Jumps inside a panel break the polynomial assumption: split the integral at the jump
Md Mahbub E Noor
Lecture 20: Numerical Integration
University of Barishal
15 / 22
Practice

Practice problem

Problem 3"Evaluate ∫01dx1+x\displaystyle\int_0^1 \frac{dx}{1+x} using the trapezoidal rule, Simpson's 1/3 rule and Simpson's 3/8 rule with h=1/6h = 1/6. Hence obtain an approximate value of ln⁡2\ln 2."
How to write a complete answer
  • Tabulate y=1/(1+x)y = 1/(1+x) at x=0,16,…,1x = 0, \tfrac16, \dots, 1: 1, 0.857143, 0.75, 0.666667, 0.6, 0.545455, 0.5
  • Trapezoidal: 0.694877; Simpson 1/3: 0.693170; Simpson 3/8: 0.693195
  • Since ∫01dx1+x=ln⁡2=0.693147\int_0^1 \frac{dx}{1+x} = \ln 2 = 0.693147, state the errors: 1.7×10−31.7\times10^{-3}, 2.3×10−52.3\times10^{-5}, 4.8×10−54.8\times10^{-5}
  • Write each rule's formula before substituting; show the grouped sums (ends, odds, evens)
Md Mahbub E Noor
Lecture 20: Numerical Integration
University of Barishal
16 / 22
When to use it

Numerical integration: method profile

Best case
  • Smooth integrands: Simpson gains a factor 16 per halving of hh
  • Tabulated data at equal steps: trapezoidal or Simpson directly
  • Romberg or Gauss: many digits from few function values
Worst case
  • Singularities or infinite derivatives at an end point
  • Fast oscillations with too few samples (aliasing)
  • Jumps inside a panel
Best suited for
  • Areas, distances, energies and probabilities from measurements
  • Statistics and ML: AUC, expectations, normalising constants
  • Inside ODE methods (Heun = trapezoid, Milne = Simpson: Lectures 21–22)
Avoid when
  • Simpson with odd nn (use 3/8 for one block)
  • Very high-order Newton–Cotes rules (negative weights, Runge-like trouble)
  • Integrals over many dimensions (use Monte Carlo)

Rule of thumb: tables → Simpson (trapezoidal for unequal or rough data); formulas → Gauss or an adaptive library routine; always compare two step sizes.

Md Mahbub E Noor
Lecture 20: Numerical Integration
University of Barishal
17 / 22

Comparative analysis: ∫01e−x2dx\int_0^1 e^{-x^2}dx with about 5 function values

RuleFunction valuesExact for degreeError orderValueError
Trapezoidal, n=4n = 451h2h^20.7429843.8×10−33.8\times10^{-3}
Simpson 1/3, n=4n = 453h4h^40.7468553.1×10−53.1\times10^{-5}
Gauss 2-point, 2 panels43 per panelh4h^40.7468032.1×10−52.1\times10^{-5}
Gauss 3-point, 1 panel35—0.7468159.5×10−69.5\times10^{-6}
Romberg R2,2R_{2,2}55h6h^60.7468349.6×10−69.6\times10^{-6}

Same small budget of function values, accuracy differing by a factor of about 400. Gauss needs ff at irrational points; the others work on tables.

Md Mahbub E Noor
Lecture 20: Numerical Integration
University of Barishal
18 / 22
Roadmap

End of Module 5: from areas to differential equations

Integration rules are accurate and stable. Their next job is to step differential equations forward in time:

What it does well (keep this)
  • Trapezoidal O(h2)O(h^2), Simpson O(h4)O(h^4)
  • Error estimate by halving hh
  • Romberg and Gauss for high accuracy
  • Stable: errors in data are averaged, not amplified
LimitationHow it is overcomeWhere
y′=f(x,y)y' = f(x, y): the integrand depends on the unknown
→
Euler (rectangle) and Heun (trapezoid)
Lecture 21 →
Higher accuracy for ODEs
→
Runge–Kutta (Simpson-like weights); Milne uses Simpson
Lecture 22 →
Integrals of noisy data with a model
→
Fit first, then integrate: least squares
Lecture 23 →
Long simulations must not blow up
→
Stability of ODE methods; stiff solvers
Lecture 22 →
Md Mahbub E Noor
Lecture 20: Numerical Integration
University of Barishal
19 / 22
Applications

Where this method runs in the real world

Click a card to see the integration rule at work.

Md Mahbub E Noor
Lecture 20: Numerical Integration
University of Barishal
20 / 22

Checkpoint

Checkpoint · think, then discuss with your neighbourSimpson's 1/3 rule with h=1h = 1 and y=1,3,9y = 1, 3, 9. What is the value?
A223\tfrac{22}{3}
B88
C163\tfrac{16}{3}
D1313
A. 13[1+4(3)+9]=223≈7.333\tfrac13[1 + 4(3) + 9] = \tfrac{22}{3} \approx 7.333. (Trapezoidal gives 12[1+2(3)+9]=8\tfrac12[1 + 2(3) + 9] = 8.)
Md Mahbub E Noor
Lecture 20: Numerical Integration
University of Barishal
21 / 22

Summary

Key points
  • Newton–Cotes: integrate the interpolating polynomial; n=1,2,3n = 1, 2, 3 give trapezoidal, Simpson 1/3, 3/8
  • Trapezoidal error −(b−a)h212f′′-\frac{(b-a)h^2}{12}f''; Simpson −(b−a)h4180f(4)-\frac{(b-a)h^4}{180}f^{(4)} (exact for cubics)
  • Romberg: Richardson on the trapezoidal rule; Gauss: optimal points, nn points exact to degree 2n−12n - 1
  • Check by halving hh; beware singularities, aliasing, jumps and odd nn
Opening problem: solved

Trapezoidal: 1610 m; Simpson: 1620 m. They differ by 0.6%, so the distance is about 1.6 km to within about 10 m.

Homework
  1. Evaluate ∫0πsin⁡x dx\int_0^{\pi} \sin x\,dx with n=6n = 6 by the trapezoidal and Simpson rules; compare with 2.
  2. Derive Simpson's 1/3 rule from the general quadrature formula with n=2n = 2.
  3. Integrate the speedometer table with Simpson's 3/8 rule (two blocks). Compare with 1610 m and 1620 m.
  4. Build the Romberg table for ∫0141+x2dx\int_0^1 \frac{4}{1+x^2}dx up to R2,2R_{2,2}. Which famous number appears?
Next lecture

ODEs I: Taylor, Picard and Euler. When the rate of change depends on the unknown itself: stepping y′=f(x,y)y' = f(x, y) forward.

Md Mahbub E Noor
Lecture 20: Numerical Integration
University of Barishal
22 / 22