Implementing numerical integration from scratch — trapezoidal, Simpson's: Visual walkthrough

5.4.24 · D2 · Coding › Scientific Computing (Python) › Implementing numerical integration from scratch — trapezoida

Hum sirf yahi assume karenge ki tum:

ek graph padh sako (ek curve, ek $x$ -axis left-to-right jaati hui, ek $y$ -axis upar jaati hui),
ek rectangle ka area paane ke liye width ko height se multiply kar sako.

Baaki sab kuch yahan build kiya gaya hai.

Step 1 — "Area under a curve" ka matlab kya hai

KYA. Hum goal naam dete hain: figure mein shaded area.

KYUN. Bahut saare curves ke liye is area ka koi clean formula nahi hota (koi nice antiderivative nahi), ya humhare paas formula nahi hota balki sirf data points hote hain. Isliye hum ise kabhi exactly compute nahi kar paate — hum us region ko un shapes se cover karenge jinका area hum pehle se jaante hain aur unhe add karenge.

PICTURE. Green region exactly woh area hai jisko hum dhund rahe hain. Notice karo ki uski curved lid hai — woh curve hi poori difficulty hai.

Figure — Implementing numerical integration from scratch — trapezoidal, Simpson's

Step 2 — Region ko equal strips mein chop karo

KYA. Hum $[a, b]$ ko $n$ vertical strips mein kaatते hain jinkaa equal width $h$ hai.

KYUN. Equal width ka matlab hai ki har strip ko same tarike se treat kiya jaata hai — ek strip ke liye ek formula, phir bas add karo. Isse shared-edge bookkeeping (Step 5 mein aane wali) bhi clean rehti hai.

PICTURE. Dotted verticals har node se neeche girte hain; base $n$ equal chunks mein split hai jinkaa width $h$ hai.

Step 3 — EK strip ko straight roof (trapezoid) se cover karo

Ek trapezoid ka area hai do parallel heights ka average, times base: $A_{i} = average height \frac{f ( x _{i} ) + f ( x _{i + 1} )}{2} \times width h .$ Term by term: $f (x_{i})$ left height hai, $f (x_{i + 1})$ right height hai; add karke aur half karne se slanted lid ki middle height milti hai; width $h$ se multiply karne par area milta hai — exactly ek rectangle jis ki height woh average hai.

KYA. Ek strip ka area, trapezoid se approximate kiya gaya.

KYUN yeh shape aur rectangle nahi? Ek flat rectangle sirf ek height use karta hai aur curve ki tilt ignore karta hai; slanted lid dono edge heights use karta hai, isliye curve ki slope follow karta hai. Yeh same effort mein strictly better guess hai.

PICTURE. Teal trapezoid; uski slanted lid dono node dots par curve ko touch karti hai. Lid aur curve ke beech chhota green gap is strip par hone wali error hai.

Step 4 — Saare trapezoids add karo → composite Trapezoidal rule

Hum har strip ka area add karte hain aur dekhte hain ki har node kya contribute karta hai: $\sum_{i = 0}^{n - 1} \frac{f ( x _{i} ) + f ( x _{i + 1} )}{2} h .$

KYA / KYUN. Ek single interior node dekho, maano $x_{1}$ . Yeh strip 0 ka right edge aur strip 1 ka left edge dono hai, isliye iski height $f (x_{1})$ do trapezoids mein aati hai — yeh twice count hoti hai. Dono far ends, $x_{0}$ aur $x_{n}$ , sirf ek trapezoid ke hain — once count hote hain. Like terms collect karne par:

PICTURE. Poora region trapezoids se tile kiya gaya. End nodes ek baar glow karte hain; har shared interior node do baar glow karta hai — literally weight-2.

Recall Trapezoid error

$O (h^{2})$ kyun hota hai Straight lid ek straight curve ke liye exact hai (degree-1). Ek bending curve ke liye per strip ka leftover gap $h^{3}$ ki tarah shrink karta hai, aur roughly $1/ h$ strips hoti hain, isliye total error $h^{2}$ ki tarah scale karta hai. $h$ aadha karo → error ~4× kam ho jaati hai. Dekho Taylor series.

Step 5 — Woh bend jise trapezoids pakad nahi sakte

KYA. Hum ek single hump par dikhate hain ki straight lid ka woh area sliver hamesha table par rehta hai.

KYUN aage parabola? Line ke baad next-simplest curve ek parabola hai — $α + β x + γ x^{2}$ form ki shape. Woh extra $γ x^{2}$ term exactly ek bend hai. Ek parabola teen points se fix hoti hai, isliye hum apna view ek strip se ek pair of strips (teen nodes) tak badhaate hain. Yahi wajah hai ki Simpson's rule pairs mein kaam karta hai — aur isliye $n$ even hona zaroori hai.

PICTURE. Straight lid (dashed teal) vs true curve; orange sliver ek bend par unavoidable trapezoid error hai.

Step 6 — Teen points se parabola fit karo, usse exactly integrate karo

Ek pair of strips lo: nodes $x_{0}, x_{1}, x_{2}$ , equally $h$ se spaced. Clever move: origin slide karo taaki middle node $0$ par baith jaaye. Tab teen nodes $- h, 0, + h$ par hain. Yeh symmetry hi algebra ko collapse kati hai.

Parabola ko $p (x) = α + β x + γ x^{2}$ kaho:

$α$ = centre par iski height ( $x = 0$ par value),
$β$ = iski tilt (linear part),
$γ$ = iski bend (curviness).

Ise pair ke across integrate karo: $\int_{- h}^{h} (α + β x + γ x^{2}) d x = 2 α h + \frac{2}{3} γ h^{3} .$ Term by term: $α$ integrate hokar $α \cdot (length 2 h) = 2 α h$ deta hai. Tilt term $β x$ integrate hokar zero deta hai — kyunki symmetric span $[- h, h]$ par yeh right mein utna add karta hai jitna left mein subtract karta hai. Bend $γ x^{2}$ integrate hokar $γ \cdot \frac{2 h ^{3}}{3}$ deta hai.

PICTURE. Shifted axis $- h, 0, h$ par teen dots se guzarti parabola; $β x$ tilt equal-and-opposite red/blue lobes ki tarah shaded jo cancel ho jaati hain.

Step 7 — $α, γ$ ko data mein convert karo → 1-4-1 stencil appear hota hai

Teen node positions ko $p$ mein plug karo: $f (- h) = α - β h + γ h^{2}, f (0) = α, f (h) = α + β h + γ h^{2} .$

Do facts turant nikaalte hain:

Middle equation kehta hai $α = f (0)$ — centre height hai hi middle data point.
Outer dono add karo: $β h$ cancel ho jaata hai, deta hai $f (- h) + f (h) = 2 α + 2 γ h^{2}$ , isliye $γ h^{2} = \frac{f ( - h ) + f ( h )}{2} - f (0) .$

Dono ko area $2 α h + \frac{2}{3} γ h^{3} = 2 α h + \frac{2}{3} h (γ h^{2})$ mein substitute karo:

=\frac{h}{3}\big[f(-h)+4f(0)+f(h)\big].$$ **KYA.** Fitted parabola ke neeche exact area, purely teen heights mein likha gaya. **KYUN numbers 1, 4, 1?** Common denominator $3$ par coefficients collect karo. Do **edge** heights har ek ek single $\frac{h}{3}$ ke saath bachte hain → weight **1**. **Centre** height $\alpha$ term ke bade $2$ ko *plus* correction ke saath pick up karta hai, netting $4$. Yahi celebrated stencil hai — memorised nahi, balki *manufactured*. > [!formula] Simpson's rule (ek pair) > $$\int_{x_0}^{x_2}f\,dx\approx\frac{h}{3}\big[\,\underbrace{f(x_0)}_{\times 1}+\underbrace{4f(x_1)}_{\times 4}+\underbrace{f(x_2)}_{\times 1}\,\big]$$ **PICTURE.** Teen nodes apne falling weights $1,\,4,\,1$ ke saath tagged hain; centre node sabse zyada blazes karta hai kyunki parabola usse sabse zyada lean karta hai. ![[deepdives/dd-coding-5.4.24-d2-s07.png]] --- ## Step 8 — Pairs chain karo → composite Simpson's rule Ab $[a,b]$ par yeh parabolic pairs side by side rakh do (iske liye $n$ **even** hona chahiye, taaki strips bina kisi leftover ke pair ho sakein). Track karo ki har node kitni baar use hota hai: - **Ends** $x_0,x_n$: ek pair ke hain → weight **1**. - **Odd-index** nodes ($x_1,x_3,\dots$): har ek apne pair ka **centre** hai → weight **4**. - **Even interior** nodes ($x_2,x_4,\dots$): har ek woh **shared boundary** hai jahan do pairs milte hain → weight $1+1=$ **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]$$ > Weight pattern: $1,4,2,4,2,\dots,4,1$. Error $\sim O(h^4)$ — $h$ aadha karo, error ~16× kam ho jaati hai. > [!mistake] Odd $n$ Simpson ko tod deta hai > Ek dangling strip ka koi partner nahi hota parabola form karne ke liye. 1-4-1 stencil strips ko pairs mein khaata hai. **Fix:** $n$ even require karo, ya last lone strip ko trapezoid se treat karo. **PICTURE.** Kaafi parabolic arcs curve ko tile kar rahe hain; nodes unke weight $1/4/2$ ke hisaab se coloured hain — alternating $4,2,4,2$ rhythm visible hai. ![[deepdives/dd-coding-5.4.24-d2-s08.png]] --- ## Step 9 — Parent ke example par Sanity check > [!example] $\int_0^1 x^2\,dx=\tfrac13$, $n=2$, $h=0.5$ > Nodes $0,0.5,1$ heights dete hain $0,\,0.25,\,1$. > $$S_2=\frac{0.5}{3}\big[\underbrace{0}_{\times1}+\underbrace{4(0.25)}_{\times4}+\underbrace{1}_{\times1}\big]=\frac{0.5}{3}\cdot 2=0.3333\ldots$$ > **Exact.** Kyunki $x^2$ degree 2 hai aur hamari lid *hai hi* ek parabola — yeh curve ko sirf approximate nahi karti, yeh **equal** karti hai. (Bonus: Simpson degree-3 cubics ke liye bhi exact hai — extra odd term cancel ho jaata hai jaise $\beta x$ Step 6 mein hua tha.) > Compare karo $T_2=\frac{0.5}{2}[0+1+2(0.25)]=0.375$ se — straight lids convex curve ko overshoot karti hain, exactly Step-5 wali failure. --- ## Ek-picture summary ![[deepdives/dd-coding-5.4.24-d2-s09.png]] Ek frame, poori derivation: exact green area → straight trapezoid lids (weights $1,2,\dots,2,1$, error $O(h^2)$) → bendy parabola lids over pairs (weights $1,4,2,4,\dots,4,1$, error $O(h^4)$). Same nodes, smarter lids, smaller gaps. > [!recall]- Feynman retelling — poora walk plain words mein > Tum ek wiggly hill ka area chahte ho lekin wiggle directly measure nahi kar sakte. Toh tum ise un shapes se cover karte ho jo tum *jaante* ho. Pehle base ko equal-width strips mein kaato. Har thin strip par hill almost straight hai, isliye tum ise ek slanted-roof box — trapezoid — se cap karte ho aur unhe add karte ho. Poori hill ke edges sirf ek box ko touch karti hain (ek baar count), lekin har inside cut-line do boxes ke beech share hoti hai (do baar count) — yahi $1,2,2,\dots,2,1$ pattern hai. Trapezoids quick hain lekin unke flat roofs bend nahi kar sakte, isliye ek hump par woh hamesha ek chhota sliver miss karte hain. Fix: ek aisi roof use karo jo bend kare. Sabse simple bending shape parabola hai, aur ek parabola ko pin down karne ke liye teen points chahiye — isliye tum do strips ek saath kaam karte ho. Agar tum cleverly middle point ko zero par rakh do, to parabola ki tilt khud hi cancel ho jaati hai jab tum iska area nikalte ho, sirf middle height aur bend reh jaate hain. Unhe apni teen data heights mein wapas convert karne par, numbers $1,4,1$ khud nikaal aate hain — middle chaar baar count hota hai kyunki parabola sabse zyada usse lean karta hai. Inhe parabola-pairs ko hill par chain karo aur shared boundaries $1,4,2,4,2,\dots,4,1$ ka rhythm banati hain. Kyunki jab curve ek parabola (ya cubic bhi) hai toh parabola actually *hai hi* woh curve, isliye Simpson unhe exactly nail karta hai — yahi wajah hai ki woh same amount of work ke liye trapezoids ko beat karta hai. > [!mnemonic] Poora page ek saath mein > **Straight roofs = $1,2,\dots,2,1$. Bendy roofs, in pairs = $1,4,2,4,\dots,4,1$. Even $n$ ya koi partner nahi.** ## Connections - [[Riemann sums]] — woh limit-of-sums jise dono rules secretly refine kar rahe hain. - [[Newton-Cotes formulas]] — trapezoid aur Simpson is family ke pehle do members hain. - [[Polynomial interpolation]] — Step 6–7 *hai hi* interpolating parabola fit karna. - [[Taylor series]] — jahan se $O(h^2)$ aur $O(h^4)$ error orders aate hain. - [[Richardson extrapolation]] — trapezoid results blend karo taaki Simpson *manufacture* ho, phir Romberg. - [[scipy.integrate.quad]] — is poore idea ka adaptive, production version.