Exercises — Systems of ODEs — RK4 for systems

Throughout, one RK4 step from $t_n$ to $t_{n+1}=t_n+h$ uses the recipe you already own:

\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).$$ Here **bold** means "a list of numbers" (a vector); $\mathbf y_n=(y_{1,n},\,y_{2,n},\dots)$ is our estimate of every unknown at step $n$. Each $\mathbf k_i$ is a **list of the same length**. **Figure below — read it before Exercise 2.1.** The picture is the *map of a single RK4 step*: it shows exactly **where on the time axis** the four slope samples $\mathbf k_1,\mathbf k_2,\mathbf k_3,\mathbf k_4$ are taken, and reminds you that they are combined with the weights $1,2,2,1$ (÷6). $\mathbf k_1$ sits at the start $t_n$, $\mathbf k_2$ and $\mathbf k_3$ both sit at the midpoint $t_n+h/2$ (the double weight is why the midpoint counts most), and $\mathbf k_4$ sits at the end $t_n+h$. Keep this in view: every numeric exercise on this page is just "evaluate these four dots, then weight-average them." ![[deepdives/dd-maths-4.8.27-d4-s01.png]] > [!definition] Figure caption > **One RK4 step, viewed on the time axis.** Yellow dot = start slope $\mathbf k_1$ at $t_n$; blue dot = first midpoint slope $\mathbf k_2$; pink dot = refined midpoint slope $\mathbf k_3$ (both at $t_n+h/2$); yellow dot = end slope $\mathbf k_4$ at $t_n+h$. The two midpoint dots carry double weight in the final $\tfrac h6(1,2,2,1)$ average. --- ## L1 — Recognition ### Exercise 1.1 The system below is written out component-by-component. Write it as a single **vector** equation $\mathbf y'=\mathbf f(t,\mathbf y)$ and state the length $n$ of the vector. $$y_1'=3y_1-2y_2,\qquad y_2'=y_1+t.$$ > [!recall]- Solution > Stack the unknowns into $\mathbf y=(y_1,y_2)$, so $n=\boxed{2}$. > The right-hand sides become the components of $\mathbf f$: > $$\mathbf f(t,\mathbf y)=\big(\,3y_1-2y_2,\ \ y_1+t\,\big).$$ > **What we did:** just grouped. **Why:** RK4-for-systems runs on the whole vector at once, so we must see it as one object. The $t$ in the second component is allowed — $\mathbf f$ may depend on time explicitly. ### Exercise 1.2 In the RK4 recipe, which stage argument is $\mathbf y_n+\tfrac h2\mathbf k_2$ and at what time is $\mathbf f$ evaluated there? > [!recall]- Solution > That is $\mathbf k_3$, evaluated at time $t_n+\tfrac h2$ (the midpoint again). It is a **refined** midpoint slope: same time as $\mathbf k_2$ but using the better intermediate state built from $\mathbf k_2$. ### Exercise 1.3 Convert the 2nd-order ODE $y''-5y'+6y=0$ into a first-order system by naming the derivatives. > [!recall]- Solution > Set $u_1=y$ and $u_2=y'$. Then $u_1'=u_2$ (by definition of $u_2$), and from the ODE $y''=5y'-6y=5u_2-6u_1$: > $$u_1'=u_2,\qquad u_2'=5u_2-6u_1.$$ > **Why this works:** one 2nd-order equation carries two independent pieces of information ($y$ and $y'$), so it is equivalent to two 1st-order equations. See [[Reducing higher-order ODEs to first-order systems]]. > [!mistake] L1 trap — "$t$ in $\mathbf f$ means it isn't a system." > **Why it feels right:** we picture systems as purely predator–prey couplings with no clock. **Why it's wrong:** $\mathbf f(t,\mathbf y)$ is *allowed* to depend on $t$ directly; that's the whole reason $t$ appears in the definition. **Fix:** treat time as an ordinary input you plug in — at $\mathbf k_2,\mathbf k_3$ use $t_n+\tfrac h2$, at $\mathbf k_4$ use $t_n+h$. --- ## L2 — Application ### Exercise 2.1 For $\mathbf f(t,\mathbf y)=(y_2,\,-y_1)$ with $\mathbf y_0=(1,0)$ and $h=0.2$, compute **only** $\mathbf k_1$ and $\mathbf k_2$. (Use the figure above: $\mathbf k_1$ is the start-of-step dot, $\mathbf k_2$ the first midpoint dot.) > [!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)$. > **Look at the figure's blue dot:** it is at $t_0+h/2$, and to build it both components got advanced *together* — the first slot used $-0.1$, not the stale $0$. ### Exercise 2.2 Take **one full RK4 step** of the system in 2.1 ($h=0.2$). Give $\mathbf y_1$ to 5 decimals, then compare to the true solution $(\cos t,-\sin t)$ at $t=0.2$. > [!recall]- Solution > $\mathbf k_1=(0,-1)$. > $\mathbf k_2=(-0.1,-1)$ (from 2.1). > State for $\mathbf k_3$: $(1,0)+0.1(-0.1,-1)=(0.99,\,-0.1)$, so $\mathbf k_3=(-0.1,\,-0.99)$. > State for $\mathbf k_4$: $(1,0)+0.2(-0.1,-0.99)=(0.98,\,-0.198)$, so $\mathbf k_4=(-0.198,\,-0.98)$. > Combine with weights $1,2,2,1$: > $$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$. Matches to ~5 and ~5 decimals — one step! ### Exercise 2.3 For $y'=y$ ($f=y$, scalar), take one RK4 step from $y_0=1$ with $h=0.5$. Compare to $e^{0.5}$. > [!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, tiny error even at $h=0.5$. > [!mistake] L2 trap — advancing one component through all four stages before touching the other. > **Why it feels right:** each equation "looks" separate, so you want to finish $y_1$'s whole RK4 first. **Why it's wrong:** $\mathbf k_2$ of component 1 needs the *advanced* value of $y_2$ — if $y_2$ is still stale, you feed old data and lose accuracy (or worse, a wrong answer). **Fix:** compute the **entire vector** $\mathbf k_1$, then form the full intermediate state, then the entire $\mathbf k_2$. Stages are synchronized across components. --- ## L3 — Analysis ### Exercise 3.1 For the linear system $\mathbf y'=A\mathbf y$ with $A=\begin{pmatrix}0&1\\-1&0\end{pmatrix}$ and $\mathbf y_0=(1,0)$, the true energy $E=y_1^2+y_2^2$ is exactly conserved ($E=1$ for all $t$). After **one** RK4 step with $h=0.2$ (Exercise 2.2 gave $\mathbf y_1=(0.980067,-0.198667)$), compute $E_1=y_{1,1}^2+y_{2,1}^2$. Is energy exactly conserved by RK4? Explain. > [!recall]- Solution > Using the exact rational one-step result, $E_1\approx 0.99999911$ — i.e. $E_1$ differs from $1$ by about $9\times10^{-7}$. > Energy is **almost** conserved but not exactly: RK4's step map is a polynomial approximation of the true rotation, so it slightly changes $E$. The drift is tiny here ($\sim10^{-6}$) because the local error is $O(h^5)$. Over *many* steps this drift accumulates. See [[Stiff systems and stability]] for when such structure-preservation matters. ### Exercise 3.2 The scalar test problem $y'=\lambda y$ with a single RK4 step multiplies $y$ by a fixed factor $R(h\lambda)$. Derive $R(z)$ where $z=h\lambda$. > [!recall]- Solution > With $f=\lambda y$: $k_1=\lambda y$. Then $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)$. > Now $y_1=y+\tfrac h6(k_1+2k_2+2k_3+k_4)$. Collecting powers of $z$: > $$R(z)=1+z+\frac{z^2}{2}+\frac{z^3}{6}+\frac{z^4}{24}.$$ > **What it means:** this is exactly the first 5 terms of $e^z$ — RK4 reproduces the true growth $e^{h\lambda}$ up to and including $z^4$, and *that* is where "order 4" comes from. See [[Local vs Global Truncation Error]]. ### Exercise 3.3 Using $R(z)$ from 3.2 with $z=0.5$, confirm it reproduces the answer of Exercise 2.3. > [!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$. > Identical to $y_1=1.648438$ from 2.3. **Why:** a single scalar RK4 step *is* multiplication by $R(z)$, so the two computations must agree exactly. > [!mistake] L3 trap — "order 4 means the answer is exact." > **Why it feels right:** matching four Taylor terms feels like matching everything. **Why it's wrong:** $R(z)$ stops at $z^4$; the true $e^z$ has $\tfrac{z^5}{120}+\cdots$ which RK4 *misses*. That missing tail is the $O(h^5)$ local error. **Fix:** treat RK4 as "correct to 4 terms of the exponential" — excellent but not exact. --- ## L4 — Synthesis ### Exercise 4.1 Reduce $y''+0.4y'+y=0$, $y(0)=2$, $y'(0)=0$ to a first-order system, then take **one** RK4 step with $h=0.1$. Report $\mathbf y_1=(u_{1,1},u_{2,1})$ to 5 decimals. > [!recall]- Solution > **Reduce.** $u_1=y,\ u_2=y'$. From $y''=-0.4y'-y=-0.4u_2-u_1$: > $$\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.$$ > Bracket sum for $u_2$: $-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.$$ > So $\mathbf y_1\approx(1.99014,\,-0.19573)$. The negative $u_2$ (velocity) shows the mass starting to swing back; the $0.4u_2$ term is friction — this is a **damped** oscillator, unlike the frictionless Example 1. ### Exercise 4.2 For the coupled system $y_1'=y_1-y_1 y_2,\ \ y_2'=-y_2+y_1 y_2$ (a Lotka–Volterra style model) with $\mathbf y_0=(1,1)$, compute $\mathbf k_1$ and explain what it says about the motion at $t=0$. > [!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)$ is an **equilibrium** — both slopes are zero, so to first order nothing moves. RK4 will predict $\mathbf y_1$ very close to $(1,1)$; the tiny drift comes only from the nonlinearity felt at the midpoint stages. This is a **degenerate/fixed-point** case you must recognize: a zero $\mathbf k_1$ does not mean you're done, but it does mean the start-slope contributes nothing. > [!mistake] L4 trap — freezing the *time* between stages when reducing a driven ODE. > **Why it feels right:** in Example 4.1 the system had no explicit $t$, so beginners assume $t$ never matters. **Why it's wrong:** if the ODE is $y''+y=\cos t$ (a forced oscillator), $\mathbf f$ depends on $t$; then $\mathbf k_2,\mathbf k_3$ must use $\cos(t_n+\tfrac h2)$ and $\mathbf k_4$ must use $\cos(t_n+h)$. Using $\cos t_n$ everywhere silently drops accuracy. **Fix:** always advance BOTH the state *and* the time slot at each stage. --- ## L5 — Mastery ### Exercise 5.1 A projectile with linear air drag obeys $\mathbf r''=-g\,\hat{\mathbf j}-c\,\mathbf r'$ in 2D. Write the **full first-order system** in the four unknowns $(x,y,v_x,v_y)$, with $g=9.8$, $c=0.1$. Then, for $\mathbf y_0=(0,0,10,10)$ and $h=0.1$, compute $\mathbf k_1$ and $\mathbf k_2$ (all four components each). > [!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.$$ > So $\mathbf f=(v_x,\ v_y,\ -0.1v_x,\ -9.8-0.1v_y)$, an $n=4$ vector. > **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).$$ > **Read it:** components 1–2 are the *velocities* (they feed positions), components 3–4 are the *accelerations*. Every one of the four moved together before $\mathbf k_2$ — the synchronization rule at full scale. ### Exercise 5.2 Complete one full RK4 step of Exercise 5.1 and report $\mathbf y_1=(x_1,y_1,v_{x,1},v_{y,1})$ to 5 decimals. > [!recall]- Solution > Continue from 5.1. > **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** with $\tfrac{0.1}{6}$ and weights $1,2,2,1$. For $x_1$ the bracket sum is > $10+2(9.95)+2(9.95025)+9.9004975=59.7009975$, so > $$x_1=0+\tfrac{0.1}{6}(59.7009975)=0.995017.$$ > For $y_1$ the bracket sum is $10+2(9.46)+2(9.4627)+8.925373=56.770873$, so > $$y_1=0+\tfrac{0.1}{6}(56.770873)=0.946181.$$ > For the velocities: > $$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.$$ > So $\mathbf y_1\approx(0.99502,\ 0.94618,\ 9.90050,\ 8.92536)$. Gravity has already pulled $v_y$ down from $10$ toward $8.93$, and drag nudged $v_x$ below $10$ — physically correct. ### Exercise 5.3 Explain, using $R(z)=1+z+\tfrac{z^2}{2}+\tfrac{z^3}{6}+\tfrac{z^4}{24}$ from 3.2, why RK4 with a large step $h$ on the drag component $v_x'=-0.1v_x$ stays stable, but why for a **stiff** eigenvalue like $\lambda=-1000$ you would need $|h\lambda|$ small. > [!recall]- Solution > Stability of one component requires $|R(h\lambda)|\le 1$. For the drag component $\lambda=-0.1$, even $h=1$ gives $z=-0.1$ and $R(-0.1)\approx0.90484$, safely inside $[-1,1]$ — stable. > For $\lambda=-1000$, taking any moderate $h$ (say $h=0.01$) gives $z=-10$; then $R(-10)=1-10+50-\tfrac{1000}{6}+\tfrac{10000}{24}=1-10+50-166.6\overline{6}+416.6\overline{6}=291$, whose magnitude is $\gg 1$ — the numerical solution **explodes** even though the true solution decays. To stay stable you must shrink $h$ so $|z|$ is small (roughly $|z|\lesssim 2.78$ for real negative $z$). That forced tiny step is the definition of **stiffness**. See [[Stiff systems and stability]]. > [!mistake] L5 trap — "smaller error always means smaller step is safe; big $h$ only costs accuracy." > **Why it feels right:** for smooth non-stiff problems, larger $h$ merely raises the $O(h^4)$ error. **Why it's wrong:** for a stiff eigenvalue, once $|R(h\lambda)|>1$ the error doesn't just grow — it **blows up geometrically**, turning a decaying true solution into numerical garbage. **Fix:** for stiff systems check the *stability* condition $|R(h\lambda)|\le1$, not just the accuracy order; or switch to an implicit method. --- ## Connections - [[4.8.27 Systems of ODEs — RK4 for systems (Hinglish)|Parent topic — RK4 for systems]] - [[RK4 for a single ODE]] — every exercise here reduces to this on each component. - [[Reducing higher-order ODEs to first-order systems]] — used in 1.3, 4.1, 5.1. - [[Euler's method for systems]] — compare the accuracy you'd lose. - [[Local vs Global Truncation Error]] — the $R(z)$ story in L3. - [[Stiff systems and stability]] — L5.3. - [[Simpson's Rule]] — the quadrature the weights $1,2,2,1$ imitate.