Exercises — Systems of ODEs — RK4 for systems

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

Throughout, $t_n$ se $t_{n+1}=t_n+h$ tak ek RK4 step mein woh recipe use hoti hai jo tumhare paas already hai:

\mathbf{k}_2=\mathbf f\!\left(t_n+\tfrac h2,\ \mathbf y_n+\tfrac h2\mathbf k_1\right),$$ $$\mathbf{k}_3=\mathbf f\!\left(t_n+\tfrac h2,\ \mathbf y_n+\tfrac h2\mathbf k_2\right),\quad \mathbf{k}_4=\mathbf f\!\left(t_n+h,\ \mathbf y_n+h\,\mathbf k_3\right),$$ $$\mathbf y_{n+1}=\mathbf y_n+\frac h6\big(\mathbf k_1+2\mathbf k_2+2\mathbf k_3+\mathbf k_4\big).$$ Yahan **bold** ka matlab hai "numbers ki ek list" (ek vector); $\mathbf y_n=(y_{1,n},\,y_{2,n},\dots)$ step $n$ par har unknown ka hamaara estimate hai. Har $\mathbf k_i$ **usi length ki ek list** hai. **Neeche wali Figure — Exercise 2.1 se pehle ise padho.** Yeh picture ek single RK4 step ka *map* hai: yeh exactly dikhati hai ki time axis par **kahan** chaar slope samples $\mathbf k_1,\mathbf k_2,\mathbf k_3,\mathbf k_4$ liye jaate hain, aur yaad dilati hai ki inhe weights $1,2,2,1$ (÷6) se combine kiya jaata hai. $\mathbf k_1$ start $t_n$ par hai, $\mathbf k_2$ aur $\mathbf k_3$ dono midpoint $t_n+h/2$ par hain (double weight isliye hai kyunki midpoint sabse zyada count karta hai), aur $\mathbf k_4$ end $t_n+h$ par hai. Ise dhyan mein rakho: is page par har numeric exercise bas "in chaar dots ko evaluate karo, phir inhe weight-average karo" hai. ![[deepdives/dd-maths-4.8.27-d4-s01.png]] > [!definition] Figure caption > **Time axis par dekha gaya ek RK4 step.** Yellow dot = $t_n$ par start slope $\mathbf k_1$; blue dot = pehla midpoint slope $\mathbf k_2$; pink dot = refined midpoint slope $\mathbf k_3$ (dono $t_n+h/2$ par); yellow dot = $t_n+h$ par end slope $\mathbf k_4$. Final $\tfrac h6(1,2,2,1)$ average mein dono midpoint dots double weight carry karte hain. --- ## L1 — Recognition ### Exercise 1.1 Neeche wala system component-by-component likha gaya hai. Ise ek single **vector** equation $\mathbf y'=\mathbf f(t,\mathbf y)$ ke roop mein likho aur vector ki length $n$ batao. $$y_1'=3y_1-2y_2,\qquad y_2'=y_1+t.$$ > [!recall]- Solution > Unknowns ko $\mathbf y=(y_1,y_2)$ mein stack karo, toh $n=\boxed{2}$. > Right-hand sides $\mathbf f$ ke components ban jaate hain: > $$\mathbf f(t,\mathbf y)=\big(\,3y_1-2y_2,\ \ y_1+t\,\big).$$ > **Humne kya kiya:** bas group kiya. **Kyun:** RK4-for-systems poore vector par ek saath run karta hai, isliye hume ise ek object ki tarah dekhna hoga. Doosre component mein $t$ allowed hai — $\mathbf f$ explicitly time par depend kar sakta hai. ### Exercise 1.2 RK4 recipe mein, kaun sa stage argument $\mathbf y_n+\tfrac h2\mathbf k_2$ hai aur wahan $\mathbf f$ kis time par evaluate hota hai? > [!recall]- Solution > Yeh $\mathbf k_3$ hai, time $t_n+\tfrac h2$ (wahi midpoint) par evaluate hota hai. Yeh ek **refined** midpoint slope hai: $\mathbf k_2$ ke saath same time, lekin $\mathbf k_2$ se build ki gayi better intermediate state use karta hai. ### Exercise 1.3 2nd-order ODE $y''-5y'+6y=0$ ko derivatives ko naam dekar first-order system mein convert karo. > [!recall]- Solution > $u_1=y$ aur $u_2=y'$ set karo. Tab $u_1'=u_2$ ($u_2$ ki definition se), aur ODE se $y''=5y'-6y=5u_2-6u_1$: > $$u_1'=u_2,\qquad u_2'=5u_2-6u_1.$$ > **Yeh kyun kaam karta hai:** ek 2nd-order equation mein do independent pieces of information hoti hain ($y$ aur $y'$), isliye yeh do 1st-order equations ke equivalent hai. Dekho [[Reducing higher-order ODEs to first-order systems]]. > [!mistake] L1 trap — "$\mathbf f$ mein $t$ ka matlab hai yeh system nahi hai." > **Kyun sahi lagta hai:** hum systems ko sirf predator–prey couplings ke roop mein sochte hain jisme koi clock nahi hota. **Kyun galat hai:** $\mathbf f(t,\mathbf y)$ ko $t$ par directly depend karne ki *permission* hai; isliye toh definition mein $t$ appear karta hai. **Fix:** time ko ek ordinary input ki tarah treat karo jise tum plug in karte ho — $\mathbf k_2,\mathbf k_3$ mein $t_n+\tfrac h2$ use karo, $\mathbf k_4$ mein $t_n+h$ use karo. --- ## L2 — Application ### Exercise 2.1 $\mathbf f(t,\mathbf y)=(y_2,\,-y_1)$ ke liye $\mathbf y_0=(1,0)$ aur $h=0.2$ ke saath, **sirf** $\mathbf k_1$ aur $\mathbf k_2$ compute karo. (Upar wali figure use karo: $\mathbf k_1$ step ka start-of-step dot hai, $\mathbf k_2$ pehla midpoint dot hai.) > [!recall]- Solution > $\mathbf k_1=\mathbf f(t_0,\mathbf y_0)=(y_2,-y_1)=(0,\,-1)$. > Intermediate state: $\mathbf y_0+\tfrac h2\mathbf k_1=(1,0)+0.1(0,-1)=(1,\,-0.1)$. > $\mathbf k_2=\mathbf f(\cdot,(1,-0.1))=(-0.1,\,-1)$. > **Figure ka blue dot dekho:** woh $t_0+h/2$ par hai, aur ise build karne ke liye dono components ko *saath mein* advance kiya gaya — pehle slot mein $-0.1$ use hua, na ki stale $0$. ### Exercise 2.2 Exercise 2.1 ke system ka **ek poora RK4 step** lo ($h=0.2$). $\mathbf y_1$ ko 5 decimals tak do, phir $t=0.2$ par true solution $(\cos t,-\sin t)$ se compare karo. > [!recall]- Solution > $\mathbf k_1=(0,-1)$. > $\mathbf k_2=(-0.1,-1)$ (2.1 se). > $\mathbf k_3$ ke liye state: $(1,0)+0.1(-0.1,-1)=(0.99,\,-0.1)$, toh $\mathbf k_3=(-0.1,\,-0.99)$. > $\mathbf k_4$ ke liye state: $(1,0)+0.2(-0.1,-0.99)=(0.98,\,-0.198)$, toh $\mathbf k_4=(-0.198,\,-0.98)$. > Weights $1,2,2,1$ se combine karo: > $$y_{1,1}=1+\tfrac{0.2}{6}\big(0+2(-0.1)+2(-0.1)+(-0.198)\big)=1+\tfrac{0.2}{6}(-0.598)=0.980067.$$ > $$y_{2,1}=0+\tfrac{0.2}{6}\big(-1+2(-1)+2(-0.99)+(-0.98)\big)=\tfrac{0.2}{6}(-5.96)=-0.198667.$$ > True: $\cos 0.2=0.980067$, $-\sin 0.2=-0.198669$. ~5 aur ~5 decimals tak match karta hai — ek hi step mein! ### Exercise 2.3 $y'=y$ ($f=y$, scalar) ke liye, $y_0=1$ se $h=0.5$ ke saath ek RK4 step lo. $e^{0.5}$ se compare karo. > [!recall]- Solution > $k_1=1$. > $k_2=y_0+\tfrac h2 k_1=1+0.25(1)=1.25$. > $k_3=1+0.25(1.25)=1.3125$. > $k_4=1+0.5(1.3125)=1.65625$. > Bracket sum $=k_1+2k_2+2k_3+k_4=1+2(1.25)+2(1.3125)+1.65625=1+2.5+2.625+1.65625=7.78125$. > $$y_1=1+\tfrac{0.5}{6}(7.78125)=1+0.648438=1.648438.$$ > True $e^{0.5}=1.648721$; absolute error $\approx 2.8\times10^{-4}$. Order-4 method hai, $h=0.5$ par bhi tiny error. > [!mistake] L2 trap — ek component ko saare chaar stages se advance karna doosre ko touch karne se pehle. > **Kyun sahi lagta hai:** har equation "alag" lagti hai, toh tum pehle $y_1$ ka poora RK4 khatam karna chahte ho. **Kyun galat hai:** component 1 ka $\mathbf k_2$ $y_2$ ki *advanced* value chahiye — agar $y_2$ abhi bhi stale hai, toh tum purana data feed kar rahe ho aur accuracy lose karte ho (ya worse, galat answer). **Fix:** poora vector $\mathbf k_1$ compute karo, phir poori intermediate state banao, phir poora $\mathbf k_2$. Stages components ke across synchronized hote hain. --- ## L3 — Analysis ### Exercise 3.1 Linear system $\mathbf y'=A\mathbf y$ ke liye $A=\begin{pmatrix}0&1\\-1&0\end{pmatrix}$ aur $\mathbf y_0=(1,0)$ ke saath, true energy $E=y_1^2+y_2^2$ exactly conserved hai ($E=1$ sabhi $t$ ke liye). **Ek** RK4 step ke baad $h=0.2$ ke saath (Exercise 2.2 ne $\mathbf y_1=(0.980067,-0.198667)$ diya), $E_1=y_{1,1}^2+y_{2,1}^2$ compute karo. Kya energy RK4 se exactly conserved hoti hai? Explain karo. > [!recall]- Solution > Exact rational one-step result use karke, $E_1\approx 0.99999911$ — yaani $E_1$, $1$ se lagbhag $9\times10^{-7}$ differ karta hai. > Energy **almost** conserved hai lekin exactly nahi: RK4 ka step map true rotation ka polynomial approximation hai, toh yeh $E$ ko slightly change karta hai. Drift yahan tiny hai ($\sim10^{-6}$) kyunki local error $O(h^5)$ hai. *Bahut saare* steps ke baad yeh drift accumulate hota hai. Jab aise structure-preservation matter karta hai, dekho [[Stiff systems and stability]]. ### Exercise 3.2 Scalar test problem $y'=\lambda y$ mein ek single RK4 step $y$ ko ek fixed factor $R(h\lambda)$ se multiply karta hai. $R(z)$ derive karo jahan $z=h\lambda$. > [!recall]- Solution > $f=\lambda y$ ke saath: $k_1=\lambda y$. Phir $k_2=\lambda(y+\tfrac h2\lambda y)=\lambda y(1+\tfrac z2)$. > $k_3=\lambda(y+\tfrac h2 k_2)=\lambda y\big(1+\tfrac z2+\tfrac{z^2}{4}\big)$. > $k_4=\lambda(y+h k_3)=\lambda y\big(1+z+\tfrac{z^2}{2}+\tfrac{z^3}{4}\big)$. > Ab $y_1=y+\tfrac h6(k_1+2k_2+2k_3+k_4)$. $z$ ki powers collect karke: > $$R(z)=1+z+\frac{z^2}{2}+\frac{z^3}{6}+\frac{z^4}{24}.$$ > **Iska matlab:** yeh exactly $e^z$ ke pehle 5 terms hain — RK4 true growth $e^{h\lambda}$ ko $z^4$ tak reproduce karta hai, aur *wahi* jagah se "order 4" aata hai. Dekho [[Local vs Global Truncation Error]]. ### Exercise 3.3 3.2 se $R(z)$ use karke $z=0.5$ ke saath, confirm karo ki yeh Exercise 2.3 ka answer reproduce karta hai. > [!recall]- Solution > $R(0.5)=1+0.5+\tfrac{0.25}{2}+\tfrac{0.125}{6}+\tfrac{0.0625}{24}=1+0.5+0.125+0.0208333+0.0026042=1.6484375$. > 2.3 se $y_1=1.648438$ ke identical. **Kyun:** ek single scalar RK4 step *hi* $R(z)$ se multiplication hai, toh dono computations exactly agree karni chahiye. > [!mistake] L3 trap — "order 4 ka matlab answer exact hai." > **Kyun sahi lagta hai:** chaar Taylor terms match karna sab kuch match karne jaisa lagta hai. **Kyun galat hai:** $R(z)$, $z^4$ par ruk jaata hai; true $e^z$ mein $\tfrac{z^5}{120}+\cdots$ hai jo RK4 *miss* karta hai. Woh missing tail hi $O(h^5)$ local error hai. **Fix:** RK4 ko "exponential ke 4 terms tak correct" samjho — excellent hai lekin exact nahi. --- ## L4 — Synthesis ### Exercise 4.1 $y''+0.4y'+y=0$, $y(0)=2$, $y'(0)=0$ ko first-order system mein reduce karo, phir $h=0.1$ ke saath **ek** RK4 step lo. $\mathbf y_1=(u_{1,1},u_{2,1})$ ko 5 decimals tak report karo. > [!recall]- Solution > **Reduce.** $u_1=y,\ u_2=y'$. $y''=-0.4y'-y=-0.4u_2-u_1$ se: > $$\mathbf f(u_1,u_2)=\big(u_2,\ -u_1-0.4u_2\big),\qquad \mathbf y_0=(2,0).$$ > **k1:** $\mathbf f(2,0)=(0,\,-2)$. > **k2 state:** $(2,0)+0.05(0,-2)=(2,\,-0.1)$; $\mathbf f=( -0.1,\ -2+0.04)=(-0.1,\,-1.96)$. > **k3 state:** $(2,0)+0.05(-0.1,-1.96)=(1.995,\,-0.098)$; $\mathbf f=(-0.098,\ -1.995+0.0392)=(-0.098,\,-1.9558)$. > **k4 state:** $(2,0)+0.1(-0.098,-1.9558)=(1.9902,\,-0.19558)$; $\mathbf f=(-0.19558,\ -1.9902+0.078232)=(-0.19558,\,-1.911968)$. > **Combine.** > $$u_{1,1}=2+\tfrac{0.1}{6}\big(0+2(-0.1)+2(-0.098)+(-0.19558)\big)=2+\tfrac{0.1}{6}(-0.59158)=1.990140.$$ > $u_2$ ke liye bracket sum: $-2+2(-1.96)+2(-1.9558)+(-1.911968)=-2-3.92-3.9116-1.911968=-11.743568$. > $$u_{2,1}=0+\tfrac{0.1}{6}(-11.743568)=-0.195726.$$ > Toh $\mathbf y_1\approx(1.99014,\,-0.19573)$. Negative $u_2$ (velocity) dikhata hai ki mass wapas swing karna shuru ho gaya hai; $0.4u_2$ term friction hai — yeh ek **damped** oscillator hai, frictionless Example 1 se alag. ### Exercise 4.2 Coupled system $y_1'=y_1-y_1 y_2,\ \ y_2'=-y_2+y_1 y_2$ (ek Lotka–Volterra style model) ke liye $\mathbf y_0=(1,1)$ ke saath, $\mathbf k_1$ compute karo aur explain karo ki yeh $t=0$ par motion ke baare mein kya kehta hai. > [!recall]- Solution > $\mathbf k_1=\big(y_1-y_1y_2,\ -y_2+y_1y_2\big)=\big(1-1,\ -1+1\big)=(0,\,0)$. > **Interpretation:** $(1,1)$ ek **equilibrium** hai — dono slopes zero hain, toh first order mein kuch move nahi karta. RK4 $(1,1)$ ke bahut close $\mathbf y_1$ predict karega; tiny drift sirf midpoint stages par feel hone wali nonlinearity se aata hai. Yeh ek **degenerate/fixed-point** case hai jise tumhe recognize karna chahiye: zero $\mathbf k_1$ ka matlab yeh nahi ki tum done ho, lekin iska matlab hai ki start-slope kuch contribute nahi karta. > [!mistake] L4 trap — ek driven ODE reduce karte waqt stages ke beech *time* freeze karna. > **Kyun sahi lagta hai:** Example 4.1 mein system mein koi explicit $t$ nahi tha, toh beginners assume karte hain ki $t$ kabhi matter nahi karta. **Kyun galat hai:** agar ODE $y''+y=\cos t$ hai (ek forced oscillator), toh $\mathbf f$, $t$ par depend karta hai; phir $\mathbf k_2,\mathbf k_3$ mein $\cos(t_n+\tfrac h2)$ use karna hoga aur $\mathbf k_4$ mein $\cos(t_n+h)$ use karna hoga. Har jagah $\cos t_n$ use karne se silently accuracy drop ho jaati hai. **Fix:** har stage par DONO state *aur* time slot advance karo. --- ## L5 — Mastery ### Exercise 5.1 Linear air drag wala ek projectile 2D mein $\mathbf r''=-g\,\hat{\mathbf j}-c\,\mathbf r'$ obey karta hai. Chaar unknowns $(x,y,v_x,v_y)$ mein **poora first-order system** likho, $g=9.8$, $c=0.1$ ke saath. Phir, $\mathbf y_0=(0,0,10,10)$ aur $h=0.1$ ke liye, $\mathbf k_1$ aur $\mathbf k_2$ compute karo (chaar components mein se har ek). > [!recall]- Solution > **Reduce.** Position $\to$ velocity, velocity $\to$ acceleration: > $$x'=v_x,\quad y'=v_y,\quad v_x'=-c\,v_x,\quad v_y'=-g-c\,v_y.$$ > Toh $\mathbf f=(v_x,\ v_y,\ -0.1v_x,\ -9.8-0.1v_y)$, ek $n=4$ vector hai. > **k1** at $(0,0,10,10)$: > $$\mathbf k_1=(10,\ 10,\ -0.1(10),\ -9.8-0.1(10))=(10,\ 10,\ -1,\ -10.8).$$ > **k2 state** $=\mathbf y_0+0.05\,\mathbf k_1=(0.5,\ 0.5,\ 9.95,\ 9.46)$: > $$\mathbf k_2=\big(9.95,\ 9.46,\ -0.1(9.95),\ -9.8-0.1(9.46)\big)=(9.95,\ 9.46,\ -0.995,\ -10.746).$$ > **Padho:** components 1–2 *velocities* hain (yeh positions ko feed karte hain), components 3–4 *accelerations* hain. $\mathbf k_2$ se pehle charon saath mein move kiye — poore scale par synchronization rule. ### Exercise 5.2 Exercise 5.1 ka ek poora RK4 step complete karo aur $\mathbf y_1=(x_1,y_1,v_{x,1},v_{y,1})$ ko 5 decimals tak report karo. > [!recall]- Solution > 5.1 se continue karo. > **k3 state** $=\mathbf y_0+0.05\,\mathbf k_2=(0.4975,\ 0.473,\ 9.95025,\ 9.4627)$: > $$\mathbf k_3=(9.95025,\ 9.4627,\ -0.995025,\ -10.74627).$$ > **k4 state** $=\mathbf y_0+0.1\,\mathbf k_3=(0.995025,\ 0.94627,\ 9.9004975,\ 8.925373)$: > $$\mathbf k_4=(9.9004975,\ 8.925373,\ -0.99004975,\ -10.6925373).$$ > **Combine** $\tfrac{0.1}{6}$ aur weights $1,2,2,1$ ke saath. $x_1$ ke liye bracket sum hai > $10+2(9.95)+2(9.95025)+9.9004975=59.7009975$, toh > $$x_1=0+\tfrac{0.1}{6}(59.7009975)=0.995017.$$ > $y_1$ ke liye bracket sum hai $10+2(9.46)+2(9.4627)+8.925373=56.770873$, toh > $$y_1=0+\tfrac{0.1}{6}(56.770873)=0.946181.$$ > Velocities ke liye: > $$v_{x,1}=10+\tfrac{0.1}{6}\big(-1+2(-0.995)+2(-0.995025)+(-0.99004975)\big)=9.900498.$$ > $$v_{y,1}=10+\tfrac{0.1}{6}\big(-10.8+2(-10.746)+2(-10.74627)+(-10.6925373)\big)=8.925362.$$ > Toh $\mathbf y_1\approx(0.99502,\ 0.94618,\ 9.90050,\ 8.92536)$. Gravity ne already $v_y$ ko $10$ se $8.93$ ki taraf pull kar liya hai, aur drag ne $v_x$ ko $10$ se thoda neeche nudge kiya — physically correct hai. ### Exercise 5.3 3.2 se $R(z)=1+z+\tfrac{z^2}{2}+\tfrac{z^3}{6}+\tfrac{z^4}{24}$ use karke explain karo ki drag component $v_x'=-0.1v_x$ par large step $h$ ke saath RK4 kyun stable rehta hai, lekin **stiff** eigenvalue jaise $\lambda=-1000$ ke liye tumhe $|h\lambda|$ small kyun rakhna padega. > [!recall]- Solution > Ek component ki stability ke liye $|R(h\lambda)|\le 1$ chahiye. Drag component ke liye $\lambda=-0.1$, $h=1$ bhi lene par $z=-0.1$ milta hai aur $R(-0.1)\approx0.90484$, jo safely $[-1,1]$ ke andar hai — stable. > $\lambda=-1000$ ke liye, koi bhi moderate $h$ lene par (jaise $h=0.01$) $z=-10$ milta hai; phir $R(-10)=1-10+50-\tfrac{1000}{6}+\tfrac{10000}{24}=1-10+50-166.6\overline{6}+416.6\overline{6}=291$, jiska magnitude $\gg 1$ hai — numerical solution **explode** kar jaata hai jabki true solution decay karta hai. Stable rehne ke liye $h$ ko itna chhota karna hoga ki $|z|$ small ho (real negative $z$ ke liye roughly $|z|\lesssim 2.78$). Yahi forced tiny step **stiffness** ki definition hai. Dekho [[Stiff systems and stability]]. > [!mistake] L5 trap — "chhota error matlab chhota step safe hai; bada $h$ sirf accuracy cost karta hai." > **Kyun sahi lagta hai:** smooth non-stiff problems ke liye, bada $h$ sirf $O(h^4)$ error raise karta hai. **Kyun galat hai:** ek stiff eigenvalue ke liye, jab ek baar $|R(h\lambda)|>1$ ho jaata hai toh error sirf grow nahi karta — woh **geometrically blow up** karta hai, ek decaying true solution ko numerical garbage mein badal deta hai. **Fix:** stiff systems ke liye *stability* condition $|R(h\lambda)|\le1$ check karo, sirf accuracy order nahi; ya phir implicit method pe switch karo. --- ## Connections - [[4.8.27 Systems of ODEs — RK4 for systems (Hinglish)|Parent topic — RK4 for systems]] - [[RK4 for a single ODE]] — yahan har exercise har component par isi mein reduce hoti hai. - [[Reducing higher-order ODEs to first-order systems]] — 1.3, 4.1, 5.1 mein use hua. - [[Euler's method for systems]] — compare karo ki tum kitni accuracy kho dete. - [[Local vs Global Truncation Error]] — L3 mein $R(z)$ ki kahani. - [[Stiff systems and stability]] — L5.3. - [[Simpson's Rule]] — woh quadrature jise weights $1,2,2,1$ imitate karte hain.