# HW3 Companion — Adaptive Quadrature & Numerical ODEs > [!info]- How to use this companion > This is a scaffold, not a solution key. For each problem you get: > - **What to hand in** — what is a *written* calculation, what is a *code* calculation, and what to *report*. > - **Hints** — a ladder of collapsed callouts. Open **Hint 1** first; open the next only if still stuck. > - **Code guidance** — a callout per language. You already built the engines in lecture (your adaptive > Simpson's, RK4, and RKF); the code here shows how to *drive* them for each problem, with `# TODO` > marks where the real setup goes. Problem 5 uses a built-in solver. > > Submit **one written PDF** (calculations, tables, plots, short reflections) plus your **code files** and **figures**. Label every axis. --- ## Deliverables at a glance | # | Topic | By hand | In code | Report | |---|-------|---------|---------|--------| | 1 | Adaptive Simpson's | the exact integral (antiderivative) | adaptive + uniform Simpson, count points | value + agreement, the two point-counts, node picture | | 2 | Time-varying logistic | your hypotheses first | RK4, both cases | trajectory plots, hypothesis vs. result | | 3 | Harvested logistic | (none required) | RK4 sweeps of $P(0)$ and $H$ | fates of trajectories, threshold vs. attractor, sudden vs. gradual | | 4 | Mass-spring (RKF) | the four regime conditions | RKF: 4 undriven + 2 driven | phase-plane portraits, resonance & beating | | 5 | Duffing | energy / double-well intuition | built-in ODE solver | $y(t)$ and phase plots, periodicity, chaos | --- ## Problem 1 — Adaptive Simpson's on a spike **What to hand in.** - **By hand:** the exact value $\int_{-1000}^{1000}\frac{dx}{1+x^2}=2\arctan(1000)$ (an antiderivative away). - **In code:** run your adaptive Simpson's to $10^{-6}$; count its function evaluations; then a *uniform* composite Simpson's, refined until it matches that accuracy, counting its points too. - **Report:** the value and agreement, the two point-counts and their ratio, and a picture of where the adaptive nodes landed. > [!hint]- Hint 1 — counting the points > The simplest honest count: wrap the integrand so every call ticks a counter, then hand the *wrapped* function to your method. In Python, a tiny class or a global works; in MATLAB a `persistent` counter in a nested function; in R an environment; in Mathematica a `Module` variable. The count you get is exactly "points used." > [!hint]- Hint 2 — the uniform comparison > Run your Day-14 composite Simpson's with $n$ panels, doubling $n$ until $|\,S_n - 2\arctan(1000)\,|\le10^{-6}$. A uniform rule evaluates $f$ at $n+1$ points. Compare that to the adaptive count — expect the uniform rule to need roughly an order of magnitude more here. > [!hint]- Hint 3 — seeing where the work went > If your wrapper also *stores* each $x$ it was called at, a histogram (or a rug plot) of those $x$-values shows the nodes piling up near $x=0$ and thinning out in the tails. That picture is the whole point of adaptivity. > [!hint]- Hint 4 — the explanation > Simpson's error on a panel scales like $h^5 f^{(4)}$. Out in the flat tails $f^{(4)}\approx0$, so huge panels already beat the tolerance; near the peak $f^{(4)}$ is large, so the method must chop $h$ down. A uniform grid is forced to use the *smallest* needed $h$ everywhere, wasting almost all of its points on tails that never needed them. > [!example]- Python > > ```python > import numpy as np > class Counted: > def __init__(self, f): self.f, self.n, self.xs = f, 0, [] > def __call__(self, x): self.n += 1; self.xs.append(x); return self.f(x) > > g = Counted(lambda x: 1/(1 + x**2)) > I = adaptiveSimpsons(-1000, 1000, g, 1e-6) # your Day-15 routine > print(I, 2*np.arctan(1000), g.n) # value, exact, points used > # uniform: your composite Simpson simp(f,a,b,n); double n until |I-exact|<1e-6, count n+1 > ``` > *Look for:* `g.n` far below the uniform `n+1`; `g.xs` clustered near 0. > [!example]- R > > ```r > make_counter <- function(f) { e <- environment(); n <- 0; xs <- c() > list(f = function(x){ n <<- n + 1; xs <<- c(xs, x); f(x) }, > count = function() n, nodes = function() xs) } > g <- make_counter(function(x) 1/(1 + x^2)) > I <- adaptiveSimpsons(-1000, 1000, g$f, 1e-6) > c(I, 2*atan(1000), g$count()) > ``` > *Look for:* `g$count()` well below the uniform point count. > [!example]- MATLAB > > ```matlab > function y = counted(x) % put in its own file, or use a nested fn > persistent n xs > if isempty(n), n = 0; xs = []; end > n = n + 1; xs(end+1) = x; assignin('base','evalCount',n); assignin('base','nodes',xs); > y = 1./(1 + x.^2); > end > % I = adaptiveSimpsons(-1000, 1000, @counted, 1e-6); disp([I, 2*atan(1000), evalCount]) > ``` > *Look for:* `evalCount` far below the uniform `n+1`. > [!example]- Mathematica *(numerical only)* > > ```wolfram > cnt = 0; nodes = {}; > g[x_?NumericQ] := (cnt++; AppendTo[nodes, x]; 1/(1 + x^2)); > I = adaptiveSimpsons[-1000., 1000., g, 10^-6]; > {I, 2 ArcTan[1000.], cnt} > ``` > *Look for:* `cnt` far below the uniform point count; `Histogram[nodes]` peaks at 0. --- ## Problem 2 — The time-varying logistic (RK4) **What to hand in.** - **By hand:** write your hypothesis for each case *before* you simulate — that is the point of the exercise. - **In code:** integrate both cases with your RK4 and plot $P(t)$ from a few initial populations. - **Report:** the trajectory plots next to your hypotheses, and whether the numerics agreed. > [!hint]- Hint 1 — non-autonomous means $t$ enters the slope > Here $f(t,P)$ depends on $t$ *explicitly* (through $k(t)$ or $N(t)$), so make sure your RK4 actually passes the time into every stage: $k_2,k_3$ use $t+\tfrac h2$ and $k_4$ uses $t+h$. A solver that ignores $t$ will silently solve the wrong problem. > [!hint]- Hint 2 — case one, $k=1+\sin t,\ N=1$ > The growth rate breathes between $0$ and $2$ but never goes negative, so $N=1$ is still the ceiling every trajectory rises to — just with a rippling approach rather than a smooth one. Predict "settles to 1 with wiggles," then check. > [!hint]- Hint 3 — case two, $k=1,\ N=2+\sin t$ > Now the *carrying capacity itself* oscillates between $1$ and $3$. The population chases a moving target and settles into a sustained oscillation around $2$ (roughly the interval $[1.2,\,2.6]$), not a fixed value. Ask: does it lag the capacity? > [!hint]- Hint 4 — a free correctness check > Set $k=N=1$ (constants) and you must recover the ordinary logistic curve rising to $1$. If that fails, your RHS or your RK4 wiring is off — fix it before trusting the time-varying runs. > [!example]- Python > > ```python > import numpy as np > k = lambda t: 1 + np.sin(t); N = lambda t: 1.0 # case 1 > rhs = lambda t, P: k(t) * P * (1 - P / N(t)) > # drive with YOUR rk4 over t in [0, 40], h ~ 0.01, looping several P(0) > 0 > ``` > *Look for:* every $P(0)\in(0,1)$ (and above) rises to $\approx1$ with ripples; for case 2 (`N = lambda t: 2 + np.sin(t)`) it oscillates around 2. > [!example]- R > > ```r > k <- function(t) 1 + sin(t); N <- function(t) 1 # case 1 > rhs <- function(t, P) k(t) * P * (1 - P / N(t)) > # feed to your rk4 on [0,40]; case 2 uses N <- function(t) 2 + sin(t) > ``` > [!example]- MATLAB > > ```matlab > k = @(t) 1 + sin(t); N = @(t) 1; % case 1 > rhs = @(t,P) k(t) .* P .* (1 - P ./ N(t)); > % drive with your rk4 on [0 40]; case 2: N = @(t) 2 + sin(t) > ``` > [!example]- Mathematica *(numerical only)* > > ```wolfram > k[t_] := 1 + Sin[t]; N0[t_] := 1; (* case 1 *) > rhs[t_, P_] := k[t] P (1 - P/N0[t]); > (* drive with your rk4 on [0,40]; case 2: N0[t_] := 2 + Sin[t] *) > ``` --- ## Problem 3 — Harvesting: two ways an equilibrium can vanish (RK4) **What to hand in.** - **By hand:** nothing required — the equilibria and critical $H$ are given in the problem. - **In code:** for each model, RK4 trajectories from a spread of $P(0)>0$ that straddle the given equilibria, at each specified $H$. - **Report:** where each trajectory ends up, which equilibrium is the attractor and which the survival threshold, and how the *sudden* (constant-yield) versus *gradual* (constant-effort) loss looks in your plots. > [!hint]- Hint 1 — the two right-hand sides > Constant yield: `rhs = k*P*(1 - P/N) - H`. Constant effort: `rhs = k*P*(1 - P/N) - H*P`. Take $k=N=1$. Everything else is choosing $P(0)$ and $H$ and watching. > [!hint]- Hint 2 — choose initial populations that bracket the equilibria > For constant yield at $H=3/16$ (equilibria $\tfrac14,\tfrac34$), run $P(0)$ below $\tfrac14$, between $\tfrac14$ and $\tfrac34$, and above $\tfrac34$. You will see $\tfrac34$ pull trajectories in while $\tfrac14$ is the knife-edge: start below it and the population dies. > [!hint]- Hint 3 — cross the critical value > Raise $H$ to $\tfrac14$ (equilibria merge at $\tfrac12$), then to $0.3$. Above critical there is *no* positive equilibrium at all, so every trajectory falls to zero — the collapse is abrupt. (Once $P$ reaches $0$, stop; negative $P$ is unphysical.) > [!hint]- Hint 4 — contrast with constant effort > Constant effort at $H=\tfrac12$ has equilibria $0$ and $\tfrac12$; every positive start rises or falls to $\tfrac12$. As $H\to1$ that stable level slides continuously to $0$. Same endpoint — extinction — reached with warning, not without. > [!example]- Python > > ```python > k = N = 1.0 > yield_rhs = lambda t, P, H: k*P*(1 - P/N) - H > effort_rhs = lambda t, P, H: k*P*(1 - P/N) - H*P > # for each H, loop P0 in e.g. [0.1, 0.2, 0.3, 0.5, 0.9]; drive your rk4 on [0, 40] > # clip at 0: once P <= 0, hold it at 0 (extinct) > ``` > *Look for:* yield $H=3/16$ splits at $P_0=\tfrac14$; yield $H=0.3$ sends everyone to 0; effort $H=\tfrac12$ pulls all to $\tfrac12$. > [!example]- R > > ```r > k <- 1; N <- 1 > yield_rhs <- function(t, P, H) k*P*(1 - P/N) - H > effort_rhs <- function(t, P, H) k*P*(1 - P/N) - H*P > # loop over H and P0; drive your rk4; clip P at 0 > ``` > [!example]- MATLAB > > ```matlab > k = 1; N = 1; > yieldRhs = @(t,P,H) k*P*(1 - P/N) - H; > effortRhs = @(t,P,H) k*P*(1 - P/N) - H*P; > % loop over H and P0; drive your rk4; clip P at 0 > ``` > [!example]- Mathematica *(numerical only)* > > ```wolfram > k = 1; N0 = 1; > yieldRhs[t_, P_, H_] := k P (1 - P/N0) - H; > effortRhs[t_, P_, H_] := k P (1 - P/N0) - H P; > (* loop over H and P0; drive your rk4; clip P at 0 *) > ``` --- ## Problem 4 — The mass-spring system in the phase plane (RKF) **What to hand in.** - **By hand:** state the condition that separates the four regimes (compare $\gamma$ to $2\sqrt{mk}$) and which $\gamma$ you use for each. - **In code:** build the first-order system; integrate the four undriven regimes with your **RKF** and the two driven cases with your fixed-step **RK4** (see Hint 4 for why), and plot each in the $(y,v)$ phase plane. - **Report:** the six phase-plane figures with a sentence each; identify which driven run resonates and which beats. > [!hint]- Hint 1 — turn it into a system your RKF can eat > Let $\vec Y=[\,y,\,v\,]$. Then $\vec Y'=[\,v,\ (f(t)-k y-\gamma v)/m\,]$. Your RKF integrates systems, so hand it this vector-valued RHS and initial vector $[y(0),v(0)]$; it returns times and a two-column solution (position and velocity). > [!hint]- Hint 2 — the four undriven regimes and their portraits > Fix $m=k=1$, so the divider is $\gamma=2\sqrt{mk}=2$: undamped $\gamma=0$, underdamped $\gamma=1$, critically damped $\gamma=2$, overdamped $\gamma=3$. In the phase plane ($v$ vs. $y$): undamped is a closed loop (energy conserved), underdamped spirals inward, critical and overdamped fall straight to the origin without looping. > [!hint]- Hint 3 — resonance vs. beating (driven, $f(t)=\cos\omega t$) > Drive the *undamped* oscillator ($\gamma=0$, natural frequency $\omega_0=\sqrt{k/m}=1$). At $\omega=\omega_0$ the amplitude grows without bound — **resonance** (the exact response is $y=\tfrac{t}{2}\sin t$). At $\omega$ near but not equal to $\omega_0$ (say $0.9$) the amplitude swells and collapses in a slow throb — **beating**. Both show clearly in $y(t)$; the phase plane spirals outward (resonance) or fills a band (beating). > [!hint]- Hint 4 — why the driven runs use RK4, not your RKF > Pure resonance grows forever, so an absolute-tolerance adaptive method (your RKF) keeps *shrinking* its step as the amplitude climbs and can quit with "min h exceeded." A **fixed-step RK4 has no step controller to starve**, so it sails right through — just pick $h$ well under the period $2\pi$ ($h=0.01$ tracks the exact $\tfrac{t}{2}\sin t$ to $\sim10^{-8}$ over $t\le50$). So: RKF for the four undriven portraits, where adaptivity is a clean win; your RK4 for the two driven runs. (This contrast — adaptive methods struggle with unbounded growth, fixed-step ones don't — is itself worth a sentence in your write-up.) > [!example]- Python > > ```python > import numpy as np > def spring(m, k, gamma, fext): > return lambda t, Y: np.array([Y[1], (fext(t) - k*Y[0] - gamma*Y[1]) / m]) > > undamped = spring(1, 1, 0, lambda t: 0.0) # undriven: gamma = 0,1,2,3 -> your RKF > resonance = spring(1, 1, 0, lambda t: np.cos(t)) # driven, omega = omega0 = 1 -> your RK4 > beating = spring(1, 1, 0, lambda t: np.cos(0.9*t))# driven, omega = 0.9 -> your RK4 > # undriven: t, Y = yourRKF(undamped, 0, 20, [1.0, 0.0], 1e-6) # adaptive, non-uniform t > # driven: t, Y = yourRK4(resonance, 0, 50, [0.0, 0.0], h=0.01) # fixed step, no min-h trouble > # phase plane: plt.plot(Y[:,0], Y[:,1]) > ``` > *Look for:* undamped traces a closed loop; resonance spirals outward; the RKF returns *non-uniform* times (fine for plotting), the RK4 a uniform grid. > [!example]- R > > ```r > spring <- function(m, k, gamma, fext) > function(t, Y) c(Y[2], (fext(t) - k*Y[1] - gamma*Y[2]) / m) > resonance <- spring(1, 1, 0, function(t) cos(t)) > # undriven -> your RKF on [0,20], tol ~ 1e-6; driven -> your RK4 on [0,50], h = 0.01 > # phase plane: plot(Y[,1], Y[,2], type = "l") > ``` > [!example]- MATLAB > > ```matlab > spring = @(m,k,g,fext) @(t,Y) [Y(2); (fext(t) - k*Y(1) - g*Y(2))/m]; > resonance = spring(1,1,0,@(t) cos(t)); % omega = omega0 = 1 > % undriven: [t,Y] = myRKFSys(undamped, [0 20], [1 0], 1e-6); % adaptive > % driven: [t,Y] = myRK4(resonance, [0 50], [0 0], 0.01); % fixed step h = 0.01 > % plot(Y(:,1), Y(:,2)) % phase plane > ``` > *Look for:* the four RKF portraits are clean; the RK4-driven resonance matches the exact $\tfrac{t}{2}\sin t$ envelope. > [!example]- Mathematica *(numerical only)* > > ```wolfram > spring[m_, k_, g_, fext_] := Function[{t, Y}, {Y[[2]], (fext[t] - k Y[[1]] - g Y[[2]])/m}]; > resonance = spring[1, 1, 0, Cos[#] &]; > (* undriven -> your RKF on {0,20}, tol ~ 10^-6; driven -> your RK4 on {0,50}, h = 0.01 *) > (* phase plane: ListLinePlot of {y, v} *) > ``` --- ## Problem 5 — The Duffing spring (built-in solver) **What to hand in.** - **By hand:** the physical intuition — soft vs. hard spring, and the double-well picture for part (c). - **In code:** integrate each part with a **built-in** ODE solver (this one is nonlinear). - **Report:** the $y(t)$ and phase-plane plots, your read on periodicity/sinusoid in (a), and the chaos discussion in (d). > [!hint]- Hint 1 — same state vector, built-in engine > Use $[\,y,\,v\,]$ again with $\vec Y'=[\,v,\ (f(t)-\gamma v-k_1 y-k_3 y^3)/m\,]$, and hand it to your language's built-in solver (`solve_ivp`, `ode45`, `deSolve::ode`, `NDSolve`). The nonlinear $k_3 y^3$ term is why we reach for the built-in rather than a fixed-step method. > [!hint]- Hint 2 — why (a) looks strange, and that's the answer > With the soft spring $k_3=-0.5$ and $y'(0)=1$, the starting energy sits *right at* the top of the potential barrier (at $y=\sqrt2\approx1.414$). So the mass creeps out toward $1.414$ and nearly freezes there rather than swinging back — it does **not** look periodic or sinusoidal, and saying so (with the plot) is the point. Drop to $k_3=-0.05$ and the barrier is far above the energy, so you get a nearly sinusoidal oscillation. > [!hint]- Hint 3 — (c) the double well > $k_1=-1,\,k_3=1$ makes a two-well potential with minima at $y=\pm1$. Under the drive $F\cos(1.2t)$, sweep $F\in\{0.2,0.37,0.5,0.65\}$ and watch the response go from tame periodic motion to well-hopping to chaotic. The phase plane reveals the structure (closed curve $\to$ doubled loops $\to$ a smeared band). > [!hint]- Hint 4 — (d) sensitive dependence > Run the two initial conditions ($y(0)=1$ and $1.01$) on the same axes. They shadow each other briefly, then peel apart and become uncorrelated while both stay bounded — that is deterministic **chaos** (sensitive dependence on initial conditions). Name it and reflect on what it means for prediction. > [!example]- Python > > ```python > import numpy as np > from scipy.integrate import solve_ivp > def duffing(m, g, k1, k3, fext): > return lambda t, Y: [Y[1], (fext(t) - g*Y[1] - k1*Y[0] - k3*Y[0]**3) / m] > > soft = duffing(1, 0, 1, -0.5, lambda t: 0.0) # part (a) > s = solve_ivp(soft, [0, 20], [0, 1], max_step=0.01, rtol=1e-9) > # y(t) = s.y[0], v(t) = s.y[1]; phase plane: plt.plot(s.y[0], s.y[1]) > well = duffing(1, 0.3, -1, 1, lambda t: 0.37*np.cos(1.2*t)) # part (c) > ``` > *Look for:* part (a) $k_3=-0.5$ crawls to $y\approx1.414$ and flattens; $k_3=-0.05$ is nearly sinusoidal. > [!example]- R > > ```r > library(deSolve) > duffing <- function(t, Y, p) list(c(Y[2], > (p$fext(t) - p$g*Y[2] - p$k1*Y[1] - p$k3*Y[1]^3) / p$m)) > p <- list(m=1, g=0, k1=1, k3=-0.5, fext=function(t) 0) # part (a) > out <- ode(y=c(0,1), times=seq(0,20,0.01), func=duffing, parms=p) > ``` > [!example]- MATLAB > > ```matlab > duffing = @(m,g,k1,k3,fext) @(t,Y) [Y(2); (fext(t) - g*Y(2) - k1*Y(1) - k3*Y(1)^3)/m]; > soft = duffing(1,0,1,-0.5,@(t) 0); % part (a) > [t,Y] = ode45(soft, [0 20], [0; 1]); % built-in > % plot(Y(:,1), Y(:,2)) % phase plane > ``` > [!example]- Mathematica *(numerical only)* > > ```wolfram > duff[m_,g_,k1_,k3_,f_, {y0_,v0_}, tf_] := NDSolveValue[ > {m y''[t] + g y'[t] + k1 y[t] + k3 y[t]^3 == f[t], > y[0] == y0, y'[0] == v0}, y, {t, 0, tf}]; > ya = duff[1,0,1,-0.5, (0 &), {0,1}, 20]; (* part (a) *) > (* Plot[ya[t],{t,0,20}]; ParametricPlot[{ya[t], ya'[t]},{t,0,20}] *) > ``` --- *Companion to `assignment/MATH307Su26-HW3.tex`. Code here drives the engines you built in lecture (adaptive Simpson's, RK4, RKF) plus a built-in solver for the nonlinear Problem 5. The `# TODO` marks and the hand calculations are yours. Submit one written PDF plus your code files and figures.*