5.4.24 · HinglishScientific Computing (Python)

Implementing numerical integration from scratch — trapezoidal, Simpson's

1,408 words6 min readRead in English

5.4.24 · Coding › Scientific Computing (Python)


WHAT are we computing?


Trapezoidal rule — derive it

HOW (single slice): ko us line se approximate karo jo aur se guzarti hai. Us trapezoid ka area: Yeh step kyun? Trapezoid ke do parallel sides function heights hain; unhe average karo, base se multiply karo.

Saare slices ka sum:

= \frac{h}{2}\Big[f(x_0)+2f(x_1)+\dots+2f(x_{n-1})+f(x_n)\Big].$$ *Factor 2 kyun?* Har interior node do adjacent trapezoids ke beech shared hota hai, isliye usse do baar count kiya jaata hai; endpoints sirf ek baar. > [!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)$: ==straight lines (degree-1 polynomials)== ke liye exact. ```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] Parabola kyun? > Ek straight line bend nahi kar sakti; curve karti hai. *Teen* points (do slices) se ek **parabola** fit karo aur yeh curve se kaafi behtar chipakti hai. Cost: $n$ **even** hona chahiye (hum pairs mein kaam karte hain). **HOW (ek pair of slices, nodes $x_0,x_1,x_2$, spacing $h$):** $p(x)=\alpha + \beta x + \gamma x^2$ ko teen points se fit karo aur exactly integrate karo. Origin shift karo taaki $x_1=0$ ho (nodes $-h,0,h$ par): $$\int_{-h}^{h}(\alpha+\beta x+\gamma x^2)\,dx = 2\alpha h + \tfrac{2}{3}\gamma h^3.$$ *Yeh step kyun?* Odd term $\beta x$ symmetric limits par 0 integrate hota hai; sirf $\alpha$ aur $\gamma$ bachte hain. Ab $\alpha,\gamma$ ko data ke through express karo: $$f(-h)=\alpha-\beta h+\gamma h^2,\quad f(0)=\alpha,\quad f(h)=\alpha+\beta h+\gamma h^2.$$ Toh $\alpha=f(0)$ aur $f(-h)+f(h)=2\alpha+2\gamma h^2 \Rightarrow \gamma h^2 = \tfrac{f(-h)+f(h)}{2}-f(0)$. Substitute karo: $$\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].$$ *Yeh step kyun?* Terms ko common denominator $3$ par collect karo → 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:** pairs ko $[a,b]$ par chain karo. Odd-index nodes pair-centres hain (weight 4); even interior nodes adjacent pairs ke beech shared hain (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)$: ==cubics== (degree-3!) ke liye exact — ek 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$. *Kyun?* Lines ek convex curve ko overshoot karti hain → zyada bada answer. > **Simpson:** $S_2=\frac{0.5}{3}[0+4(0.25)+1]=\frac{0.5}{3}(2)=0.3333\ldots$ → **exact**. > *Exact kyun?* $x^2$ degree 2 hai, aur Simpson parabolas ko exactly integrate karta hai. > [!example] $\int_0^\pi \sin x\,dx = 2$ with $n=4$ ($h=\pi/4$) > $f$ nodes $0,\,\tfrac{\pi}{4},\,\tfrac{\pi}{2},\,\tfrac{3\pi}{4},\,\pi$ par = $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$. *Itna close kyun?* $O(h^4)$ error sirf 4 slices mein bhi fast shrink hoti hai. > **Trapezoid:** $\frac{h}{2}[0+0+2(0.7071+1+0.7071)]=1.896$ — same kaam ke liye visibly worse. --- ## Common mistakes > [!mistake] "Simpson kisi bhi $n$ ke liye kaam karta hai." > *Kyun sahi lagta hai:* yeh toh sirf ek weighted sum hai, parity se kyun farq padta hai? **Fix:** 1-4-1 stencil nodes ko *pairs of slices* mein consume karta hai. Odd $n$ mein ek dangling slice bachti hai jiske liye koi parabola nahi. Hamesha $n$ even require karo (ya last slice ko alag handle karo). > [!mistake] Endpoints ko alag weight milna bhool jaana. > *Kyun sahi lagta hai:* "bas sabhi nodes par $f$ sum karo aur scale karo." **Fix:** Endpoints sirf ek trapezoid/parabola mein belong karte hain, isliye trapezoid unhe $\tfrac12$ weight deta hai (interior ke 1 ke mukable mein), Simpson unhe 1 weight deta hai (4 ya 2 ke mukable mein). Endpoints ko galat weight dena #1 silent bug hai. > [!mistake] $n$ (subintervals) aur points ki sankhya ko confuse karna. > Points $=n+1$. *Fix:* interior nodes par loop karo `range(1, n)`, `range(1, n+1)` nahi. --- ## #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)]$ Har interior trapezoid node ko weight 2 kyun milta hai? ::: Yeh do adjacent trapezoids ke beech shared hota hai, isliye do baar count hota hai. Simpson's single-pair stencil aur uska naam ::: $\frac{h}{3}[f_0+4f_1+f_2]$ — the 1-4-1 pattern. Composite Simpson ke liye $n$ even kyun hona chahiye? ::: Har parabola ek pair of subintervals span karti hai; pairs ke liye even count chahiye. Trapezoid vs Simpson ki error order ::: Trapezoid $O(h^2)$, Simpson $O(h^4)$. Har rule jitnay degree tak polynomial ko exactly integrate karta hai ::: Trapezoid: degree 1; Simpson: degree 3. Node spacing $h$ in terms of $a,b,n$ ::: $h=(b-a)/n$. $n$ subintervals ke liye evaluation points ki sankhya ::: $n+1$. Simpson mein kaunse nodes ko weight 4 milta hai vs 2? ::: Odd-index nodes → 4; even interior nodes → 2. --- > [!recall]- Feynman: 12-saal ke bachche ko explain karo > Socho tum ek wiggly pahaad ka area jaanna chahte ho. Tum ise directly measure nahi kar sakte, isliye tum ise un shapes se dhakko jo tumhe pata hai unka area. **Trapezoid** wala tarika — slanted roof wale boxes side by side rakhta hai — quick hai lekin slanted roofs curve ko thoda miss karti hain. **Simpson's** wala tarika — har pair of boxes par bendy roofs (parabolas) use karta hai — yeh pahaad se bahut closely chipakti hain, isliye teri area guess same number of boxes ke saath bahut behtar hoti hai. Dono mein trick: heights count karo, lekin **edges** par heights **middle** wali heights se alag count hoti hain. > [!mnemonic] Simpson's weights yaad karo > **"One, four, two... ride a parabolic boat across pairs."** > Edges = **1**, odd middles = **4** (woh peak jis par parabola bend karti hai), even shared middles = **2**. Aur **"Even n for Simpson — pairs need partners."** ## Connections - [[Riemann sums]] — woh limit definition jise dono rules approximate karte hain. - [[Newton-Cotes formulas]] — trapezoid (deg 1) aur Simpson (deg 2) pehle members hain. - [[Polynomial interpolation]] — Simpson = interpolating parabola ko integrate karo. - [[Taylor series]] — $O(h^2)$/$O(h^4)$ error terms ka source. - [[Richardson extrapolation]] — $T_n$ combine karo Simpson pane ke liye; Simpson combine karo Romberg ke liye. - [[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 ```