Worked examples — Implementing ODE solvers from scratch — Euler, RK4

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

Tumne Euler aur RK4 parent note mein already dekha hai. Is page ka ek hi kaam hai: in dono update rules ko har tarah ki situation se guzaarna jo ek step-marching solver face kar sakta hai — growth, decay, oscillation, ek curve jo "galat" direction mein bend karti hai, itna bada step ki method blow up ho jaaye, ek real-world cooling problem, aur ek exam trap. Agar tum in sab ko follow kar sako, toh koi bhi ODE-solver question tumhe surprise nahi kar sakta.

Shuru karne se pehle, do tools ko plain words mein re-anchor karte hain taaki neeche koi bhi symbol aisa na ho jo tumhe diya na gaya ho.

"Slope" har example ka hero hai. Yeh picture apne dimaag mein rakho: ek hiker fog mein apne paon ke paas compass padh raha hai (Euler), versus ek hiker jo commit karne se pehle teen baar aage jhankta hai (RK4).

---

The scenario matrix

Har case class jo ek ODE-stepper se puchhi ja sakti hai, aur woh example jo use hit karta hai:

Cell	Kya ise distinct banata hai	Example
A. Growth $y^{\prime }=+ky$	slope badhti hai, curve upar bend karti hai → Euler undershoot karta hai	Ex 1
B. Decay $y^{\prime }=-ky$	slope shrink hoti hai, curve zero ki taraf bend karti hai → error ka sign flip hota hai	Ex 2
C. Non-autonomous $y^{\prime }=f(t)$	slope $t$ par depend karti hai — $k_2,k_3,k_4$ ke andar $t$-shifts matter karte hain	Ex 3
D. Oscillation $y^{\prime }=\cos t$	curve dono taraf turn karti hai; slope zero se guzarti hai	Ex 4
E. Degenerate / zero slope $y^{\prime }=0$	flat line; dono methods ko constant return karna chahiye	Ex 5
F. Big-step instability large $h$ on decay	ek critical $h$ ke baad dono methods blow up ho jaate hain	Ex 6
G. Real-world word problem	Newton cooling, $f$ set up karna aur units padhna zaroori hai	Ex 7
H. Exam twist	reuse-$k_1$ trap + error-scaling reasoning	Ex 8

Hum neeche aathon ko work out karte hain. Har ek mein pehle Forecast karo — steps padhne se pehle actually guess karo.

---

Example 1 — Cell A: exponential growth $y^{\prime }=y$, $y(0)=1$

Forecast: curve upar bend karti hai (slope badhti rehti hai). Euler slope sirf low-left start par padhta hai, isliye woh bahut-chhoti slope use karta hai → undershoot karna chahiye. Kaun sa method closer landega? Abhi guess karo.

1. Euler. $y_1=y_0+hf(t_0,y_0)=1+0.5\cdot (1)=1.5$. Yeh step kyun? Euler assume karta hai ki starting slope $f(0,1)=1$ poore half-step tak hold karti hai. Kyunki real slope step ke dauran badhti hai, $1.5$ ka value $1.6487$ se kam padh jaata hai.
2. RK4 $k_1$. $k_1=f(0,1)=1$. Kyun? Bilkul wahan slope jahan hum khade hain — wohi reading jo Euler ne use ki.
3. RK4 $k_2$. $k_2=f(0.25,1+0.25\cdot 1)=1.25$. Kyun? Midpoint $t=0.25$ par peek karo, $y$ ko $k_1$ se aage nudge karke. Wahan slope badi hai — yeh undershoot ko already correct kar raha hai.
4. RK4 $k_3$. $k_3=f(0.25,1+0.25\cdot 1.25)=1.3125$. Kyun? Midpoint dobara karo lekin $y$ ko better slope $k_2$ se nudge karo. Har stage pichle ko refine karta hai.
5. RK4 $k_4$. $k_4=f(0.5,1+0.5\cdot 1.3125)=1.65625$. Kyun? End-of-step slope, $k_3$ use karke wahan pahuncha.
6. RK4 combine. $y_1=1+\frac{0.5}{6}(1+2(1.25)+2(1.3125)+1.65625)=1+\frac{0.5}{6}(7.71875)\approx 1.64323$. $1,2,2,1$ kyun? Weighted mean slope; midpoints double count karti hain kyunki woh step ke average ko best represent karti hain. Inner sum note karo: $1+2.5+2.625+1.65625=7.78125$… ruko — dhyan se count karo: $1+2(1.25)+2(1.3125)+1.65625=1+2.5+2.625+1.65625=7.78125$.

Arithmetic clearly bataate hain taaki kuch bhi fudge na ho: weighted sum hai $k_1+2k_2+2k_3+k_4=1+2(1.25)+2(1.3125)+1.65625=7.78125,$ isliye $y_1=1+\frac{0.5}{6}(7.78125)=1+0.6484375=\mathrm{1.6484375.}$

Verify: true $e^{0.5}=1.64872$. Euler error $=0.1487$; RK4 error $=0.00028$ — RK4 ek single step mein ~530× closer hai. Euler ke liye undershoot forecast confirmed. ✓

Figure padho: Figure s01 mein white curve true $e^t$ hai. Amber line follow karo — yeh Euler ka straight hop hai, aur yeh white curve ke neeche end hota hai: woh gap woh undershoot hai jo tumne forecast ki thi. Cyan square RK4 ka endpoint hai, practically white curve par baith raha hai. Visual takeaway: ek straight Euler segment ek bending curve ko hug nahi kar sakta, lekin RK4 ka four-slope average kar sakta hai.

---

Example 2 — Cell B: exponential decay $y^{\prime }=-y$, $y(0)=1$

Forecast: ab curve zero ki taraf upar se bend karti hai — yeh convex hai, girte girte flat hoti jaati hai. Euler steep starting slope $-1$ ko poore step ke liye use karta hai. Kya yeh true value ke upar ya neeche land karega?

1. Euler. $y_1=1+0.5\cdot (-1)=0.5$. Yeh step kyun? Starting slope $f(0,1)=-1$ woh sabse steep descent hai jo curve yahan kabhi karti hai; use poore step ke liye use karne se hum bahut zyada neeche chale jaate hain.
2. Compare. True $e^{-0.5}=0.6065$, isliye Euler $0.5$ par undershoot karta hai (neeche land karta hai). Yeh note kyun karte hain? Cell A (growth) mein Euler curve ke neeche land kiya tha, lekin wahan "neeche" matlab tha "bahut kam rise". Yahan "neeche" matlab hai "bahut tezi se gira". Lesson: Euler hamesha starting slope ki taraf error karta hai, aur woh high hai ya low yeh curvature par depend karta hai — lekin growth aur is decay dono ke liye yeh true curve ke neeche land karta hai kyunki dono convex hain.

Verify: $|0.6065-0.5|=0.1065$ actual error. Leading-order Taylor estimate hai $\frac{h^2}{2}|y^{\prime \prime }|=\frac{0.25}{2}\cdot 0.6065\approx 0.076$ — sahi order of magnitude; expansion mein higher-order terms true $0.107$ tak ke gap ka hisaab lagate hain. ✓

---

Example 3 — Cell C: non-autonomous $y^{\prime }=2t$, $y(0)=0$

Forecast: yahan $f$ bilkul $y$ ko ignore karta hai — slope purely time se set hoti hai. Yeh exactly woh case hai jahan $k_2,k_3,k_4$ ke andar $t$-shifts bhoolna tumhe silently barbad kar deta hai. Guess karo: kya RK4 ek parabola ke liye exact ya approximate hoga?

Is poore example mein step size $h=1$ hai, isliye jahan bhi $\frac{h}{2}$ dikhe iska matlab $0.5$ hai aur $\frac{h}{6}$ ka matlab $\frac{1}{6}$ hai — main formulas mein $h$ rakhta hoon aur end mein substitute karta hoon taaki general shape visible rahe.

1. $k_1$. $k_1=f(0,0)=2\cdot 0=0$. Kyun? Bilkul start par slope, $t=0$, zero hai.
2. $k_2$. $k_2=f(0+\frac{h}{2},0+\frac{h}{2}\cdot 0)=f(0.5,0)=2\cdot 0.5=1$. Kyun? Hume zaroor $t=0+\frac{h}{2}=0.5$ par evaluate karna hai, $t=0$ par nahi. $y$-nudge irrelevant hai kyunki $f$ $y$ nahi padhta, lekin $t$-shift sab kuch hai.
3. $k_3$. $k_3=f(0.5,0+\frac{h}{2}\cdot 1)=f(0.5,0.5)=2\cdot 0.5=1$. Kyun? Same $t=0.5$ → same slope $1$. Confirm karta hai ki jab $f$ $y$-independent ho toh do midpoint slopes agree karti hain.
4. $k_4$. $k_4=f(0+h,\dots )=f(1,\dots )=2\cdot 1=2$. Kyun? End-of-step slope $t=0+h=1$ par.
5. Combine. $h$ factor ko explicit rakhte hue: $y_1=0+\frac{h}{6}(k_1+2k_2+2k_3+k_4)=\frac{1}{6}(0+2(1)+2(1)+2)=\frac{1}{6}(6)=1.$ Kyun? Chaar slopes ka weighted average, $\frac{h}{6}=\frac{1}{6}$ se scale kiya kyunki $h=1$, exactly land karta hai.

Verify: true $t^2$ at $t=1$ is $1$ → RK4 error $=0$. Exact kyun? RK4 Taylor series ko $h^4$ tak match karta hai; ek quadratic mein $h^2$ ke baad koi terms nahi hain, isliye RK4 uske liye exact hai. ✓

---

Example 4 — Cell D: oscillation $y^{\prime }=\cos t$, $y(0)=0$

Forecast: slope $\cos 0=1$ se shuru hoti hai (tezi se badhti) lekin is quarter-period ke end tak $0$ tak ghatti hai. Euler interval ki sabse badi slope ko lock in kar leta hai. Overshoot ya undershoot?

1. Euler. $y_1=0+\frac{\pi }{2}\cdot \cos 0=\frac{\pi }{2}\cdot 1\approx 1.5708$. Yeh step kyun? Euler poore quarter-turn mein start slope $1$ par trust karta hai. Lekin real slope ghatti rehti hai zero tak, isliye constant slope $1$ badly overshoot karta hai.
2. Interpret. True value $1$, Euler deta hai $1.5708$ → overshoot of $0.5708$. Ex 1 se flip kyun? Growth mein missed slope baad mein badi thi (undershoot). Yahan missed slope baad mein chhhoti hai, isliye same "trust-the-start" rule ab overshoot karta hai. Isliye curvature sign, method nahi, error ki direction decide karta hai.

Verify: $\frac{\pi }{2}-1=0.5708$. Overshoot confirmed. ✓

Figure padho: Figure s02 mein white curve $\sin t$ hai jo $1$ tak badhti hai phir flat hoti hai. Amber straight line Euler ka step hai — yeh constant steep slope $1$ par chadhti hai aur white curve ke upar shoot karti hai, $\pi /2\approx 1.57$ par end hoti hai jabki true endpoint (cyan square) $1$ par baitha hai. s01 se contrast karo: same "trust the start slope" rule, lekin yahan curve flat hoti hai isliye error opposite direction mein jaata hai. Woh sign-flip hi poora lesson hai.

---

Example 5 — Cell E: degenerate zero slope $y^{\prime }=0$, $y(0)=7$

Forecast: slope har jagah zero hai. Ek correct solver ko $7$ unchanged return karna chahiye. Agar tumhara drift kare, toh bug hai.

1. Euler. $y_1=7+10\cdot 0=7$. Kyun? Zero slope × koi bhi step = koi movement nahi.
2. RK4. $k_1=k_2=k_3=k_4=0$, isliye $y_1=7+\frac{10}{6}(0)=7$. Kyun? Har peek zero padhta hai; zeros ka weighted average zero hai.

Verify: dono kisi bhi $h$ ke liye exactly $7$ return karte hain. ✓ Yeh sabse sasta lekin sabse important unit test hai — koi bhi solver jo tum likhte ho usp pehle run karo.

---

Example 6 — Cell F: big-step instability $y^{\prime }=-3y$, $y(0)=1$

Pehle, ek symbol ki zaroorat hai. Jab ek ODE ki form $y^{\prime }=\lambda y$ ho, toh number $\lambda$ (Greek "lambda") sirf growth/decay coefficient hai — wahi role jo matrix ki rows mein $\pm k$ play karta hai "$y^{\prime }=+ky$" aur "$y^{\prime }=-ky$". Hamare problem $y^{\prime }=-3y$ ke liye hum simply $\lambda =-3$ padh lete hain (ek decay, kyunki negative hai). Neeche sab kuch $\lambda$ ke saath likha hai taaki rule reusable rahe, lekin hamare liye $\lambda =-3$.

Forecast: true solution monotonically zero ki taraf decay karti hai. Lekin dekho bada step kya karta hai. Guess karo: kya Euler ke numbers chhhote aur positive rahenge?

1. Step 1. $y_1=1+1\cdot (-3\cdot 1)=1-3=-2$. Disaster kyun? Euler step ko slope $-3$ se multiply karta hai. $h=1$ ke saath correction $-3$ zero se past negative territory mein overshoot kar jaata hai — ek decay ke liye physically impossible.
2. Step 2. $y_2=-2+1\cdot (-3\cdot -2)=-2+6=4$. Kyun? Overshoot sign flip karta hai, slope bhi flip hoti hai, aur error har step badhta hai.
3. Step 3. $y_3=4+1\cdot (-3\cdot 4)=4-12=-8$. Kyun? Har step magnitude $|1+h\lambda |=|1+1\cdot (-3)|=|-2|=2$ se multiply hoti hai. Yeh hamesha double hota rehta hai — unstable.

Euler ki stability rule: forward Euler $y^{\prime }=\lambda y$ par bounded tabhi rehta hai jab $|1+h\lambda |\le 1$. Yahan $\lambda =-3$, isliye hume chahiye $|1-3h|\le 1\Rightarrow 0\le h\le \frac{2}{3}$. Hamara $h=1$ ise violate karta hai.

Kya RK4 bach jaata hai? Nahi — higher order ka matlab unconditionally stable nahi hota. $y^{\prime }=\lambda y$ ke liye har RK4 step $y$ ko polynomial se multiply karta hai $R(z)=1+z+\frac{z^2}{2}+\frac{z^3}{6}+\frac{z^4}{24},z=h\lambda ,$ ($e^z$ ki Taylor series $z^4$ par truncated — exactly chaar stages). RK4 bounded tabhi hai jab $|R(z)|\le 1$. $h=1,\lambda =-3$ ke saath $z=-3$ milta hai aur $R(-3)=1-3+4.5-4.5+3.375=1.375,$ jiska magnitude $>1$ hai, isliye RK4 bhi $h=1$ ke saath is problem par blow up hoga. RK4 ki stability limit yahan roughly $h\lesssim 0.93$ hai (jahan $|R(h\lambda )|=1$) — Euler ke $\frac{2}{3}$ se bada, lekin phir bhi finite. Lesson: stiff decay ke liye dono methods ke liye $h$ chhota karna padega; Numerical Stability and Stiff ODEs dekho.

Verify: safe $h=0.5$ ke saath: Euler $y_1=1+0.5(-3)=-0.5$ — phir bhi negative lekin $|1-1.5|=0.5\le 1$, isliye magnitudes ab shrink karti hain, badhti nahi. Sequence $-0.5,0.25,-0.125,\cdots \to 0$. Bounded. RK4 ke liye $h=0.5$, $z=-1.5$ deta hai $R(-1.5)=1-1.5+1.125-0.5625+0.2109=0.2734$, $|R|<1$, isliye RK4 bhi bounded hai. ✓

Figure padho: Figure s03 mein white curve true decay $e^{-3t}$ hai jo zero se chipki hui hai. Amber zig-zag Euler hai $h=1$ ke saath: yeh white dashed zero-line ke upar aur neeche fling hota hai, har swing pichle se bada — woh doubling jo tumne compute ki. Cyan markers Euler hain safe $h=0.5$ ke saath: same alternating signs lekin har swing chhhoti, zero mein spiral hoti. Picture "unstable vs bounded" ko zig-zag ke badhne ya ghattne ka matter bana deta hai.

---

Example 7 — Cell G: real-world word problem (Newton's cooling)

Forecast: cup pehle tezi se cool hoti hai (bada gap → steep slope), phir slowly. Poore 10-minute step mein Euler sirf initial steep slope use karega. Kya woh true temperature ko overshoot ya undershoot karega?

Setup: $f(t,T)=-0.1(T-20)$. Units: $[{}^{\circ }\text{C}/\text{min}]$. Start $T_0=90$.

1. Euler. slope $=-0.1(90-20)=-7{}^{\circ }\text{C/min}$. $T_1=90+10\cdot (-7)=90-70=20$. Suspicious kyun? Euler predict karta hai ki coffee room temperature par exactly pahunch jaati hai, initial fast cooling ko poore 10 minutes ke liye use karke. Real cooling slow ho jaati hai, isliye itna cool nahi ho sakta.
2. RK4 $k_1$. $k_1=-0.1(90-20)=-7$.
3. RK4 $k_2$. $T$-nudge: $90+5\cdot (-7)=55$; $k_2=-0.1(55-20)=-3.5$. Kyun? midpoint slope gentler hai kyunki cup already cool ho gayi.
4. RK4 $k_3$. $90+5\cdot (-3.5)=72.5$; $k_3=-0.1(72.5-20)=-5.25$. Kyun? midpoint ko $k_2$ se refine karo.
5. RK4 $k_4$. $90+10\cdot (-5.25)=37.5$; $k_4=-0.1(37.5-20)=-1.75$. Kyun? end slope, sab mein sabse gentle.
6. Combine. Slopes ka weighted sum hai $k_1+2k_2+2k_3+k_4=-7+2(-3.5)+2(-5.25)+(-1.75)=-7-7-10.5-1.75=-26.25$. Isliye $T_1=90+\frac{10}{6}(-26.25)=90-43.75=46.25{}^{\circ }\text{C}$.

Verify: exact solution $T(t)=20+70e^{-0.1t}$, isliye $T(10)=20+70e^{-1}\approx 20+25.7517=45.7517{}^{\circ }\text{C}$. Euler ka $20^{\circ }$ ek gross overshoot in cooling hai ($25.75^{\circ }$ off); RK4 ka $46.25^{\circ }$ sirf $0.498^{\circ }$ off hai ek giant 10-minute step ke bawajood. ✓ Units sab ${}^{\circ }\text{C}$ mein. Forecast (Euler bahut zyada cool karta hai) confirmed.

---

Example 8 — Cell H: exam twist (the reuse-$k_1$ trap + scaling)

Forecast (a): $k_1$ reuse karne se $k_2$ aur $k_3$ identical ho jaate hain — method ek independent slope kho deta hai. Kya $y_1$ correct $2.7083$ se bada hoga ya chhota?

Part (a):

1. Correct stages (parent se): $k_1=1,k_2=1.5,k_3=1.75,k_4=2.75$, correct $y_1=1+\frac{1}{6}(1+3+3.5+2.75)=1+\frac{10.25}{6}\approx 2.70833$ deta hai.
2. Buggy $k_3$. $k_3=f(0.5,1+0.5\cdot k_1)=f(0.5,1.5)=1.5$ ($k_2$ ke same). Damage kyun? $k_3$ ab $k_2$ ko refine nahi karta; do midpoint peeks ab identical information carry karte hain, isliye averaging kamzor hai.
3. Buggy $k_4$. $k_4=f(1,1+1\cdot k_3)=f(1,2.5)=2.5$. Kyun changed? $k_4$ galat $k_3$ par build karta hai.
4. Combine. $y_1=1+\frac{1}{6}(1+2(1.5)+2(1.5)+2.5)=1+\frac{1}{6}(9.5)\approx 2.58333$. Kyun worse? True $e=2.71828$; correct RK4 ne $2.70833$ diya (error $0.00995$), buggy deta hai $2.58333$ (error $0.135$) — bug ~$13.5\times$ worse hai, Euler-jaisi quality ki taraf collapse karta hai.

Part (b):

1. Scaling law. RK4 error $\propto h^4$. $h=0.1$ se $h=0.025$ jaana factor $\frac{0.1}{0.025}=4$ reduction hai step mein. $h^4$ kyun? Global error order $O(h^4)$ hai (parent), isliye error $4^4=256$ se shrink hota hai.
2. Predict. $2.1\times 10^{-6}/256\approx 8.2\times 10^{-9}$.

Verify: (a) buggy inner sum $\frac{9.5}{6}=1.58333$, isliye $y_1=2.58333$; correct RK4 value $2.70833$ se iska gap $0.125$ hai. ✓ (b) $2.1\times 10^{-6}/256=8.20\times 10^{-9}$. ✓

Recall Kaun se cell ne kaun sa lesson sikhaya?

Growth vs decay :::: curvature sign, method nahi, over/undershoot direction decide karta hai (Ex 1, 2, 4). Non-autonomous :::: $k_2,k_3,k_4$ ke andar $t$-shift kabhi skip mat karo (Ex 3). Zero slope :::: dono methods ko kisi bhi $h$ ke liye constant return karna chahiye — tumhara pehla unit test (Ex 5). Big step :::: Euler bounded tabhi hai jab $|1+h\lambda |\le 1$; RK4 bounded tabhi hai jab $|R(h\lambda )|\le 1$ — dono ki finite limit hoti hai (Ex 6). Reuse-$k_1$ :::: har stage ko next ko feed karna chahiye warna tum first order ki taraf wapas aa jaate ho (Ex 8a).

---

Connections

• Taylor Series Expansion — RK4 Ex 3 mein quadratic par exact kyun hai, aur Ex 2 mein $\frac{h^2}{2}{y}^{\prime \prime }$ error ka source.
• Simpson's Rule — $1,2,2,1$ averaging philosophy jo Ex 1 ke RK4 step ke peeche hai.
• Finite Difference Approximation of Derivatives — woh forward difference jo hai Euler.
• Numerical Stability and Stiff ODEs — Ex 6 mein blow-up puri tarah explain kiya gaya hai.
• Adaptive Step Size (RK45 / Dormand–Prince) — solvers $h$ kaise choose karte hain taaki Ex 6 kabhi na ho.
• scipy.integrate.solve_ivp — woh production tool jo yeh sab wrap karta hai.