Numerical integration - Trapezoidal vs Simpson - Reversing the wheel of learning - High end computer science is plain Maths

1. The Core Idea

We want to approximate:

$$\int_a^b f(x)\,dx$$

by sampling the function at discrete points.

2. Trapezoidal Rule

Instead of rectangles, we use trapezoids.

${\int }_a^{b}f(x)dx\approx \frac{h}{2}[f(x_0)+2{\sum }_{i=1}^{n-1}f(x_i)+f(x_n)]$

Key properties:

• Error ~ O(h²)
• Assumes linear variation between points

👉 Think: “connect the dots with straight lines”

3. Simpson’s Rule

Now we approximate using parabolas (quadratic fit).

$\int_a^b f(x)\,dx \approx \frac{h}{3}\left[f(x_0) + 4\sum_{\text{odd}} f(x_i) + 2\sum_{\text{even}} f(x_i) + f(x_n)\right]$

Key properties:

• Error ~ O(h⁴) (much faster convergence)
• Requires even number of intervals

Think: “fit smooth curves instead of straight lines”

4. The Real Question: Which is Better?

Method	Accuracy	Cost	When it works best
Trapezoidal	Medium	Low	Rough/oscillatory data
Simpson	High	Medium	Smooth functions

5. A Perfect Test Function

Use something smooth but non-trivial:

$$f(x) = \sin(x), \quad \int_0^\pi \sin(x)\,dx = 2$$

6. The source code (julia)

using Plots

gr()

f(x) = sin(x)

# Trapezoidal Rule

function trapezoidal(f, a, b, n)

h = (b - a) / n

s = 0.5 * (f(a) + f(b))

for i in 1:n-1

s += f(a + i*h)

end

return h * s

end

# Simpson Rule

function simpson(f, a, b, n)

h = (b - a) / n

s = f(a) + f(b)

for i in 1:n-1

x = a + i*h

s += (i % 2 == 0 ? 2 : 4) * f(x)

end

return (h/3) * s

end

# Setup

a, b = 0, pi

exact = 2.0

#exact = exp(pi) - 1

ns = [4, 8, 16, 32, 64, 128, 256]

h_vals = Float64[]

trap_err = Float64[]

simp_err = Float64[]

for n in ns

h = (b - a) / n

push!(h_vals, h)

t = trapezoidal(f, a, b, n)

s = simpson(f, a, b, n)

push!(trap_err, abs(t - exact))

push!(simp_err, abs(s - exact))

end

# Plot

p1 = plot(h_vals, trap_err,

xscale = :log10, yscale = :log10,

marker = :circle,

label = "Trapezoidal",

xlabel = "Step size (h)",

ylabel = "Error",

title = "Error vs Step Size")

# Plot 2: Simpson

p2 = plot(h_vals, simp_err,

xscale = :log10, yscale = :log10,

marker = :square,

title = "Simpson Error",

xlabel = "Step size (h)", ylabel = "Error",

label = "Simpson")

# Combine side-by-side

p3 = plot(p1, p2,

layout = (1, 2), # 1 row, 2 columns

size = (1200, 400)) # wide figure

display(p3)

wait()

9. Deep Insight (Important)

This experiment teaches:

• Numerical methods are about approximating shape
• Higher-order ≠ always better
• Smoothness of function matters
• There is a trade-off between truncation error and floating-point error