Implementing ODE solvers from scratch — Euler, RK4: Visual walkthrough

5.4.23 · D2 · Coding › Scientific Computing (Python) › Implementing ODE solvers from scratch — Euler, RK4

Yahan sab kuch parent note ko aur gehraai deta hai: the parent topic.

Step 1 — ODE actually hame kya deta hai

KYA. Hum curve $y (t)$ nahi jaante. Hamare paas sirf machine $f$ hai jo kisi bhi chosen point par slope report karta hai.

KYUN. Zyada tar real problems mein $y (t)$ ka koi formula nahi hota — lekin $f$ hamesha kuch aisa hota hai jo tum compute kar sako (yeh bas arithmetic hai). Toh game yeh hai: "main slope har jagah jaanta hoon" ko "main curve jaanta hoon" mein badlo.

PICTURE. Neeche, faint dashed curve sahi $y (t)$ hai (jo solver nahi dekh sakta). Ek point par hum ek chhota arrow rakhte hain — woh arrow wahi hai jo $f$ hame batata hai: yahaan chalne ki direction.

Figure — Implementing ODE solvers from scratch — Euler, RK4

$f (t, y)$ ::: woh machine jo kisi bhi point $(t, y)$ par slope $y^{'}$ return karta hai.
RK4 ko kaunsi smoothness chahiye? ::: $f$ itna smooth ho ki $y (t)$ ka $h^{4}$ term tak valid Taylor expansion ho (roughly 4 continuous derivatives).

Step 2 — Grid banana, phir slope ko step mein badalna

KYA. Hum $h$ ko zero tak jaane se mana karte hain. Hum isse chhota lekin finite rakhte hain aur equation ko backwards padhke next rung solve karte hain.

KYUN yeh tool. Hum $y$ ko ek rung aage predict karna chahte hain, aur derivative definition hi ek bridge hai "slope abhi" aur "value thodi der baad" ke beech. Limit drop karna woh ek approximation hai jisse sab shuru hota hai.

KAISE — algebra, ek move at a time. Finite- $h$ version se shuru karo aur limit drop karo: $f (t_{n}, y_{n}) \approx \frac{y _{n + 1} - y _{n}}{h}$ Dono sides ko $h$ se multiply karo (fraction clear karo): $h f (t_{n}, y_{n}) \approx y_{n + 1} - y_{n}$ Dono sides mein $y_{n}$ add karo (unknown $y_{n + 1}$ isolate karo): $y_{n + 1} = y_{n} + h jahan khade hain wahan slope f (t_{n}, y_{n})$

PICTURE. Neeche ka right triangle rise-over-run ko literal banata hai: run $= h$ bottom par, rise $= h \cdot f (t_{n}, y_{n})$ side par. Slope line follow karo aur $y_{n + 1}$ par land karo.

$h$ ::: step size — grid rungs $t_{n}$ aur $t_{n + 1}$ ke beech horizontal gap.
$y_{n}$ ::: time $t_{n} = t_{0} + nh$ par $y$ ka hamara current best estimate.
Discrete times kaise define hote hain? ::: $t_{n + 1} = t_{n} + h$ , toh $t_{n} = t_{0} + nh$ .

Step 3 — Euler's method: ek slope par trust karo, poora hop lo

KYA. Hum slope ek baar measure karte hain, interval ke start par, aur poore step $h$ ke liye us straight line par chalte hain.

KYUN yeh shakha hai. Sahi curve step ke dauran bend karta hai. Lekin Euler ne slope ko left edge par freeze kar diya, toh agar curve concave up hai toh Euler ki straight line sach se neeche jaati hai; concave down, toh upar jaati hai. Error us moment bak jaata hai jab hum ek slope par commit karte hain.

PICTURE. Dekho kaise red Euler segment sahi curve ko tangentially chodta hai aur phir drift karta hai jaise curve uske neeche upar curve karti hai. $t_{n + 1}$ par vertical gap is ek step ka local error hai.

Step 4 — RK4 idea: guess mat karo, aage jhaanko

KYA. Hum ek averaged hop commit karne se pehle chaar slope samples $k_{1}, k_{2}, k_{3}, k_{4}$ collect karenge.

KYUN chaar, aur kyun middle do baar. Sahi step ko $h$ mein Taylor series ke roop mein expand karo: exact answer hai $y_{n} + h y^{'} + \frac{h ^{2}}{2} y^{''} + \frac{h ^{3}}{6} y^{'''} + \frac{h ^{4}}{24} y^{''''} + \dots$ . Har slope sample $k_{i}$ , khud Taylor-expanded hoke, apne $h, h^{2}, h^{3}, h^{4}$ pieces contribute karta hai. Weights $1, 2, 2, 1$ over $6$ woh unique mixture hai jo blend ke $h^{1}, h^{2}, h^{3}, h^{4}$ pieces ko sahi series se term-by-term match karta hai; pehla mismatch $h^{5}$ par hai. Do midpoint samples $h^{2}$ aur $h^{3}$ corrections ka bulk uthate hain — wahi instinct jaise Simpson's Rule mein hai, jisme midpoint par bhi zyada weight hota hai.

PICTURE. Left panel: interval par chaar slope arrows (start $k_{1}$ , midpoint $k_{2}, k_{3}$ , end $k_{4}$ ). Right panel: ek skeletal "ledger" jo $h$ ki har power dikhata hai aur weighted blend kaise us power par true-series term cancel karta hai, first leftover $h^{5}$ par chodke.

$h$ ki kaunsi power pehli hai jo RK4 match nahi karta? ::: $h^{5}$ — $h^{4}$ tak sab exactly cancel ho jaata hai.

Step 5 — $k_{1}$ aur $k_{2}$ banana (start, phir pehli midpoint jhaanki)

KYA. $k_{1}$ hamari pehli honest reading hai. $k_{2}$ sirf midpoint scout karne ke liye $k_{1}$ use karta hai, phir us scouted spot par fresh slope measure karta hai.

KYUN. Midpoint slope curve ke bending ke liye correct karta hai jo akele $k_{1}$ miss kar gaya tha. Padhne se pehle midpoint par jaana zaroori hai (dono $t$ aur $y$ mein), warna hum sirf $k_{1}$ dobara paate.

PICTURE. Dashed blue trial-line $k_{1}$ se driven half-step hai; uski tip par solid arrow $k_{2}$ hai, visibly steeper kyunki curve chadh rahi hai.

$k_{1}$ ::: left edge par slope, same slope jo Euler use karta hai.
$k_{2}$ ::: $k_{1}$ ke saath trial half-step se pahunche midpoint par slope.

Step 6 — $k_{3}$ se refine karo aur $k_{4}$ se end tak pahuncho

KYA. $k_{3}$ ek doosra, improved midpoint estimate hai; $k_{4}$ $k_{3}$ par bana end-of-interval estimate hai.

KYUN chain $k_{1} \to k_{2} \to k_{3} \to k_{4}$ . Har stage pichle par khada hota hai taaki errors average hone par partially cancel ho jaayein. Yahi sequential dependence poora trick hai.

PICTURE. Do midpoint arrows ( $k_{2}$ pale-yellow, $k_{3}$ pink) almost ek doosre ke upar baithe hain lekin pink $k_{3}$ sach ke thoda aur nazdeek nudged hai; far-right arrow $k_{4}$ hai.

Step 7 — Weighted average: chaar jhankiyon se ek hop

KYA. Hum chaar readings ko ek averaged slope mein collapse karte hain aur uske saath ek single Euler-style hop lete hain.

KYUN $1, 2, 2, 1$ . Midpoint interval ke average behaviour ko best represent karta hai, toh ise double weight milti hai — exactly waise jaise Simpson's Rule apne middle sample ko load karta hai. Yahi exact numbers woh ek hi choice hain jo update ko true Taylor series se $h^{4}$ term tak agree karaate hain (Step 4).

PICTURE. Green averaged slope (chaar arrows se blended) ek clean hop drive karta hai jo essentially sahi curve par land karta hai — Step 3 ke mote red gap se is chhote green gap ko compare karo.

Step 8 — Edge & degenerate cases (taaki kuch surprise na kare)

KYA. Hum boundaries check karte hain: ek $f$ jo $t$ ignore kare, ek $f$ jo $y$ ignore kare, ek almost-flat step, aur ek backward step.

Autonomous $f (y)$ only ( $t$ nahi): $t + \frac{h}{2}$ aur $t + h$ shifts abhi bhi exist karte hain lekin kuch nahi karte, kyunki $f$ $t$ ignore karta hai. RK4 phir bhi kaam karta hai; shifts simply harmless hain — neeche ka figure dikhata hai do midpoint samples same slope value par land karte hain kyunki sirf $y$ matter karta hai.
Pure-time $f (t)$ only ( $y$ nahi, jaise $y^{'} = sin t$ ): ab stages ke andar har $y + \dots$ ignore ho jaata hai, aur RK4 exactly Simpson's rule mein $\int f d t$ ke liye collapse ho jaata hai. Yahan do midpoint samples $k_{2}, k_{3}$ identical ho jaate hain — woh combined weight $2 + 2 = 4$ ke saath Simpson ke ek midpoint mein fuse ho jaate hain. Yeh RK4–Simpson kinship ka sabse clear proof hai.
Flat step, $f \equiv 0$ : saare chaar $k$ 's $0$ hain, toh $y_{n + 1} = y_{n}$ . Dono methods sahi se move karne se mana karte hain. Koi division blow up nahi hota kyunki hum constant $6$ se divide karte hain, $h$ se kabhi nahi.
Backward step, $h < 0$ : har formula unchanged hai; $t$ aur $y$ simply left move karte hain. Integration reverse karne ke liye useful.

KYUN yeh dikhao. Agar tum sirf bending curve test karte ho toh tum bhool sakte ho ki $t$ -shifts ko bhoolna autonomous $f$ ke liye silently kaam karta hai aur non-autonomous $f$ ke liye toot jaata hai.

PICTURE (autonomous). Left: $f$ sirf $y$ par depend karta hai, toh dono midpoint times par read slope wahi same height-driven value hai — $t$ -shift kuch nahi badalta.

PICTURE (pure-time → Simpson). Yahan $f = sin t$ : chaar slope readings $t_{n}, t_{n} + \frac{h}{2}, t_{n} + h$ par land karti hain; do midpoint reads coincide karte hain aur $sin t$ ke neeche shaded area wahi hai jo RK4 compute karta hai — literally Simpson's rule.

Ek-picture summary

KYA. Ek figure, dono methods, same start, same $h$ : Euler ka red straight hop drift karta hua, RK4 ki chaar jhankiyaan ek green hop mein blend hoti hui jo sahi curve se chipki rehti hai. $t_{n + 1}$ par do error gaps dekho.

Recall Feynman: poora walkthrough simple words mein

Tum fog mein ek valley cross kar rahe ho. Tum aage equally spaced checkpoints mark karte ho — yeh hai grid $t_{n + 1} = t_{n} + h$ . Tumhare paas sirf ek compass hai jo, jahan bhi tum khade ho, batata hai zameen kis taraf tilts karti hai — yeh hai $f (t, y)$ . Euler apne pair ke neeche tilt ek baar read karta hai aur ek poora step us direction mein march karta hai. Lekin valley curve karti hai, toh tum trail se drift ho jaate ho; chhote steps lo aur kam drift hoga, lekin tum hamesha chalte rahoge. RK4 patient hai: yeh apne pair ke neeche tilt read karta hai ( $k_{1}$ ), aadhe step aage scout karta hai aur wahan read karta hai ( $k_{2}$ ), us behtar reading se dobara midpoint scout karta hai ( $k_{3}$ ), phir far end scout karta hai ( $k_{4}$ ). Yeh chaar readings blend karta hai — do middle wale double count karte hain, chhe se divide karte hain — ek trustworthy direction mein, aur us direction par ek hi step leta hai. Har step par chaar careful jhankiyaan tumhe almost perfectly trail par rakhti hain, isliye serious computing RK4 par lean karta hai. (Peeche chal rahe ho? Bas $h$ negative karo — har formula phir bhi hold karta hai.)

Connections

Taylor Series Expansion — $1, 2, 2, 1$ weights forced hain is series ko $h^{4}$ tak match karne se.
Simpson's Rule — Step 8 ka pure-time case dikhata hai RK4 literally Simpson ban jaata hai.
Finite Difference Approximation of Derivatives — Step 2 ka dropped limit forward difference hai.
Numerical Stability and Stiff ODEs — kya hota hai jab RK4 ka chhota step bhi safe nahi hota.
Adaptive Step Size (RK45 / Dormand–Prince) — solver ko $h$ khud choose karne do sensed error se.
scipy.integrate.solve_ivp — production tool jo yeh sab wrap karta hai.