# Administrative Information *Plan for today: "**Euler's method for systems.**" (1) **Recall** scalar forward Euler $y_{i+1}=y_i+h\,f(t_i,y_i)$, now pulled into its own function `myEuler`. (2) **Go vector-valued** — most interesting dynamics are systems or higher-order equations. Rewrite a second-order ODE as a first-order system: $y''=-y$ becomes $\vec{y}=\begin{bmatrix}y_1\\y_2\end{bmatrix},\qquad \vec{y}\,'=f(t,\vec{y})=\begin{bmatrix}y_2\\-y_1\end{bmatrix},\qquad \vec{y}(0)=\begin{bmatrix}1\\0\end{bmatrix}.$ (3) **The key point** — the Euler update does not change at all: $\vec{y}_{i+1}=\vec{y}_i+h\,f(t_i,\vec{y}_i)$, only now the state $\vec{y}$ is a vector and $f$ returns a vector. (4) **Code it** (`LiveIVP02`) and compare the two components against the exact $(\cos t,-\sin t)$; watch the phase-plane orbit, a circle that Euler slowly **spirals outward** because a first-order method does not conserve the energy.* > [!info]- On the horizon (announced) > - **Homework 3** is posted, **due July 31** (adaptive quadrature and numerical ODEs — you will use your adaptive Simpson's method plus RK4/RKF and a built-in solver). > - **Next:** the Runge-Kutta midpoint method (a second-order stepper), then higher-order RK. > - 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 keep `myEuler` from last time and add **`LiveIVP02`** (`codes/16/day16_live/`): forward Euler on the **system** $y_1'=y_2,\ y_2'=-y_1$ (the harmonic oscillator $y''=-y$), with initial state $(1,0)$. The state `y` becomes a two-column array (each row is $(y_1,y_2)$ at that time), `f` returns a two-vector, and the **update line is unchanged**. Plotted against $\cos t$ and $-\sin t$ the components match; the phase-plane orbit $(y_1,y_2)$ is very nearly the unit circle but drifts outward (radius $1.00\to1.05$ over $[0,10]$) — first-order Euler leaks energy. Swapping the vector field to $f=[\,y_2,\,-\sin y_1\,]$ turns the very same driver into a nonlinear **pendulum** (see the finale below). Four languages below; Python is run-verified. > [!example]- Python - LiveIVP02 (systems Euler, run-verified) > ```python > import numpy as np > f = lambda t, y: np.array([y[1], -y[0]]) # vector field: y1'=y2, y2'=-y1 > 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]) # SAME Euler step, now on a vector > # y[:,0] ~ cos t, y[:,1] ~ -sin t; phase plane (y[:,0], y[:,1]) drifts outward > ``` > [!example]- MATLAB - LiveIVP02 (systems Euler) > ```matlab > f = @(t,y) [y(2), -y(1)]; % vector field: 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,:)); % SAME Euler step, now on a vector > 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: a circle Euler spirals out of > ``` > [!example]- R - LiveIVP02 (systems Euler) > ```r > f <- function(t, y) c(y[2], -y[1]) # vector field > 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, ]) # Euler step > # matplot(t, cbind(y[,1], cos(t), y[,2], -sin(t)), type="l"); plot(y[,1], y[,2]) > ``` > [!example]- Mathematica - LiveIVP02 (systems Euler) > ```wolfram > f[t_, y_] := {y[[2]], -y[[1]]}; (* vector field *) > ic = {1., 0.}; t0 = 0.; tF = 10.; h = 0.01; > t = Range[t0, tF, h]; > y = ConstantArray[0., {Length[t], 2}]; y[[1]] = ic; (* each row is the state *) > Do[y[[i + 1]] = y[[i]] + h f[t[[i]], y[[i]]], {i, 1, Length[t] - 1}]; (* Euler step *) > ``` # Lecture Boards + Transcript + GenAI ## Recall: scalar Euler as a function For a scalar initial value problem $y'=f(t,y)$ (the **rule of evolution**) with $y(t_0)=y_0\in\mathbb{R}$, forward Euler steps the solution forward one $h$ at a time, $y_{n+1}=y_n+h\,f(t_n,y_n),\qquad n=0,1,2,\dots,$ with global accuracy $O(h)$. That single line, now living in `myEuler`, is the entire method. The point of today is that it needs almost no change to handle far richer problems. > [!example]- Board — recall the scalar IVP and Euler, $O(h)$ > ![[day16_board1.png]] ## Second-order ODEs as first-order systems Interesting dynamics rarely come as a single scalar equation. They arrive as **systems**, and in two common ways. **Directly, as coupled equations.** Population models are a clean example. The linear law $x'=kx$ alone gives pure exponential growth $x=Ce^{kt}$ — the function whose derivative is a multiple of itself — fine at first but unphysical forever (Malthus's worry: an exponentially growing population outrunning food that grows only linearly). The fix is a **nonlinear correction**, the **logistic** law $x'=kx(1-x/N)$, whose $(1-x/N)$ factor bends growth down to a stable carrying capacity $N$. Couple two logistic species with interaction terms and you get a **competing-species** model, $\frac{dx}{dt}=2x\Big(1-\frac{x}{2}\Big)-xy,\qquad \frac{dy}{dt}=3y\Big(1-\frac{y}{3}\Big)-2xy.$ The $2x(1-x/2)$ and $3y(1-y/3)$ pieces are each species' self-limiting growth; the $-xy$ and $-2xy$ terms are the competition (each eats at the other). This is a **first-order, two-dimensional, nonlinear** system — and the honest reason numerical methods matter here is that the symbolic tools you would reach for (equilibria, linearization, phase lines) get hard fast once the equations are nonlinear. **Or from many interacting objects.** Take a small box of gas and attach to each molecule a position vector and a velocity vector; Newton's laws then give one enormous coupled system. You will not simulate $10^{23}$ particles, but the *structure* is identical — many states marched by a single rule. **By reduction of a higher-order linear ODE.** The damped, driven mass-spring collects three forces — inertia $my''$, damping $by'$ (proportional to velocity, bleeding off kinetic energy until the mass turns around), and the Hookean spring $ky$ (restoring, proportional to displacement): $m y''+b y'+k y=f(t),\qquad y(t_0)=y_0,\ \ y'(t_0)=y_0',$ It becomes first order by naming the velocity $v=y'$ as a second state. Then $v'=y''=-\tfrac{k}{m}y-\tfrac{b}{m}v+\tfrac{f(t)}{m}$, so $y'=F(t,y,v)=v,\qquad v'=g(t,y,v)=-\frac{k}{m}\,y-\frac{b}{m}\,v+\frac{f(t)}{m}.$ The state now **carries two things, a displacement and a velocity**. Any $n$th-order ODE collapses this way into $n$ first-order equations by carrying the lower derivatives as extra state variables. **Key point.** A good stepper should **not care** whether the system is linear or nonlinear, nor how many dynamic quantities are evolving. All it needs is a state and a rule for its rate of change. > [!example]- Board — competing species (nonlinear system) and higher-order reduction $v=y'$ > ![[day16_board2.png]] ## Euler for systems: the update is unchanged Stack the state into a vector $\vec{y}$ and let $f$ return the vector of rates. Then forward Euler is *character-for-character the same*: $\vec{y}_{\,n+1}=\vec{y}_{\,n}+h\,f(t_n,\vec{y}_{\,n}).$ Take the harmonic oscillator $y''=-y$ (the mass-spring with $m=k=1,\ b=0,\ f=0$). Writing $y_1=y,\ y_2=y'$, $\vec{y}\,'=f(t,\vec{y})=\begin{bmatrix}y_2\\-y_1\end{bmatrix},\qquad \vec{y}(0)=\begin{bmatrix}1\\0\end{bmatrix}\ \Rightarrow\ \vec{y}(t)=\begin{bmatrix}\cos t\\-\sin t\end{bmatrix}.$ In `LiveIVP02` the state `y` is a two-column array (one row per time), `f` returns a two-vector, and the update line never changes. The two components track $\cos t$ and $-\sin t$; the exact orbit in the phase plane $(y_1,y_2)$ is the unit circle. Euler's orbit, however, **spirals slowly outward** — the radius grows from $1.00$ to about $1.05$ over $[0,10]$ — because a first-order method does not conserve the oscillator's energy. Shrinking $h$ slows the drift (it is $O(h)$) but never removes it; that is exactly the motivation for the higher-order steppers coming next. ![[day16_systems_euler.png]] ## From linear to nonlinear: the pendulum, and why swings are fun The same code solves a genuinely **nonlinear** system with a one-line change of the vector field. A pendulum obeys $y''=-\sin y$ (with $y$ the angle), i.e. the system $y_1'=y_2,\ y_2'=-\sin y_1$. For **small** swings the small-angle approximation $\sin y\approx y$ turns this back into the mass-spring $y''=-y$, whose phase-plane orbit is a circle. Start from a **larger** displacement, though, and the approximation breaks: the orbits stop being circles and their sides **pinch inward** toward the turning points, becoming lens (or "eye") shapes as the amplitude approaches a full flip. That pinch is something you have felt. Near the top of a big swing the pendulum slows — its velocity bleeds to zero at the turning point and lingers before returning — which is the slow, hang-time part of a playground swing (live, this was the vector-field demo where the styrofoam peanuts crawl through the pinched region and race through the bottom). The linear mass-spring gives you only the boring circle; the nonlinearity is what makes swings fun. It is also the whole point of the method: nothing in the Euler loop changed, yet it handled a nonlinear system the symbolic toolkit would struggle with. ![[day16_pendulum_phase.png]] # Check Your Understanding (CYU) 1. **Reduce to a system.** Write $y''+0.1\,y'+4y=\cos(2t)$, $y(0)=1,\ y'(0)=0$, as a first-order system in $(y,v)$ with $v=y'$. What are $y'$ and $v'$? 2. **Read the model.** In $x'=2x(1-x/2)-xy$, what behavior does the $2x(1-x/2)$ term encode on its own, and what does the $-xy$ term add? 3. **Why the drift?** Euler's phase-plane orbit for $y''=-y$ spirals outward. Does halving $h$ eliminate the outward drift or only slow it, and why? 4. **Shapes in code.** If the state $\vec{y}$ has two components, what is the shape of the value returned by `f(t, y)`, and of the full solution array over all time steps?