5.4.24Scientific Computing (Python)
Implementing numerical integration from scratch — trapezoidal, Simpson's
1,592 words7 min readdifficulty · medium
WHAT are we computing?
Trapezoidal rule — derive it
HOW (single slice): approximate by the line through and . Area of that trapezoid: Why this step? The trapezoid's two parallel sides are the function heights; average them, multiply by base .
Sum all slices:
= \frac{h}{2}\Big[f(x_0)+2f(x_1)+\dots+2f(x_{n-1})+f(x_n)\Big].$$ *Why the factor 2?* Each interior node is shared by two adjacent trapezoids, so it's counted twice; endpoints only once. > [!formula] Composite Trapezoidal Rule > $$T_n=\frac{h}{2}\left[f(x_0)+f(x_n)+2\sum_{i=1}^{n-1}f(x_i)\right],\quad h=\frac{b-a}{n}$$ > Error $\sim O(h^2)$: exact for ==straight lines (degree-1 polynomials)==. ```python def trapezoid(f, a, b, n): h = (b - a) / n s = 0.5 * (f(a) + f(b)) # endpoints get weight 1/2 for i in range(1, n): # interior nodes weight 1 s += f(a + i * h) return s * h ``` --- ## Simpson's rule — derive it > [!intuition] WHY a parabola? > A straight line can't bend; a curve does. Fit a **parabola** through *three* points (two slices) and it hugs the curve far better. Cost: $n$ must be **even** (we work in pairs). **HOW (one pair of slices, nodes $x_0,x_1,x_2$, spacing $h$):** Fit $p(x)=\alpha + \beta x + \gamma x^2$ through the three points and integrate exactly. Shift origin so $x_1=0$ (nodes at $-h,0,h$): $$\int_{-h}^{h}(\alpha+\beta x+\gamma x^2)\,dx = 2\alpha h + \tfrac{2}{3}\gamma h^3.$$ *Why this step?* Odd term $\beta x$ integrates to 0 over symmetric limits; only $\alpha$ and $\gamma$ survive. Now express $\alpha,\gamma$ via the data: $$f(-h)=\alpha-\beta h+\gamma h^2,\quad f(0)=\alpha,\quad f(h)=\alpha+\beta h+\gamma h^2.$$ So $\alpha=f(0)$ and $f(-h)+f(h)=2\alpha+2\gamma h^2 \Rightarrow \gamma h^2 = \tfrac{f(-h)+f(h)}{2}-f(0)$. Substitute: $$\int = 2 f(0) h + \tfrac{2}{3}h\big(\tfrac{f(-h)+f(h)}{2}-f(0)\big) = \frac{h}{3}\big[f(-h)+4f(0)+f(h)\big].$$ *Why this step?* Collect terms over common denominator $3$ → the famous **1-4-1** pattern. > [!formula] Simpson's rule (one pair) > $$\int_{x_0}^{x_2} f\,dx \approx \frac{h}{3}\big[f(x_0)+4f(x_1)+f(x_2)\big]$$ **Composite:** chain pairs across $[a,b]$. Odd-index nodes are pair-centres (weight 4); even interior nodes are shared between adjacent pairs (weight 2): > [!formula] Composite Simpson's Rule ($n$ even) > $$S_n=\frac{h}{3}\Big[f(x_0)+f(x_n)+4\!\!\sum_{i\,\text{odd}}\!\!f(x_i)+2\!\!\sum_{i\,\text{even, interior}}\!\!f(x_i)\Big]$$ > Error $\sim O(h^4)$: exact for ==cubics== (degree-3!), a free bonus. ```python def simpson(f, a, b, n): if n % 2 == 1: raise ValueError("Simpson needs even n") h = (b - a) / n s = f(a) + f(b) for i in range(1, n): s += (4 if i % 2 == 1 else 2) * f(a + i * h) return s * h / 3 ``` ![[5.4.24-Implementing-numerical-integration-from-scratch-—-trapezoidal,-Simpson's.png]] --- ## Worked examples > [!example] $\int_0^1 x^2\,dx=\tfrac13$ with $n=2$ (h=0.5) > Nodes: $x_0=0,\,x_1=0.5,\,x_2=1$; $f=0,\,0.25,\,1$. > **Trapezoid:** $T_2=\frac{0.5}{2}[0+1+2(0.25)] = 0.25\cdot1.5=0.375$. *Why?* lines overshoot a convex curve → too big. > **Simpson:** $S_2=\frac{0.5}{3}[0+4(0.25)+1]=\frac{0.5}{3}(2)=0.3333\ldots$ → **exact**. > *Why exact?* $x^2$ is degree 2, and Simpson integrates parabolas exactly. > [!example] $\int_0^\pi \sin x\,dx = 2$ with $n=4$ ($h=\pi/4$) > $f$ at nodes $0,\,\tfrac{\pi}{4},\,\tfrac{\pi}{2},\,\tfrac{3\pi}{4},\,\pi$ = $0,\,0.7071,\,1,\,0.7071,\,0$. > **Simpson:** $\frac{h}{3}[0 + 4(0.7071) + 2(1) + 4(0.7071) + 0]=\frac{\pi/4}{3}(7.657)=2.0046$. *Why so close?* $O(h^4)$ error shrinks fast even with 4 slices. > **Trapezoid:** $\frac{h}{2}[0+0+2(0.7071+1+0.7071)]=1.896$ — visibly worse for the same work. --- ## Common mistakes > [!mistake] "Simpson works for any $n$." > *Why it feels right:* it's just a weighted sum, why care about parity? **Fix:** the 1-4-1 stencil consumes nodes in *pairs of slices*. An odd $n$ leaves one dangling slice with no parabola. Always require $n$ even (or handle the last slice separately). > [!mistake] Forgetting endpoints get a different weight. > *Why it feels right:* "just sum $f$ at all nodes and scale." **Fix:** endpoints belong to only one trapezoid/parabola, so trapezoid weights them $\tfrac12$ (vs 1 interior), Simpson weights them 1 (vs 4 or 2). Mis-weighting endpoints is the #1 silent bug. > [!mistake] Confusing $n$ (subintervals) with the number of points. > Points $=n+1$. *Fix:* loop interior nodes `range(1, n)`, not `range(1, n+1)`. --- ## #flashcards/coding Trapezoidal rule formula for $n$ subintervals ::: $T_n=\frac{h}{2}[f(x_0)+f(x_n)+2\sum_{i=1}^{n-1}f(x_i)]$ Why does each interior trapezoid node get weight 2? ::: It is shared by two adjacent trapezoids, so counted twice. Simpson's single-pair stencil and its name ::: $\frac{h}{3}[f_0+4f_1+f_2]$ — the 1-4-1 pattern. Why must $n$ be even for composite Simpson? ::: Each parabola spans a pair of subintervals; pairs need an even count. Error order of trapezoid vs Simpson ::: Trapezoid $O(h^2)$, Simpson $O(h^4)$. Highest polynomial degree each rule integrates exactly ::: Trapezoid: degree 1; Simpson: degree 3. Node spacing $h$ in terms of $a,b,n$ ::: $h=(b-a)/n$. Number of evaluation points for $n$ subintervals ::: $n+1$. In Simpson, which nodes get weight 4 vs 2? ::: Odd-index nodes → 4; even interior nodes → 2. --- > [!recall]- Feynman: explain to a 12-year-old > Imagine you want the area of a wiggly hill. You can't measure it directly, so you cover it with shapes you *know* the area of. The **trapezoid** way puts slanted-roof boxes side by side — quick but the slanted roofs miss the curve a little. **Simpson's** way uses bendy roofs (parabolas) over every pair of boxes — they hug the hill much closer, so your area guess is much better with the same number of boxes. The trick in both: count the heights, but the heights on the **edges** count differently from the ones in the **middle**. > [!mnemonic] Remember Simpson's weights > **"One, four, two... ride a parabolic boat across pairs."** > Edges = **1**, odd middles = **4** (the peak the parabola bends to), even shared middles = **2**. And **"Even n for Simpson — pairs need partners."** ## Connections - [[Riemann sums]] — the limit definition both rules approximate. - [[Newton-Cotes formulas]] — trapezoid (deg 1) & Simpson (deg 2) are the first members. - [[Polynomial interpolation]] — Simpson = integrate the interpolating parabola. - [[Taylor series]] — source of the $O(h^2)$/$O(h^4)$ error terms. - [[Richardson extrapolation]] — combine $T_n$ to get Simpson; combine Simpson for Romberg. - [[scipy.integrate.quad]] — production-grade adaptive integration. ## 🖼️ Concept Map ```mermaid flowchart TD I[Definite integral area under curve] -->|approximated as| W[Weighted sum sum wi f xi] W -->|split a,b into n slices| H[Node spacing h = b-a over n] W -->|straight-line fit| TR[Trapezoidal rule] W -->|parabola fit| SI[Simpson's rule] TR -->|one slice| TA[Area = avg height x h] TA -->|sum all slices| TC[Composite Tn formula] TC -->|interior nodes shared| TW[Weights 1,2,...,2,1 x h/2] TC -->|exact for degree-1| TE[Error O h^2] SI -->|fit 3 points| SP[p x = alpha + beta x + gamma x^2] SP -->|integrate exactly| SF[h/3 f-h + 4f0 + fh] SI -->|requires| SN[n must be even] SF -->|smaller error than| TR ``` ## 🔊 Hinglish (regional understanding) > [!intuition]- Hinglish mein samjho > Dekho, integration ka matlab hai curve ke neeche ka **area** nikalna. Lekin har function ka clean antiderivative nahi milta, ya kabhi humein sirf data points milte hain. Toh hum curve ko aisi shapes se replace karte hain jinka area hum already jaante hain. Yeh hai numerical integration ka core idea. > > **Trapezoidal rule** mein hum do consecutive points ko ek seedhi line se jodte hain, jisse ek trapezium ban jaata hai — uska area = (do heights ka average) × width $h$. Saare slices add karo, toh interior points do baar count hote hain (isliye weight 2) aur endpoints sirf ek baar (weight 1). Simple aur fast, par seedhi line curve ko thoda miss kar deti hai, isliye error $O(h^2)$ hota hai. > > **Simpson's rule** smarter hai: do slices ke upar ek **parabola** fit karta hai jo curve ko bahut better hug karti hai. Derive karne par 1-4-1 ka famous pattern aata hai (composite mein 1-4-2-4-1...). Iska error $O(h^4)$ hai — matlab same number of points par bahut zyada accurate. Ek bonus: parabola hone ke baad bhi yeh degree-3 (cubic) tak ke functions ko bilkul exact integrate kar deta hai! > > Do cheezein yaad rakhna: Simpson ke liye $n$ **even** hona chahiye (kyunki parabola pairs mein kaam karti hai), aur endpoints ka weight hamesha alag hota hai. Yeh do silent bugs sabse zyada log karte hain. Bas weights sahi lagao, aur tumhara from-scratch integrator scipy jaisa kaam karega. ![[audio/5.4.24-Implementing-numerical-integration-from-scratch-—-trapezoidal,-Simpson's.mp3]]