Exercises — Implementing ODE solvers from scratch — Euler, RK4

Euler slope ko sirf step ke left endpoint par evaluate karta hai, point $(t_n,y_n)$ par, aur phir maan leta hai ki woh ek slope poore interval $[t_n,t_n+h]$ mein constant rehta hai. $y_{n+1}=y_n+hf(t_n,y_n)$ Socho ek straight line $(t_n,y_n)$ se us slope ke saath nikal rahi hai; tum us par horizontal distance $h$ slide karte ho aur nayi height $y_{n+1}$ padhte ho.

Start par slope: $f(t_0,y_0)=f(0,0)=2\cdot 0=0$. $y_1=y_0+hf(0,0)=0+0.5\cdot 0=0$ Yahan $f$ sirf $t$ par depend karta hai, aur $t=0$ par slope genuinely $0$ hai, toh Euler pehle step mein koi rise predict nahi karta (true solution $y=t^2$ bhi origin par flat hai). $y_1=0$.

Step 1 ($t=0$ se): slope $=f(0,1)=1$. $y_1=1+0.5\cdot 1=1.5(t_1=0.5)$ Step 2 ($t=0.5$ se): slope $=f(0.5,1.5)=1.5$. $y_2=1.5+0.5\cdot 1.5=2.25(t_2=1.0)$ Error $=2.71828-2.25=0.46828$. Curve $e^t$ upar ki taraf bend karta hai, toh Euler — har baar lower left-endpoint slope use karke — undershoot karta hai.

Yahan $f(t,y)=\sin t$ hai, toh har $k_i$ sirf shifted time ka $\sin$ hai:

Euler ka global error $O(h)$ hai — $h$ ko adha karne se error lagbhag adha ho jaata hai (ratio $0.125/0.065\approx 1.9\approx 2$ yeh confirm karta hai). $0.05$ se $0.025$ tak dobara adha karte hue: $\text{error}(0.025)\approx \frac{0.065}{2}\approx 0.0325$ $O(h)$ kyun? Har step ek $O(h^2)$ size ka Taylor term phek deta hai (local truncation error). $N=T/h$ steps mein yeh accumulate hote hain: $N\cdot O(h^2)=\frac{T}{h}\cdot O(h^2)=O(h)$. Dekho Taylor Series Expansion.

RK4 ka global error $O(h^4)$ hai, toh $h$ adha karne par error $2^4=16$ se ghatta hai: $\text{error}(0.05)\approx \frac{2.1\times 10^{-6}}{16}\approx 1.31\times 10^{-7}$ $\text{error}(0.025)\approx \frac{1.31\times 10^{-7}}{16}\approx 8.2\times 10^{-9}$ Yahi wajah hai ki RK4 dominate karta hai: per step same effort seize karne par accuracy sixteen-fold badhti hai har halving mein, jabki Euler ka sirf doubling hota hai.

Euler left endpoint par khinchi tangent line follow karta hai. Aisi curve ke liye jo concave up hai (upar ki taraf bend karti hai, jaise $e^t$), tangent curve ke neeche rehti hai, toh Euler undershoot karta hai (jaise L2.1 mein: $2.25<2.718$). Aisi curve ke liye jo concave down hai, tangent curve ke upar rehti hai, toh Euler overshoot karta hai. Error ka sign $y^{\prime \prime }$ ke sign se set hota hai — precisely woh $\frac{h^2}{2}{y}^{\prime \prime }$ Taylor term jo Euler discard karta hai. RK4 ise largely cancel karta hai kyunki woh commit karne se pehle midpoint aur end par peek karta hai.

Steps ki sankhya $N=1/h$, aur error $\sim C{h}^p$, toh $\epsilon$ tak pahunchne ke liye $h\sim (\epsilon /C{)}^{1/p}$ chahiye, yani $N\sim {\epsilon }^{-1/p}$ (constant $C$ ignore karte hue).

• Euler $p=1$: $N\sim {\epsilon }^{-1}=10^8$ steps — ek sau million.
• RK4 $p=4$: $N\sim {\epsilon }^{-1/4}=10^2=100$ steps. Bhalecha RK4 4 function calls per step cost karta hai, uska total work $\approx 4\times 100=400$ evaluations Euler ke $10^8$ ko crush karta hai. Verdict: kisi bhi demanding accuracy par, RK4 orders of magnitude se jeetta hai — yahi poori wajah hai ki yeh default workhorse hai.

$f(t,y)=\lambda y$ ke saath ($z=\lambda h$ likho):

Har step $R(z)$ se multiply karta hai, toh $n$ steps baad $y_n=R(z{)}^n{y}_0$. Numerical solution ke liye true solution ki tarah decay karne ke liye ($e^{\lambda t}\to 0$ jab $\lambda <0$), hamen chahiye $|R(z)|\le 1.$ Real $z=\lambda h<0$ ke liye yeh tab tak hold karta hai jab $z\gtrsim -2.79$ (RK4 ki real stability boundary, sabse negative $z$ jahan $|R(z)|=1$). Toh strongly negative $\lambda$ ke liye tumhe $h\lesssim 2.79/|\lambda |$ rakhna padega — ek bahut chhota step — ya phir numbers blow up ho jaate hain bhalecha true answer shrink kar raha ho. Yeh Numerical Stability and Stiff ODEs ka darwaza hai: stiff problems mein huge $|\lambda |$ hota hai, jo tiny $h$ force karta hai aur implicit solvers ko motivate karta hai.

$\mathbf{y}$ aur har $k_i$ ko vectors ki tarah treat karo; $f$ ek vector return karta hai. Python mein NumPy arrays ke saath code identical hai — arithmetic y + h/2 * k1, weighted sum, sab componentwise broadcast hota hai:

def rk4_system(f, t0, y0, h, n):   # y0 is a numpy array
    t, y = t0, y0.copy()
    ts, ys = [t], [y.copy()]
    for _ in range(n):
        k1 = f(t,       y)
        k2 = f(t + h/2, y + h/2 * k1)
        k3 = f(t + h/2, y + h/2 * k2)
        k4 = f(t + h,   y + h   * k3)
        y = y + (h/6) * (k1 + 2*k2 + 2*k3 + k4)
        t = t + h
        ts.append(t); ys.append(y.copy())
    return ts, ys

Update $k_i$ mein linear hai, toh yeh har component par independently apply hota hai — wahi weighted average of slopes, ab ek slope per dimension. Isi tarah 2nd-order ODE $y^{\prime \prime }=g$ solve hota hai: system $\mathbf{y}=(y,v)$, ${\mathbf{y}}^{\prime }=(v,g)$ ke roop mein rewrite karo.

$\mathbf{y}=(y,v)$ lo jahan $v=y^{\prime }$. Toh ${\mathbf{y}}^{\prime }=(v,-y)$, start ${\mathbf{y}}_0=(1,0)$, $h=\frac{\pi }{2}\approx 1.5708$. Har stage ek vector hai; hum dono components track karte hain aur end mein position padhte hain.

• (a) $h=0$: har update $y_{n+1}=y_n+0=y_n$ hoga aur $t$ kabhi advance nahi karta. Loop chalta hai lekin $t_0$ par constant list produce karta hai — mathematically vacuous (koi time pass nahi hota). Code mein iske against guard lagao.
• (b) $n=0$: for loop body kabhi execute nahi hoti; tum sirf initial point $([t_0],[y_0])$ return karte ho. Correct — zero steps maangne par.
• (c) $f\equiv 0$: saare $k_i=0$, toh $y_{n+1}=y_n$ hamesha — exactly true constant solution. Euler aur RK4 dono yahan exact hain (zero error), kyunki ek constant degree-0 polynomial hai jo RK4 ki order-4 exactness ke andar aata hai.
• (d) discontinuous $f$: RK4 ki order-4 accuracy assume karti hai ki $f$ smooth hai (yeh $y$ ki Taylor expansion par rely karti hai). $[t_n,t_n+h]$ ke andar ek jump us assumption ko todta hai; tum silently accuracy kho dete ho. Fix: grid point ko exactly discontinuity par rakho taaki koi step ise straddle na kare — wahi event-handling jo scipy.integrate.solve_ivp provide karta hai.