Worked examples — Systems of ODEs — RK4 for systems

4.8.27 · D3 · Maths › Numerical Methods › Systems of ODEs — RK4 for systems

Kisi bhi example se pehle, ek reminder us machinery ka jo hum baar baar reuse karte hain.

Yahan "vector" ka matlab sirf numbers ki ek ordered list hai, ek number har unknown ke liye. Agar hum do quantities $y_1,y_2$ track kar rahe hain, to $\mathbf{y}=(y_1,y_2)$ ek plane mein ek point hai, aur $\mathbf{f}$ ek aisa pair return karta hai — direction aur speed jisme har quantity abhi move kar rahi hai. Notation ka yahi matlab hai.

---

Scenario matrix

Is topic ka har problem in cells mein se kisi ek mein aata hai. Neeche ke worked examples mein se har ek ko us cell ke saath tag kiya gaya hai jo wo cover karta hai, taaki milke poora grid fill ho jaye.

Cell	Kya alag banata hai	Example
A. Coupled, mixed signs	$f_1,f_2$ dono doosre variable pe depend karte hain, opposite signs (rotation)	Ex 1
B. Decoupled system	equations ek doosre ko reference nahi karte — ek sanity check	Ex 2
C. Reduced 2nd-order	ek single higher-order ODE ko system mein convert kiya	Ex 3
D. Nonlinear word problem	unknowns ke products, real units (predator–prey)	Ex 4
E. Degenerate / steady state	initial state jahan $\mathbf{f}=0$ (kuch nahi hona chahiye)	Ex 5
F. All-positive growth	dono variables grow karte hain, positive slopes throughout	Ex 6
G. Limiting / stiff twist	ek component itni tez decay karta hai ki step size matter karta hai	Ex 7
H. Non-autonomous	$\mathbf{f}$ explicitly $t$ pe depend karta hai, sirf $\mathbf{y}$ pe nahi	Ex 8

"Sign behaviour" ke columns (dono slopes negative, dono positive, mixed, zero) aur "structure" ke rows (coupled / decoupled / reduced / nonlinear / stiff / time-dependent) sab cover ho jaate hain. Neeche ki figure in chaar cells mein se chaar ko state space mein motions ke roop mein plot karti hai, taaki aap arithmetic grind karne se pehle dekh sako ki scenarios kitne alag hain.

Figure ko ek map ki tarah padhiye: blue arc (Cell A, Ex 1) ek point hai jo circle pe curve kar raha hai — mixed signs isse rotate karte hain; orange arrow (Cell B, Ex 2) dikhata hai ki ek coordinate shrink ho rahi hai jabki doosri grow kar rahi hai — ek decoupled system; green curve (Cell F, Ex 6) mein dono coordinates milke baahir race kar rahe hain — pure growth; red dot (Cell E, Ex 5) ek fixed point pe frozen baitha hai — kuch nahi hil raha. Jab neeche har example in paths mein se ek trace kare to is picture ko apne zehan mein rakhna.

---

Ex 1 — Cell A: coupled rotation, mixed signs

Forecast: ye ek point ki equations hain jo circle pe spin kar raha hai (position aur velocity ek doosre ko opposite signs ke saath feed kar rahe hain). Ek chhote step ke baad $y_1$ ko $1$ se thoda neeche jaana chahiye aur $y_2$ ko thoda negative hona chahiye. Kyunki RK4 order 4 hai, ~4 correct digits expect karo.

Yahan $\mathbf{f}(t,\mathbf{y})=(y_2,-y_1)$. Yeh figure mein blue arc hai.

1. ${\mathbf{k}}_1=(y_2,-y_1)=(0,-1)$. Yeh step kyun? Current state $(1,0)$ pe slope — yeh hai jahan $\mathbf{y}$ abhi ja raha hai.
2. Midpoint state $=\mathbf{y}+\frac{h}{2}{\mathbf{k}}_1=(1,0)+0.1(0,-1)=(1,-0.1)$, to ${\mathbf{k}}_2=(-0.1,-1)$. Yeh step kyun? Step ke predicted middle mein slope dobara measure karo — dono components advanced state use karte hain.
3. $\mathbf{y}+\frac{h}{2}{\mathbf{k}}_2=(1,0)+0.1(-0.1,-1)=(0.99,-0.1)$, to ${\mathbf{k}}_3=(-0.1,-0.99)$. Yeh step kyun? Ek refined midpoint slope, ${\mathbf{k}}_2$ ko correct karta hua.
4. $\mathbf{y}+h{\mathbf{k}}_3=(1,0)+0.2(-0.1,-0.99)=(0.98,-0.198)$, to ${\mathbf{k}}_4=(-0.198,-0.98)$. Yeh step kyun? Interval ke end ka slope.
5. Combine weights $1,2,2,1$ over $6$ ke saath: $y_{1,1}=1+\frac{0.2}{6}(0+2(-0.1)+2(-0.1)+(-0.198))=1+\frac{0.2}{6}(-0.598)\approx \mathrm{0.980067.}$ $y_{2,1}=0+\frac{0.2}{6}(-1+2(-1)+2(-0.99)+(-0.98))=\frac{0.2}{6}(-5.96)\approx -\mathrm{0.198667.}$ Yeh step kyun? Ab hum charon slope vectors ko ek best estimate mein blend karte hain step ki average velocity ka (beech ke slopes double — Simpson weighting) aur $\mathbf{y}$ ko us average times $h$ se move karte hain. Yahi state ka actual advance hai.

Verify: $\cos 0.2=0.980067$ (~5 decimals tak match), $-\sin 0.2=-0.198669$ (~4 decimals tak match). Ek step mein agreement. Forecast confirmed.

---

Ex 2 — Cell B: decoupled sanity check

Forecast: dono equations kabhi ek doosre ka zikr nahi karti. Ek decoupled system pe RK4 ko, component-by-component, exactly wahi reproduce karna chahiye jo scalar RK4 har equation pe alag se dega. $y_1$ shrink hogi (negative slope), $y_2$ tez grow karegi (slope $2y_2$). Yeh figure mein orange arrow hai.

Yahan $\mathbf{f}=(-y_1,2y_2)$.

1. ${\mathbf{k}}_1=(-1,2)$. Kyun? $(1,1)$ pe slope.
2. State $=(1,1)+0.25(-1,2)=(0.75,1.5)$, ${\mathbf{k}}_2=(-0.75,3)$. Kyun? Midpoint slope.
3. State $=(1,1)+0.25(-0.75,3)=(0.8125,1.75)$, ${\mathbf{k}}_3=(-0.8125,3.5)$. Kyun? Refined midpoint.
4. State $=(1,1)+0.5(-0.8125,3.5)=(0.59375,2.75)$, ${\mathbf{k}}_4=(-0.59375,5.5)$. Kyun? End slope.
5. Combine: $y_{1,1}=1+\frac{0.5}{6}(-1+2(-0.75)+2(-0.8125)+(-0.59375))=1+\frac{0.5}{6}(-4.71875)\approx \mathrm{0.606771.}$ $y_{2,1}=1+\frac{0.5}{6}(2+2(3)+2(3.5)+5.5)=1+\frac{0.5}{6}(20.5)\approx \mathrm{2.708333.}$ Yeh step kyun? Har component ko apna $1,2,2,1$ blend milta hai — kyunki equations decoupled hain, yeh literally scalar RK4 hai jo do baar parallel mein run ho raha hai.

Verify: $e^{-0.5}=0.606531$ (error $\approx 2.4\times 10^{-4}$), $e^1=2.718282$ (error $\approx 0.01$). Bada growth rate ($2y_2$) per unit $t$ ek bada step banata hai, isliye uska error bada hai — expected. Dono wahi same numbers hain jo aapko scalar RK4 har line pe akele run karne se milte.

---

Ex 3 — Cell C: ek 2nd-order ODE ko reduce karna

Forecast: yeh reduction trick hai. $y(0)=y^{\prime }(0)=1$ ke saath exact solution $e^t$ hai, to $u_1$ ko $e^{0.5}=1.6487$ track karna chahiye.

Set karo $u_1=y,u_2=y^{\prime }$. Tab $u_1^{\prime }=u_2,u_2^{\prime }=u_1,$ to $\mathbf{f}=(u_2,u_1)$, starting at $\mathbf{u}=(1,1)$.

1. ${\mathbf{k}}_1=(1,1)$. Kyun? $(1,1)$ pe slope: $u_1^{\prime }=u_2=1$, $u_2^{\prime }=u_1=1$.
2. State $=(1,1)+0.25(1,1)=(1.25,1.25)$, ${\mathbf{k}}_2=(1.25,1.25)$. Kyun? Midpoint slope — components swap karo.
3. State $=(1,1)+0.25(1.25,1.25)=(1.3125,1.3125)$, ${\mathbf{k}}_3=(1.3125,1.3125)$. Kyun? Refined midpoint.
4. State $=(1,1)+0.5(1.3125,1.3125)=(1.65625,1.65625)$, ${\mathbf{k}}_4=(1.65625,1.65625)$. Kyun? End slope.
5. Combine (dono components yahan identical hain kyunki $u_1=u_2$ throughout hai): $u_{1,1}=1+\frac{0.5}{6}(1+2(1.25)+2(1.3125)+1.65625)=1+\frac{0.5}{6}(7.78125)\approx \mathrm{1.648438.}$ Yeh step kyun? Charon slopes ko fixed $1,2,2,1$ weights ke saath blend karo $u_1$ (aur $u_2$) ko ek step advance karne ke liye; reduction ne is averaging mein kuch nahi badla.

Verify: $e^{0.5}=1.648721$; error $\approx 2.8\times 10^{-4}$. 2nd-order ODE usi systems recipe se solve hua — reduction karne ki koi cost nahi thi.

---

Ex 4 — Cell D: nonlinear predator–prey word problem

Forecast: $RF$ ke products isko nonlinear banate hain — coupling multiplicative hai. Shuru mein $R^{\prime }=2-2=0$ aur $F^{\prime }=2-1=1$: rabbits momentarily steady, foxes rise kar rahe hain. To expect karo $R\approx 2$ (nearly flat) aur $F$ roughly $1.1$ tak.

Yahan $\mathbf{f}=(R-RF,RF-F)$, state $(2,1)$.

1. ${\mathbf{k}}_1=(2-2\cdot 1,2\cdot 1-1)=(0,1)$. Kyun? Current slopes.
2. State $=(2,1)+0.05(0,1)=(2,1.05)$; ${\mathbf{k}}_2=(2-2\cdot 1.05,2\cdot 1.05-1.05)=(-0.1,1.05)$. Kyun? Midpoint slope, dono populations advance.
3. State $=(2,1)+0.05(-0.1,1.05)=(1.995,1.0525)$; $RF=1.995\cdot 1.0525=2.09974$, to ${\mathbf{k}}_3=(1.995-2.09974,2.09974-1.0525)=(-0.104738,1.047238)$. Kyun? Refined midpoint.
4. State $=(2,1)+0.1(-0.104738,1.047238)=(1.989526,1.104724)$; $RF=1.989526\cdot 1.104724=2.197931$, to ${\mathbf{k}}_4=(1.989526-2.197931,2.197931-1.104724)=(-0.208405,1.093207)$. Kyun? End slope.
5. Combine: $R_1=2+\frac{0.1}{6}(0+2(-0.1)+2(-0.104738)+(-0.208405))=2+\frac{0.1}{6}(-0.617881)\approx \mathrm{1.989702.}$ $F_1=1+\frac{0.1}{6}(1+2(1.05)+2(1.047238)+1.093207)=1+\frac{0.1}{6}(6.287683)\approx \mathrm{1.104795.}$ Yeh step kyun? Chahe $\mathbf{f}$ nonlinear ho, averaging step unchanged hai — hum phir bhi charon slope vectors ko $1,2,2,1$ over $6$ se weight karte hain. Nonlinearity sirf affect karti hai ki har ${\mathbf{k}}_i$ kaise compute hua, unhe blend karne ke tarike ko nahi.

Verify: rabbits $\approx 1.9897$ hazaar (barely girre, jaise forecast tha — wo apne peak ke near the), foxes $\approx 1.1048$ hazaar (rise kiye, jaise forecast tha). Units throughout animals ke hazaar hain; koi negative populations nahi aayi.

---

Ex 5 — Cell E: degenerate steady state

Forecast: $(1,1)$ pe $R^{\prime }=1-1=0$ aur $F^{\prime }=1-1=0$ milte hain — field exactly zero hai. Ek steady state ka matlab hai kuch nahi hona chahiye: answer dobara $(1,1)$ hona chahiye, machine precision tak. Yeh woh degenerate case hai jo sloppy code ko pakad leta hai — figure mein red dot.

1. ${\mathbf{k}}_1=(1-1,1-1)=(0,0)$. Kyun? Fixed point pe slope definition se zero hai.
2. State $=(1,1)+0.15(0,0)=(1,1)$; ${\mathbf{k}}_2=(0,0)$. Kyun? Hum kabhi $(1,1)$ se bahar nahi gaye, to slope abhi bhi zero hai.
3. State $=(1,1)$; ${\mathbf{k}}_3=(0,0)$. Kyun? Same reasoning.
4. State $=(1,1)$; ${\mathbf{k}}_4=(0,0)$. Kyun? Same.
5. Combine: ${\mathbf{y}}_1=(1,1)+\frac{0.3}{6}(0+0+0+0)=(1,1)$. Yeh step kyun? Chaar zero vectors ka $1,2,2,1$ blend zero hai, to average velocity zero hai aur state advance nahi hota — exactly wahi jo ek equilibrium demand karta hai.

Verify: state unchanged hai — RK4 equilibria ko exactly respect karta hai (kyunki har ${\mathbf{k}}_i=0$). Agar aapke program ne "akele ek component advance kiya hota," to shayad fixed point se thoda hat jaata; synchronized recipe aisa nahi karti.

---

Ex 6 — Cell F: all-positive growth

Forecast: $\mathbf{y}(0)=(1,1)$ aur field ke saath jo "har component doosre ki value pe grow kare," dono slopes throughout positive rehte hain — pure exponential growth. Expect karo ki dono $e^{0.4}=1.4918$ reach karen. Yeh figure mein green curve hai.

Yahan $\mathbf{f}=(y_2,y_1)$, state $(1,1)$.

1. ${\mathbf{k}}_1=(1,1)$. Kyun? Start pe slope.
2. State $=(1,1)+0.2(1,1)=(1.2,1.2)$, ${\mathbf{k}}_2=(1.2,1.2)$. Kyun? Midpoint.
3. State $=(1,1)+0.2(1.2,1.2)=(1.24,1.24)$, ${\mathbf{k}}_3=(1.24,1.24)$. Kyun? Refined midpoint.
4. State $=(1,1)+0.4(1.24,1.24)=(1.496,1.496)$, ${\mathbf{k}}_4=(1.496,1.496)$. Kyun? End slope.
5. Combine (dono components identical): $y_{1,1}=1+\frac{0.4}{6}(1+2(1.2)+2(1.24)+1.496)=1+\frac{0.4}{6}(7.376)\approx \mathrm{1.491733.}$ Yeh step kyun? Chaar positive slopes ko $1,2,2,1$ over $6$ se weight karo aur aage badhao; doubled midpoints exponential ki accelerating curve ko accurate rakhte hain.

Verify: $e^{0.4}=1.491825$; error $\approx 9\times 10^{-5}$. Dono slopes forecast ke mutabiq positive rahi.

---

Ex 7 — Cell G: limiting / stiff twist

Forecast: true solution $(e^{-20t},e^{-t})$ hai; $t=0.2$ pe woh $(e^{-4},e^{-0.2})=(0.0183,0.8187)$ hai. Lekin fast component rate $20$ ke saath decay karta hai, aur $h\cdot 20=4$ bada hai. $y^{\prime }=\lambda y$ ke liye RK4 ki stability limit roughly $|h\lambda |\lesssim 2.78$ chahiye. $|h\lambda |=4$ ke saath hum iske bahar hain — expect karo ki $y_1$ estimate reason se aage blow karega (possibly wrong sign ya magnitude), jabki $y_2$ (rate $1$, $|h\lambda |=0.2$) theek rahega.

$y_1$ compute karo (ek scalar RK4 $f=-20y_1$ ke saath, $h=0.2$, to $h\lambda =-4$). Linear $y^{\prime }=\lambda y$ ke liye, ek RK4 step multiply karta hai $R(z)=1+z+\frac{z^2}{2}+\frac{z^3}{6}+\frac{z^4}{24},z=h\lambda =-4.$

1. $z=-4$. Kyun? Step ka stability argument.
2. $R(-4)=1-4+8-\frac{64}{6}+\frac{256}{24}=1-4+8-10.6667+10.6667=5.0$. Kyun? Polynomial mein plug karo.
3. To $y_{1,1}=1\cdot R(-4)=5.0$. Kyun? Ek step state ko $R(z)$ se scale karta hai — yeh factor hai ek linear field ke liye four-slope $1,2,2,1$ blend ka combined effect.
4. $y_2$ ke liye: $z=-0.2$, $R(-0.2)=1-0.2+0.02-0.001333+0.0000667\approx 0.818733$. Kyun? Chhota $|z|$, safely stable.

Verify: true $y_1(0.2)=e^{-4}=0.0183$ lekin RK4 ne $5.0$ diya — ek catastrophic overshoot (factor $\sim 270$ se galat, aur $|R(-4)|=5>1$ to yeh har step grow karega). $y_2$ ne $0.8187$ diya vs true $e^{-0.2}=0.8187$ — perfect. Fix: $h$ chhota karo taaki stiffest rate ke liye $|h\lambda |<2.78$ ho, matlab $h<2.78/20=0.139$, ya ek implicit method use karo. Explicit RK4 sabse fast component se limited hota hai.

---

Ex 8 — Cell H: non-autonomous (explicit $t$)

Forecast: ab $\mathbf{f}$ explicitly $t$ mention karta hai, to $t_n+\frac{h}{2}$ aur $t_n+h$ arguments actually matter karte hain — yeh woh case hai jahan time advance karna bhool jaane se sab galat ho jaata hai. Start slope $y_2^{\prime }=0-0=0$, $y_1^{\prime }=1$: position rise kar rahi hai, velocity flat. Kyunki exact answer straight line $y=t$ hai, RK4 ko $0.2$ aur $1$ essentially zero error ke saath return karna chahiye.

Yahan $\mathbf{f}(t,\mathbf{y})=(y_2,t-y_1)$, state $(0,1)$ at $t_0=0$.

1. ${\mathbf{k}}_1=(1,0-0)=(1,0)$. Kyun? $t=0$, state $(0,1)$ pe slopes.
2. $t_0+\frac{h}{2}=0.1$, state $=(0,1)+0.1(1,0)=(0.1,1)$; ${\mathbf{k}}_2=(1,0.1-0.1)=(1,0)$. Kyun? Midpoint — note karo ki humne $t=0.1$ use kiya, $0$ nahi.
3. State $=(0,1)+0.1(1,0)=(0.1,1)$; ${\mathbf{k}}_3=(1,0.1-0.1)=(1,0)$. Kyun? $t=0.1$ pe refined midpoint.
4. $t_0+h=0.2$, state $=(0,1)+0.2(1,0)=(0.2,1)$; ${\mathbf{k}}_4=(1,0.2-0.2)=(1,0)$. Kyun? End, $t=0.2$ pe.
5. Combine: $y_{1,1}=0+\frac{0.2}{6}(1+2(1)+2(1)+1)=\frac{0.2}{6}(6)=\mathrm{0.2.}$ $y_{2,1}=1+\frac{0.2}{6}(0+0+0+0)=1.$ Yeh step kyun? Charon slopes ka $1,2,2,1$ blend average velocity $(1,0)$ deta hai, to hum $h\cdot (1,0)=(0.2,0)$ se step karte hain — lekin sirf isliye kyunki humne har ${\mathbf{k}}_i$ ko correct time pe evaluate kiya.

Verify: exact solution $y(t)=t$ hai (check: $y^{\prime \prime }=0$, to $y^{\prime \prime }+y=0+t=t$ ✓; $y(0)=0$, $y^{\prime }(0)=1$ ✓), jisse $y(0.2)=0.2$ aur $y^{\prime }(0.2)=1$ milta hai. RK4 ne exactly $y_{1,1}=0.2$ aur $y_{2,1}=1$ return kiya — zero error. Yeh isliye hua kyunki humne har stage mein $t$ sahi se advance kiya; agar hum $t=0$ freeze karte, to $t-y_1$ term galat time pe evaluate hota aur answer galat hota.

---

Recall

Recall Quick self-test

Ex 5 exactly apni jagah kyun raha? ::: Fixed point pe $\mathbf{f}=0$ hai, to charon ${\mathbf{k}}_i$ vanish ho jaate hain aur weighted average kuch add nahi karta. Ex 7 ne $0$ ke near ki bajaye $5.0$ kyun return kiya? ::: RK4 amplification $R(h\lambda )$ at $h\lambda =-4$ equals $5>1$ hai, to step decay ki bajaye amplify karta hai — hum stability region ke bahar the. Ex 8 mein agar aap har stage mein $t=0$ rakhte to kya break hota? ::: $t-y_1$ term galat time pe evaluate hota, galat $\mathbf{k}$'s milte aur exact answer kho jaata.

---

Connections

• Parent: RK4 for systems
• RK4 for a single ODE — Ex 2 dikhata hai ki decoupled case isme reduce ho jaata hai.
• Reducing higher-order ODEs to first-order systems — Ex 3, Ex 8.
• Euler's method for systems — sasta cousin; agar accuracy gap dekhna ho to in examples ko iske saath try karo.
• Local vs Global Truncation Error — kyun Ex 2 mein bade rates ne bada error diya.
• Stiff systems and stability — Ex 7 mein failure.
• Simpson's Rule — woh quadrature jo RK4 imitate karta hai.