# Administrative Information *Plan for today: "**Runge-Kutta methods: RK2, then RK4.**" (1) **Recall** Euler's weakness — it commits to a single slope at the *left* endpoint, so it is only first order, $O(h)$ — and last class's midpoint idea. (2) **RK2** — sample a *second* slope at a trial point inside the step and combine the two; this is the midpoint/Heun family, **second order** $O(h^2)$. Derive it and implement it live. (3) **RK4** — take *four* slope samples per step, $k_1,k_2,k_3,k_4$, and step with their weighted average $(1,2,2,1)/6$; **fourth order** $O(h^4)$, the classic workhorse. Implement `myRK4` and confirm the order jumps to $\approx 4$.* > [!info]- On the horizon (announced) > - **Homework 3** (adaptive quadrature + numerical ODEs) is due July 31 — RK2/RK4 are the engines for its ODE problems. > - Reflection post due tonight. ## Current Code State [Dropbox Share of Live Codes](https://www.dropbox.com/scl/fo/3ir4tbnisnzstqia3hbep/ANyMKUyY2OnTkVgqdxYlnbE?rlkey=xonhgv5o5k90hprp6zwijn9w8&dl=0) We pick up from last class's two live scripts (`codes/17/day17_live/`): **`LiveIVP01`**, scalar forward Euler on $y'=y,\ y(0)=1$ (exact $e^{t}$), and **`LiveIVP02`**, 2D forward Euler on the oscillator $y''=-y$ written as $y_1'=y_2,\ y_2'=-y_1$ (exact $(\cos t,-\sin t)$). Both are the plain inline stepping loops. The Runge-Kutta **RK4 staging** we will build toward is at the end of this section. Four languages; Python is run-verified. > [!example]- Python - LiveIVP01 (scalar Euler, starting point) > ```python > import numpy as np > f = lambda t, y: y > ic, tInitial, tFinal, h = 1.0, 0.0, 10.0, 0.01 > t = np.arange(tInitial, tFinal + h, h) > y = np.zeros(len(t)); y[0] = ic > for i in range(len(t) - 1): > y[i+1] = y[i] + h*f(t[i], y[i]) # forward Euler step > # compare y to np.exp(t) > ``` > [!example]- MATLAB - LiveIVP01 (scalar Euler, starting point) > ```matlab > f = @(t,y) y; > ic = 1; tInitial = 0; tFinal = 10; h = 0.01; > t = tInitial:h:tFinal; > y = zeros(length(t),1); y(1) = ic; > for i = 1:length(t)-1 > y(i+1) = y(i) + h*f(t(i), y(i)); % forward Euler step > end > hold on; plot(t, y, 'r'); plot(t, exp(t), 'b'); hold off > ``` > [!example]- R - LiveIVP01 (scalar Euler, starting point) > ```r > f <- function(t, y) y > ic <- 1; tInitial <- 0; tFinal <- 10; h <- 0.01 > t <- seq(tInitial, tFinal, by = h) > y <- numeric(length(t)); y[1] <- ic > for (i in 1:(length(t)-1)) y[i+1] <- y[i] + h*f(t[i], y[i]) # forward Euler > ``` > [!example]- Mathematica - LiveIVP01 (scalar Euler, starting point) > ```wolfram > f[t_, y_] := y; > ic = 1.; tInitial = 0.; tFinal = 10.; h = 0.01; > t = Range[tInitial, tFinal, h]; > y = ConstantArray[0., Length[t]]; y[[1]] = ic; > Do[y[[i+1]] = y[[i]] + h f[t[[i]], y[[i]]], {i, 1, Length[t]-1}]; (* forward Euler *) > ``` > [!example]- Python - LiveIVP02 (2D Euler, starting point) > ```python > import numpy as np > f = lambda t, y: np.array([y[1], -y[0]]) # y1'=y2, y2'=-y1 (i.e. y''=-y) > ic = np.array([1.0, 0.0]) # y1(0)=1, y2(0)=0 > t0, tF, h = 0.0, 10.0, 0.01 > t = np.arange(t0, tF + h, h) > y = np.zeros((len(t), 2)); y[0] = ic # each row is the state (y1, y2) > for i in range(len(t) - 1): > y[i+1] = y[i] + h*f(t[i], y[i]) # forward Euler step, vector state > # y[:,0] ~ cos t, y[:,1] ~ -sin t; phase plane drifts outward > ``` > [!example]- MATLAB - LiveIVP02 (2D Euler, starting point) > ```matlab > f = @(t,y) [y(2), -y(1)]; % y1'=y2, y2'=-y1 > ic = [1, 0]; % y1(0)=1, y2(0)=0 > tInitial = 0; tFinal = 10; h = 0.01; > t = tInitial:h:tFinal; > y = zeros(length(t), 2); % each ROW is the state at that time > y(1,:) = ic; > for i = 1:length(t)-1 > y(i+1,:) = y(i,:) + h*f(t(i), y(i,:)); % forward Euler step, vector state > end > figure; hold on > plot(t, y(:,1)); plot(t, cos(t)); plot(t, y(:,2)); plot(t, -sin(t)); > legend('~cos','cos','~sin','sin'); hold off > figure; plot(y(:,1), y(:,2)) % phase plane > ``` > [!example]- R - LiveIVP02 (2D Euler, starting point) > ```r > f <- function(t, y) c(y[2], -y[1]) # y1'=y2, y2'=-y1 > ic <- c(1, 0); t0 <- 0; tF <- 10; h <- 0.01 > t <- seq(t0, tF, by = h) > y <- matrix(0, length(t), 2); y[1, ] <- ic # each row is the state > for (i in 1:(length(t)-1)) y[i+1, ] <- y[i, ] + h*f(t[i], y[i, ]) # forward Euler > ``` > [!example]- Mathematica - LiveIVP02 (2D Euler, starting point) > ```wolfram > f[t_, y_] := {y[[2]], -y[[1]]}; (* y1'=y2, y2'=-y1 *) > ic = {1., 0.}; t0 = 0.; tF = 10.; h = 0.01; > t = Range[t0, tF, h]; > y = ConstantArray[0., {Length[t], 2}]; y[[1]] = ic; > Do[y[[i+1]] = y[[i]] + h f[t[[i]], y[[i]]], {i, 1, Length[t]-1}]; (* forward Euler *) > ``` **The slope picture.** RK2 probes one extra slope at a midpoint and steps with it; RK4 samples four slopes across the step ($k_1$ at the start, $k_2$ and $k_3$ at the middle, $k_4$ at the end) and steps with their weighted average. ![[day17_rk2_slopes.png]] ![[day17_rk4_slopes.png]] > [!quote]- External visualizations (attributed) > Two outside resources give nice slope-by-slope picture sequences; linked here rather than reproduced, per their licenses. > - Harold Serrano, *Visualizing the Runge-Kutta Method* — a step-by-step picture sequence of the four RK4 slopes. https://www.haroldserrano.com/blog/visualizing-the-runge-kutta-method > - Autar K. Kaw, *Runge-Kutta 2nd Order Method for Solving ODEs*, in **Numerical Methods with Applications** (Holistic Numerical Methods, mathforcollege.com), Ch. 08.03 — including the RK2 slope figure. Licensed CC BY-NC-ND 4.0. https://nm.mathforcollege.com/NumericalMethodsTextbookUnabridged/chapter-08.03-runge-kutta-2nd-order-method-for-solving-ordinary-differential-equations.html (the specific slope graphic: `08.03.graph16.png`). ### The RK4 staging (where we are headed) Each step samples four slopes and steps with their weighted average: $k_1=h\,f(t_i,y_i),\quad k_2=h\,f\!\big(t_i+\tfrac{h}{2},\,y_i+\tfrac12 k_1\big),\quad k_3=h\,f\!\big(t_i+\tfrac{h}{2},\,y_i+\tfrac12 k_2\big),\quad k_4=h\,f\!\big(t_i+h,\,y_i+k_3\big),$ $y_{i+1}=y_i+\tfrac16\big(k_1+2k_2+2k_3+k_4\big).$ On $y'=y$ over $[0,2]$, RK4's error at $t=2$ falls by about $16\times$ each time $h$ is halved (order $\approx4$): $1.7\times10^{-4}\to1.1\times10^{-5}\to7.4\times10^{-7}\to4.7\times10^{-8}$. > [!example]- Python - myRK4 (staging, run-verified) > ```python > import numpy as np > def myRK4(tInitial, tFinal, h, f, ic): > t = np.arange(tInitial, tFinal + h, h) > y = np.zeros(len(t)); y[0] = ic > for i in range(len(t) - 1): > k1 = h*f(t[i], y[i]) > k2 = h*f(t[i]+h/2, y[i]+0.5*k1) > k3 = h*f(t[i]+h/2, y[i]+0.5*k2) > k4 = h*f(t[i]+h, y[i]+k3) > y[i+1] = y[i] + (k1 + 2*k2 + 2*k3 + k4)/6 > return y, t > ``` > [!example]- MATLAB - myRK4 (staging) > ```matlab > function [y, t] = myRK4(tInitial, tFinal, h, f, ic) > t = tInitial:h:tFinal; > y = zeros(length(t), 1); y(1) = ic; > for i = 1:length(t)-1 > k1 = h*f(t(i), y(i)); > k2 = h*f(t(i)+h/2, y(i)+1/2*k1); > k3 = h*f(t(i)+h/2, y(i)+1/2*k2); > k4 = h*f(t(i)+h, y(i)+k3); > y(i+1) = y(i) + 1/6*(k1 + 2*k2 + 2*k3 + k4); > end > end > ``` > [!example]- R - myRK4 (staging) > ```r > myRK4 <- function(tInitial, tFinal, h, f, ic) { > t <- seq(tInitial, tFinal, by = h) > y <- numeric(length(t)); y[1] <- ic > for (i in 1:(length(t)-1)) { > k1 <- h*f(t[i], y[i]) > k2 <- h*f(t[i]+h/2, y[i]+0.5*k1) > k3 <- h*f(t[i]+h/2, y[i]+0.5*k2) > k4 <- h*f(t[i]+h, y[i]+k3) > y[i+1] <- y[i] + (k1 + 2*k2 + 2*k3 + k4)/6 > } > list(y = y, t = t) > } > ``` > [!example]- Mathematica - myRK4 (staging) > ```wolfram > myRK4[tInitial_, tFinal_, h_, f_, ic_] := Module[{t, y, k1, k2, k3, k4}, > t = Range[tInitial, tFinal, h]; > y = ConstantArray[0., Length[t]]; y[[1]] = ic; > Do[ > k1 = h f[t[[i]], y[[i]]]; > k2 = h f[t[[i]]+h/2, y[[i]]+1/2 k1]; > k3 = h f[t[[i]]+h/2, y[[i]]+1/2 k2]; > k4 = h f[t[[i]]+h, y[[i]]+k3]; > y[[i+1]] = y[[i]] + 1/6 (k1 + 2 k2 + 2 k3 + k4), > {i, 1, Length[t]-1}]; > {y, t}]; > ``` # Lecture Boards + Transcript + GenAI > [!note] Transcripts pending > Today's audio was not yet transcribed at build time; the narrative below is reconstructed from the boards and the live code, and can be enriched once the transcript lands. ## Recall: Euler, and the value of a better slope Euler comes from truncating a Taylor step: $y(t_i+h)=y(t_i)+y'(t_i)h+O(h^2)$, and since the ODE hands us $y'(t_i)=f(t_i,y_i)$, dropping the remainder gives $y_{i+1}\approx y_i+h\,f(t_i,y_i)$. That single step is *quadratically* accurate **locally**, but summed over the $\sim 1/h$ steps needed to cross a fixed interval it is only *linearly* accurate **globally** ($O(h)$). Today's driving question: **can we get a more accurate time-update?** > [!example]- Board — recall Euler from Taylor; locally $O(h^2)$, globally $O(h)$ > ![[day17_board1.png]] ## RK2: a second slope, second order **Retain more of the Taylor series.** Keep the quadratic term this time: $y_{i+1}=y_i+y'(t_i)\,h+\frac{y''(t_i)}{2}h^2+O(h^3).$ We already know $y'=f(t_i,y_i)$. For $y''$, differentiate $f(t,y(t))$ by the chain rule, using $\tfrac{dy}{dt}=f$: $y''=\frac{d}{dt}f(t,y)=\frac{\partial f}{\partial t}+\frac{\partial f}{\partial y}\frac{dy}{dt}=f_t+f_y\,f.$ Substituting, $y_{i+1}=y_i+h\Big[\,f(t_i,y_i)+\tfrac{h}{2}f_t(t_i,y_i)+\tfrac{h}{2}f_y(t_i,y_i)\,f(t_i,y_i)\,\Big]+O(h^3).$ **Recognize the bracket.** The 2D Taylor expansion of $f$ about $(t_i,y_i)$ is $f(t,y)=f(t_i,y_i)+f_t\,(t-t_i)+f_y\,(y-y_i)+\dots$. Choose $t-t_i=\tfrac{h}{2}$ and $y-y_i=\tfrac{h}{2}f(t_i,y_i)$: then $f\big(t_i+\tfrac h2,\ y_i+\tfrac h2 f(t_i,y_i)\big)=f(t_i,y_i)+\tfrac h2 f_t+\tfrac h2 f_y f+\dots$ — *exactly* the bracket. Therefore $\boxed{\,y_{i+1}=y_i+h\,f\!\Big(t_i+\tfrac{h}{2},\ y_i+\tfrac{h}{2}\,f(t_i,y_i)\Big)\,}$ the **midpoint** method (RK2): take a half Euler step to the middle, **resample the slope there**, and step the full $h$ with that better slope. Because it matches the Taylor series through $O(h^3)$, it is **locally $O(h^3)$, globally $O(h^2)$** — one order better than Euler. Verified in the sandbox: recovered order $\approx 2$ (doubling the resolution quarters the error). > [!example]- Boards — retain the $h^2$ term, $y''=f_t+f_yf$, and recognize the midpoint resample > ![[day17_board2.png]] > ![[day17_board3.png]] > ![[day17_board4.png]] > [!example]- Python - myRK2 (midpoint, run-verified) > ```python > import numpy as np > def myRK2(tInitial, tFinal, h, f, ic): > t = np.arange(tInitial, tFinal + h, h) > y = np.zeros((len(t),) + np.shape(ic)); y[0] = ic > for i in range(len(t) - 1): > y[i+1] = y[i] + h*f(t[i]+h/2, y[i] + (h/2)*f(t[i], y[i])) # midpoint step > return y, t > ``` > [!example]- MATLAB - myRK2 (midpoint) > ```matlab > function [y, t] = myRK2(tInitial, tFinal, h, f, ic) > t = tInitial:h:tFinal; > y = zeros(length(t), 1); y(1) = ic; > for i = 1:length(t)-1 > y(i+1) = y(i) + h*f(t(i)+h/2, y(i) + (h/2)*f(t(i), y(i))); % midpoint step > end > end > ``` ## RK4: four slopes, fourth order RK4 keeps going: sample **four** slopes per step ($k_1$ at the start, $k_2,k_3$ at the middle, $k_4$ at the end) and step with the weighted average $(1,2,2,1)/6$ — the staging in the Current Code State above. Matching the Taylor series through $O(h^5)$ makes it **globally $O(h^4)$**. Everything generalizes verbatim to **systems** (`myRK2Sys`, `myRK4Sys`): the state $\vec y$ is a vector, $f$ returns a vector, and the update lines are character-for-character the same. **The payoff.** `LiveIVP04` races Euler vs RK2 vs RK4 on the oscillator $y''=-y$ (state $(y,y')$, start $(3,0)$) with a *deliberately coarse* $h=0.25$ out to $t=20$. In the phase plane the three separate cleanly: Euler pumps in so much energy the orbit spirals from radius $3$ out to $\approx 34$; RK2 nearly closes ($3\to3.12$); RK4 sits on the circle ($3\to3.000$). Recovered orders on $y'=y$ confirm the story: Euler $\approx1$, RK2 $\approx2$, RK4 $\approx4$. ![[day17_methods_phase.png]] # Check Your Understanding (CYU) 1. **The chain rule step.** In the RK2 derivation, why is $y''=f_t+f_y f$ (not just $f_t$)? Which rule brings in the $f_y f$ term? 2. **One RK2 step by hand.** For $y'=y$, $y(0)=1$, $h=0.5$: compute $y_1$ with the midpoint method and compare to Euler's $y_1=1.5$ and the exact $e^{0.5}=1.6487$. 3. **Order counting.** RK2 matches Taylor through $O(h^3)$ locally. Why is its *global* order $O(h^2)$ rather than $O(h^3)$? 4. **Why RK4 stays on the circle.** In the $h=0.25$ comparison, Euler's orbit blew up to radius $\approx34$ while RK4 held $3.000$. What quantity is each method (failing to) conserve, and why does higher order help?