MATH307Su26 - From Linear Springs to Nonlinear Pendulums

# From Linear Springs to Nonlinear Pendulums > [!abstract] The project in one line > Solve the linear oscillator exactly by hand, make code reproduce its period to machine precision, then point the same verified instrument at the nonlinear pendulum and measure how its period grows with amplitude, ending with a comparison of integrators (Euler versus RK4) judged by how badly each one leaks energy. ## Introduction You are adjusting a [grandfather clock](https://en.wikipedia.org/wiki/Grandfather_clock), trying to get it to keep perfect time. You notice something strange: when the pendulum swings with small angles, it keeps excellent time, but when you accidentally give it a large push, the clock runs slow. The mathematics explains why. You have crossed from the linear world, where everything behaves predictably, into the nonlinear one, where familiar rules break down. A mass on a spring oscillates at the same frequency regardless of how far you stretch it. A [pendulum](https://en.wikipedia.org/wiki/Pendulum)'s frequency depends on its amplitude. This fundamental difference is not a curiosity: it decided the design of every precision clock for three centuries, and it marks the boundary where the linear models of this course stop telling the whole truth. This project builds that boundary with numbers. The linear oscillator has an exact pencil-and-paper solution, which makes it the perfect ground truth for testing numerical machinery. Once your quadrature and your integrators reproduce the linear answers to [machine precision](https://en.wikipedia.org/wiki/Machine_epsilon), you point them at the pendulum, where no elementary closed form exists, and let them measure the nonlinearity honestly. ![Pendulum phase portraits at 5°, 90°, and 150° with the separatrix: small swings live on near-ellipses, large swings distort toward the eye-shaped boundary](Media/pendulum_phase.png) ## Project Description The project moves in four passes. First, solve the linear [harmonic oscillator](https://en.wikipedia.org/wiki/Harmonic_oscillator) exactly by hand and run the worked example to small checkable numbers (a period of exactly $\pi$ seconds, an energy of exactly $2$ J). Second, mirror the hand work in code: derive the pendulum's exact period as an integral, evaluate it by quadrature, and confirm that the linear limit reproduces the hand period at machine precision (the anchor the numerics must pass before they are trusted). Third, run the real experiment: build [Euler](https://en.wikipedia.org/wiki/Euler_method) and [RK4](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods) integrators from scratch, simulate the pendulum, measure its period from the trajectory, and sweep the amplitude to watch the period grow; then generalize by comparing how the two integrators leak energy. Finally, run the whole study on three oscillators you choose, picked for three different nonlinearity personalities. **Foundation task:** Solve the linear oscillator by hand (solution, period, energy), derive the pendulum period integral, and make the code reproduce the hand numbers: the linear period at machine precision through the quadrature, and RK4's fourth-order convergence against the exact solution. (This is Steps 1 to 3 below.) **Application task:** Use the verified machinery to measure the pendulum's period at amplitudes from $5°$ to $170°$, check every measurement against the elliptic-integral ground truth, and quantify where the small-angle prediction fails. Then compare Euler and RK4 by their energy drift over ten periods. (This is Steps 4 and 5 plus the Generalization section below.) **Key deliverable:** The same study on three oscillators you choose with three nonlinearity personalities (linear, weakly nonlinear, strongly nonlinear): the period predicted before any computing, each measurement verified against an exact value, and a cross-oscillator ranking of nonlinearity defended from the measured numbers, not the pictures. Two companions support the milestone forms: the [[MATH307Su26 - From Linear Springs to Nonlinear Pendulums (Milestone Map)|Milestone Map]] says which artifact from this handout answers each form field, and the [[MATH307Su26 - From Linear Springs to Nonlinear Pendulums (Verification Table)|Verification Table]] is filled in as you verify and submitted with Milestone 2. ## Mathematical Background: energy, phase space, and the exact period **Idea.** The linear oscillator can be solved completely by hand, so it supplies exact numbers to test code against; the pendulum cannot, but [conservation of energy](https://en.wikipedia.org/wiki/Conservation_of_energy) still turns its period into a definite integral we can evaluate by quadrature. One system calibrates the instrument, the other one needs it. ### The linear foundation **Recall.** [Hooke's law](https://en.wikipedia.org/wiki/Hooke%27s_law) says a spring stretched by $y$ pulls back with force $-ky$. [Newton's second law](https://en.wikipedia.org/wiki/Newton%27s_laws_of_motion) then gives the equation of motion of a mass $m$ on that spring, $m\ddot{y} + ky = 0 \quad\Longleftrightarrow\quad \ddot{y} = -\omega_0^2\, y, \qquad \omega_0 = \sqrt{\frac{k}{m}}.$ This is a linear equation, and it rewards us with a complete solution. Trying $y=\cos(\omega t)$ gives $-\omega^2\cos(\omega t) = -\omega_0^2\cos(\omega t)$, which works exactly when $\omega=\omega_0$; the same for $\sin(\omega_0 t)$; and linearity lets us combine them to match any start. With initial position $y_0$ and initial velocity $v_0$, $y(t) = y_0\cos(\omega_0 t) + \frac{v_0}{\omega_0}\sin(\omega_0 t),$ which repeats when $\omega_0 t$ advances by $2\pi$. The period is therefore $\boxed{\ T = \frac{2\pi}{\omega_0} = 2\pi\sqrt{\frac{m}{k}}\,,\ }$ and here is the point that the whole project turns on: **the amplitude does not appear**. This is [simple harmonic motion](https://en.wikipedia.org/wiki/Simple_harmonic_motion), and every swing, large or small, takes the same time. Total mechanical energy is conserved along the motion: $E = \frac{1}{2}m\dot{y}^2 + \frac{1}{2}ky^2 = \text{constant},$ (differentiate and substitute the equation of motion: $\dot E = \dot y\,(m\ddot y + ky) = 0$). Fixing $E$ fixes an ellipse in the [phase plane](https://en.wikipedia.org/wiki/Phase_space) $(y, \dot{y})$, $\frac{y^2}{2E/k} + \frac{\dot{y}^2}{2E/m} = 1,$ so the motion traces perfect closed ellipses at the single frequency $\omega_0$. **Note.** Ask yourself before moving on (this was the Su25 check and it is still the right one): why must undamped oscillator trajectories be closed curves in phase space, and what would an open, spiraling trajectory mean physically? (Answer in energy language: a closed curve is a level set of $E$; a spiral means $E$ is changing, which is exactly the diagnostic the Generalization section uses against numerical integrators.) ### Worked example, run to completion Take $m=1$ kg, $k=4$ N/m, released from rest at $y_0=1$ m (so $v_0=0$). Then $\omega_0=\sqrt{4/1}=2\ \text{rad/s},\qquad T=\frac{2\pi}{2}=\pi\approx3.141593\ \text{s},\qquad E=\tfrac12 k y_0^2 = 2\ \text{J},$ and the solution is $y(t)=\cos(2t)$, $v(t)=-2\sin(2t)$. One quarter period after release, $t=T/4=\pi/4$: $y(T/4)=\cos(\pi/2)=0,\qquad v(T/4)=-2\sin(\pi/2)=-2\ \text{m/s},$ and the energy check closes by hand: all the energy is kinetic there, $\tfrac12(1)(-2)^2 = 2$ J $=E$. These four numbers ($\omega_0=2$, $T=\pi$, the quarter-period state $(0,-2)$, and $E=2$) are what the warm-up code must reproduce. ![The worked linear oscillator: sinusoidal time evolution with period pi seconds, a closed phase-plane ellipse, and kinetic/potential energy trading places under a constant total of 2 J](Media/pendulum_linear.png) ### The nonlinear reality: the pendulum A mass on a rigid rod of length $L$ swinging through angle $\theta$ feels gravity's restoring [torque](https://en.wikipedia.org/wiki/Torque) $-mgL\sin\theta$ against its rotational inertia $mL^2$: $mL^2\ddot{\theta} = -mgL\sin\theta \quad\Longrightarrow\quad \ddot{\theta} + \frac{g}{L}\sin\theta = 0.$ The restoring term is $\sin\theta$, not $\theta$: the equation is nonlinear. Its conserved energy (kinetic plus gravitational, measuring height as $L(1-\cos\theta)$) is $E = \frac{1}{2}mL^2\dot{\theta}^2 + mgL\,(1-\cos\theta).$ **Recall.** The [Taylor series](https://en.wikipedia.org/wiki/Taylor_series) $\sin\theta = \theta - \theta^3/6 + \cdots$ says that for small angles $\sin\theta\approx\theta$ (the [small-angle approximation](https://en.wikipedia.org/wiki/Small-angle_approximation)). Making that replacement turns the pendulum equation into the harmonic oscillator with $\omega_0 = \sqrt{\frac{g}{L}}, \qquad T_0 = 2\pi\sqrt{\frac{L}{g}},$ amplitude-independent again. For our worked pendulum ($g=9.81$ m/s², $L=1$ m): $\omega_0 = 3.132092$ rad/s and $T_0 = 2.006067$ s. But the replacement is an approximation, and the whole question is what its error does as the amplitude grows. **The amplitude dependence.** For a swing of amplitude $\theta_0$ the true period expands as (see [pendulum (mechanics)](https://en.wikipedia.org/wiki/Pendulum_(mechanics)) for the derivation via the elliptic integral below) $\boxed{\ T(\theta_0) = T_0\left(1 + \frac{\theta_0^2}{16} + \frac{11\,\theta_0^4}{3072} + \cdots\right),\ }$ so the leading correction $\theta_0^2/16$ acts as this project's **dial**: read it first and it predicts the fractional period increase before you compute anything. At $\theta_0 = 90° = \pi/2$ the dial reads $(\pi/2)^2/16 = 0.154213$, predicting a period about $15\%$ long; the exact answer (below) is $18.03\%$ longer, i.e., the frequency drops by $15.28\%$. This amplitude-frequency coupling destroys the concept of a single "natural frequency," and it is why a hard-pushed grandfather clock runs slow. > [!note] Where this comes from > Nothing here is new machinery: the equations of motion are supplied physics, the small-angle step is the Taylor series from class, the exact period below is a definite integral evaluated by the quadrature rules from class, and Euler/RK4 are the ODE steppers whose error orders come from Taylor analysis. The project is the assembly. ### The exact period: energy conservation becomes an integral Release the pendulum from rest at $\theta_0$. Conservation of energy between the release point and a later angle $\theta$ gives $\frac{1}{2}\dot{\theta}^2 = \frac{g}{L}\left(\cos\theta - \cos\theta_0\right) \quad\Longrightarrow\quad dt = \frac{d\theta}{\omega_0\sqrt{2(\cos\theta-\cos\theta_0)}},$ and integrating over the quarter swing from $0$ to $\theta_0$ (then multiplying by four, by symmetry), $T(\theta_0) = \frac{4}{\omega_0}\int_0^{\theta_0} \frac{d\theta}{\sqrt{2\left(\cos\theta - \cos\theta_0\right)}}.$ As written this integrand blows up at $\theta=\theta_0$ (the pendulum lingers at the turning point), which is terrible for quadrature. The classical substitution $\sin\varphi = \sin(\theta/2)\,/\sin(\theta_0/2)$ absorbs the singularity and leaves the [complete elliptic integral of the first kind](https://en.wikipedia.org/wiki/Elliptic_integral): $\boxed{\ T(\theta_0) = \frac{4}{\omega_0}\,K(k), \qquad k=\sin\frac{\theta_0}{2}, \qquad K(k)=\int_0^{\pi/2}\frac{d\varphi}{\sqrt{1-k^2\sin^2\varphi}}\,.\ }$ This is the project's ground truth: a smooth, bounded integrand on a fixed interval, tailor-made for the quadrature rules from class. And it carries its own built-in exactness check: at zero amplitude $k=0$, the integrand is the constant $1$, $K(0)=\pi/2$ exactly, and the formula collapses to $T = 2\pi/\omega_0$, the linear period. Your quadrature must reproduce that limit **at machine precision**, and for the worked spring ($\omega_0=2$) that means recovering $T=\pi$ to sixteen digits. That is the mirror anchor. > [!note] Why plain trapezoid is spectacular here > The integrand of $K(k)$ is even about both endpoints ($\varphi=0$ and $\varphi=\pi/2$), so every odd derivative vanishes there and the [Euler-Maclaurin](https://en.wikipedia.org/wiki/Euler%E2%80%93Maclaurin_formula) correction terms all cancel: the humble [trapezoid rule](https://en.wikipedia.org/wiki/Trapezoidal_rule) converges far faster than its usual $O(h^2)$, reaching machine precision with a few hundred points. This is the same superconvergence surprise the Orbital Mechanics project met on full-orbit integrals; here we exploit it on purpose. **Hand warm-ups (do these before coding; the warm-up code step reproduces every one).** 1. The worked spring: $\omega_0=2$ rad/s, $T=\pi\approx3.141593$ s, $E=2$ J, quarter-period state $(y,v)=(0,-2)$. 2. The worked pendulum: $\omega_0=\sqrt{9.81}=3.132092$ rad/s, $T_0=2.006067$ s. 3. The dial at $90°$: $1+(\pi/2)^2/16 = 1.154213$. 4. The linear limit: $K(0)=\pi/2$, so $(4/\omega_0)K(0) = 2\pi/\omega_0$ exactly. ## Implementation Guidelines The workflow is the same in every language, so we state it once and then implement it step by step. > [!abstract] The workflow (state it once, run it in any language) > Given an oscillator (spring $(m,k)$ or pendulum $(g,L)$), an amplitude $\theta_0$, and a step size $h$: > 1. **Warm up**: define the verification helper, then reproduce the hand numbers ($\omega_0=2$, $T=\pi$, the quarter-period state $(0,-2)$; pendulum $\omega_0=3.132092$ rad/s, $T_0=2.006067$ s; the dial $1+\theta_0^2/16=1.154213$ at $90°$). > 2. **Period integral**: composite trapezoid for $K(k)=\int_0^{\pi/2}\bigl(1-k^2\sin^2\varphi\bigr)^{-1/2}d\varphi$; the linear limit $T=(4/\omega_0)K(0)=2\pi/\omega_0$ must come out at machine precision. > 3. **Integrators from scratch**: Euler $u_{n+1}=u_n+h\,f(u_n)$ and RK4 on the first-order system $u=(y,\dot y)$; RK4 against the exact $y=\cos 2t$, and halving $h$ must cut the error by $16$. > 4. **Simulate and measure**: release the pendulum from rest at $\theta_0$; the first two zero crossings sit at $T/4$ and $3T/4$, so $T=2(t_2-t_1)$; verify against $(4/\omega_0)K(\sin(\theta_0/2))$; sweep the amplitudes. > 5. **Render**: the period-versus-amplitude curve (exact, series, measured) and the phase portraits. > 6. **Generalize**: energy drift, Euler versus RK4 (its own section below). > 7. **Your own oscillators**: three nonlinearity personalities, prediction first (the finale). Pick your language below; each callout contains Steps 1 to 5 in order, run top to bottom. Code is four-language and portable (Python is the run-verified reference; Mathematica is numerical only). This project ships no data files: every oscillator is a couple of numbers, so there is nothing to download. > [!example]- Python (reference) > > **Step 1: define the helper, then warm up on the hand numbers.** The helper compares two values (or arrays) and reports the largest difference; every check in this project runs through it. > ```python > import numpy as np > import matplotlib.pyplot as plt > def same(name, X, Y, tol=1e-9): # verification helper > d = np.max(np.abs(np.asarray(X) - np.asarray(Y))) > print(f"{name}: max difference {d:.1e}", "OK" if d < tol else "CHECK THIS") > m, k = 1.0, 4.0 # the worked spring > w0 = np.sqrt(k/m); T_lin = 2*np.pi/w0 # 2 rad/s, pi seconds > same("omega0 equals 2", w0, 2.0) > same("quarter-period state equals (0, -2)", > [np.cos(w0*T_lin/4), -w0*np.sin(w0*T_lin/4)], [0.0, -2.0]) > g, L = 9.81, 1.0 # the worked pendulum > w0p = np.sqrt(g/L); T0p = 2*np.pi/w0p > print("pendulum omega0:", round(w0p, 6), " T0:", round(T0p, 6)) # 3.132092 2.006067 > print("period dial at 90 deg:", round(1 + (np.pi/2)**2/16, 6)) # 1.154213 > ``` > *Look for:* both helper checks OK at roundoff ($\sim10^{-16}$); pendulum $\omega_0=3.132092$ rad/s, $T_0=2.006067$ s; the dial $1.154213$. > > **Step 2: the period integral by quadrature (the mirror anchor).** Composite trapezoid on $[0,\pi/2]$; at $k=0$ the integrand is the constant $1$, so the linear period must come out at machine level, not merely close. > ```python > def K_trap(kmod, n=400): # complete elliptic integral K(k) > phi = np.linspace(0, np.pi/2, n+1) > y = 1.0/np.sqrt(1.0 - (kmod*np.sin(phi))**2) > h = (np.pi/2)/n # trapezoid superconverges here > return h*(0.5*y[0] + y[1:-1].sum() + 0.5*y[-1]) > T_num = (4/w0)*K_trap(0.0) # linear limit: K(0) = pi/2 exactly > same("linear period at machine precision", T_num, np.pi, tol=1e-13) > k90 = np.sin(np.pi/4) # modulus at theta0 = 90 deg > print("K(sin 45 deg):", round(K_trap(k90), 6)) # 1.854075 > same("quadrature self-consistency (n=400 vs 800)", K_trap(k90), K_trap(k90, 800), tol=1e-13) > from scipy.special import ellipk # library cross-check; takes m = k^2 > same("quadrature K matches scipy ellipk", K_trap(k90), ellipk(k90**2), tol=1e-12) > T_exact = lambda th0, w: (4.0/w)*K_trap(np.sin(th0/2)) > T90 = T_exact(np.pi/2, w0p) > print("T(90 deg):", round(T90, 6), " ratio T/T0:", round(T90/T0p, 6)) # 2.367842 1.180341 > ``` > *Look for:* the linear-period difference at $\sim4\times10^{-16}$ (machine precision, the anchor); $K(\sin45°)=1.854075$; all three helper checks OK; $T(90°)=2.367842$ s, ratio $1.180341$. > > **Step 3: build Euler and RK4 from scratch, then make RK4 prove its order.** The integrators live on the first-order system $u=(y,\dot y)$; the exact solution $\cos(2t)$ grades RK4, and halving $h$ must cut the error by $16$. > ```python > def euler_step(f, t, u, h): > return u + h*f(t, u) > def rk4_step(f, t, u, h): > k1 = f(t, u) > k2 = f(t + h/2, u + h/2*k1) > k3 = f(t + h/2, u + h/2*k2) > k4 = f(t + h, u + h*k3) > return u + h/6*(k1 + 2*k2 + 2*k3 + k4) > def integrate(step, f, u0, t_end, h): > n = int(round(t_end/h)); ts = np.arange(n+1)*h > us = np.empty((n+1, 2)); us[0] = u0 > u = np.array(u0, dtype=float) > for i in range(n): > u = step(f, ts[i], u, h); us[i+1] = u > return ts, us > f_lin = lambda t, u: np.array([u[1], -(k/m)*u[0]]) # the spring, u = (y, v) > for h in (0.01, 0.005): > ts, us = integrate(rk4_step, f_lin, [1.0, 0.0], T_lin, h) > print(h, np.max(np.abs(us[:, 0] - np.cos(w0*ts)))) # 6.441e-09 then 4.019e-10 > ``` > *Look for:* max errors $6.44\times10^{-9}$ at $h=0.01$ and $4.02\times10^{-10}$ at $h=0.005$: the ratio $16.03$ is the fourth order showing itself. > > **Step 4: simulate the pendulum, measure the period, sweep the amplitude.** Released from rest, $\theta(t)$ crosses zero at $T/4$ and $3T/4$, so $T=2(t_2-t_1)$; linear interpolation at the crossing is unusually accurate because $\ddot\theta = 0$ exactly where $\theta=0$. > ```python > f_pend = lambda t, u: np.array([u[1], -(g/L)*np.sin(u[0])]) # u = (theta, omega) > def measured_period(theta0, h=0.001): > ts, us = integrate(rk4_step, f_pend, [theta0, 0.0], 0.9*T_exact(theta0, w0p), h) > th = us[:, 0] > idx = np.where(np.sign(th[:-1]) != np.sign(th[1:]))[0] > cross = [ts[i] + h*th[i]/(th[i] - th[i+1]) for i in idx[:2]] > return 2*(cross[1] - cross[0]) > for deg in (5, 90): > th0 = np.deg2rad(deg) > same(f"measured period matches the exact value ({deg} deg)", > measured_period(th0), T_exact(th0, w0p), tol=1e-5) > for deg in (5, 30, 60, 90, 120, 150, 170): # the amplitude sweep > th0 = np.deg2rad(deg) > print(deg, round(T_exact(th0, w0p)/T0p, 6), round(1 + th0**2/16, 6)) > r2, r4 = (T_exact(np.deg2rad(d), w0p)/T0p for d in (2, 4)) > print("clock loses", round((r4 - r2)*86400, 1), "s/day") # 19.7 > ``` > *Look for:* measured periods $2.007022$ s ($5°$) and $2.367842$ s ($90°$) matching quadrature to $\sim10^{-10}$; exact ratios $1.000476,\ 1.017409,\ 1.073182,\ 1.180341,\ 1.372881,\ 1.762204,\ 2.439363$; the series good to $-0.03\%$ at $30°$ but off by $-18.94\%$ at $150°$; the clock number $19.7$ s/day. > > **Step 5: render the period curve and the phase portraits.** The period figure carries all three voices (exact quadrature, small-angle series, RK4 measurements); the phase figure shows the ellipses distorting toward the [separatrix](https://en.wikipedia.org/wiki/Separatrix_(mathematics)) $\dot\theta=\pm2\omega_0\cos(\theta/2)$. > ```python > degs = np.arange(1, 176) > ratios = [T_exact(np.deg2rad(d), w0p)/T0p for d in degs] > plt.plot(degs, ratios, label="exact (quadrature)") > plt.plot(degs, 1 + np.deg2rad(degs)**2/16, "--", label="series $1+\\theta_0^2/16quot;) > meas = [(d, measured_period(np.deg2rad(d))/T0p) for d in (5, 30, 60, 90, 120, 150, 170)] > plt.plot(*zip(*meas), "o", label="measured (RK4)") > plt.xlabel("amplitude theta_0 (degrees)"); plt.ylabel("period ratio T/T_0") > plt.legend(); plt.show() > for deg in (5, 90, 150): # phase portraits > th0 = np.deg2rad(deg) > ts, us = integrate(rk4_step, f_pend, [th0, 0.0], 1.02*T_exact(th0, w0p), 0.001) > plt.plot(np.rad2deg(us[:, 0]), us[:, 1], label=f"{deg} deg") > thsep = np.linspace(-np.pi, np.pi, 400) > plt.plot(np.rad2deg(thsep), 2*w0p*np.cos(thsep/2), "k--", label="separatrix") > plt.plot(np.rad2deg(thsep), -2*w0p*np.cos(thsep/2), "k--") > plt.xlabel("angle theta (degrees)"); plt.ylabel("angular velocity (rad/s)") > plt.legend(); plt.show() > ``` > *Look for:* the measured points sitting exactly on the exact curve while the dashed series peels away past $60°$; the $5°$ orbit a tiny ellipse, the $150°$ orbit hugging the eye-shaped separatrix. > [!example]- MATLAB > > **Step 1: define the helper, then warm up on the hand numbers.** MATLAB note: the helper is an anonymous function; the functions of Steps 2 to 4 are *local functions* and must sit at the END of your script file. > ```matlab > same = @(name, X, Y) fprintf('%s: max difference %.1e\n', name, max(abs(X - Y), [], 'all')); > m = 1.0; k = 4.0; % the worked spring > w0 = sqrt(k/m); T_lin = 2*pi/w0; % 2 rad/s, pi seconds > same('omega0 equals 2', w0, 2.0) > same('quarter-period state equals (0, -2)', ... > [cos(w0*T_lin/4), -w0*sin(w0*T_lin/4)], [0, -2]) > g = 9.81; L = 1.0; % the worked pendulum > w0p = sqrt(g/L); T0p = 2*pi/w0p; > fprintf('pendulum omega0: %.6f T0: %.6f\n', w0p, T0p) % 3.132092 2.006067 > fprintf('period dial at 90 deg: %.6f\n', 1 + (pi/2)^2/16) % 1.154213 > ``` > *Look for:* both helper differences at roundoff ($\sim10^{-16}$); $\omega_0=3.132092$ rad/s, $T_0=2.006067$ s; the dial $1.154213$. > > **Step 2: the period integral by quadrature (the mirror anchor).** `ellipke` takes the parameter $m=k^2$, the same convention as scipy. > ```matlab > T_num = (4/w0)*K_trap(0, 400); % linear limit: K(0) = pi/2 exactly > same('linear period at machine precision', T_num, pi) > k90 = sin(pi/4); % modulus at theta0 = 90 deg > fprintf('K(sin 45 deg): %.6f\n', K_trap(k90, 400)) % 1.854075 > same('quadrature self-consistency (n=400 vs 800)', K_trap(k90, 400), K_trap(k90, 800)) > same('quadrature K matches ellipke', K_trap(k90, 400), ellipke(k90^2)) > T90 = T_exact(pi/2, w0p); > fprintf('T(90 deg): %.6f ratio T/T0: %.6f\n', T90, T90/T0p) % 2.367842 1.180341 > > % local functions (END of the script file): > function K = K_trap(kmod, n) % complete elliptic integral K(k) > phi = linspace(0, pi/2, n+1); > y = 1 ./ sqrt(1 - (kmod*sin(phi)).^2); > h = (pi/2)/n; % trapezoid superconverges here > K = h*(0.5*y(1) + sum(y(2:end-1)) + 0.5*y(end)); > end > function T = T_exact(theta0, w) > T = (4/w)*K_trap(sin(theta0/2), 400); > end > ``` > *Look for:* the linear-period difference at $\sim10^{-16}$ (machine precision); $K(\sin45°)=1.854075$; $T(90°)=2.367842$ s, ratio $1.180341$. > > **Step 3: build Euler and RK4 from scratch, then make RK4 prove its order.** State vectors are columns; the elementwise operators (`.^`, `./`) matter. > ```matlab > f_lin = @(t, u) [u(2); -(k/m)*u(1)]; % the spring, u = (y, v) > hs = [0.01, 0.005]; errs = zeros(1, 2); > for i = 1:2 > [ts, us] = integrate_ode(@rk4_step, f_lin, [1; 0], T_lin, hs(i)); > errs(i) = max(abs(us(:, 1) - cos(w0*ts))); > fprintf('h=%g: max error %.3e\n', hs(i), errs(i)); % 6.441e-09 then 4.019e-10 > end > fprintf('ratio: %.2f\n', errs(1)/errs(2)) % 16.03 > > % local functions (END of the script file): > function u = euler_step(f, t, u, h) > u = u + h*f(t, u); > end > function u = rk4_step(f, t, u, h) > k1 = f(t, u); k2 = f(t + h/2, u + h/2*k1); > k3 = f(t + h/2, u + h/2*k2); k4 = f(t + h, u + h*k3); > u = u + h/6*(k1 + 2*k2 + 2*k3 + k4); > end > function [ts, us] = integrate_ode(step, f, u0, t_end, h) > n = round(t_end/h); ts = (0:n)'*h; > us = zeros(n+1, 2); us(1, :) = u0(:)'; u = u0(:); > for i = 1:n > u = step(f, ts(i), u, h); us(i+1, :) = u'; > end > end > ``` > *Look for:* errors $6.44\times10^{-9}$ and $4.02\times10^{-10}$, ratio $16.03$ (fourth order). > > **Step 4: simulate the pendulum, measure the period, sweep the amplitude.** > ```matlab > f_pend = @(t, u) [u(2); -(g/L)*sin(u(1))]; % u = (theta, omega) > for deg = [5, 90] > th0 = deg2rad(deg); > same(sprintf('measured period matches the exact value (%d deg)', deg), ... > measured_period(th0, 0.001, w0p, f_pend), T_exact(th0, w0p)) > end > for deg = [5, 30, 60, 90, 120, 150, 170] % the amplitude sweep > th0 = deg2rad(deg); > fprintf('%3d %.6f %.6f\n', deg, T_exact(th0, w0p)/T0p, 1 + th0^2/16); > end > r2 = T_exact(deg2rad(2), w0p)/T0p; r4 = T_exact(deg2rad(4), w0p)/T0p; > fprintf('clock loses %.1f s/day\n', (r4 - r2)*86400) % 19.7 > > % local function (END of the script file): > function T = measured_period(theta0, h, w, f_pend) > [ts, us] = integrate_ode(@rk4_step, f_pend, [theta0; 0], 0.9*T_exact(theta0, w), h); > th = us(:, 1); > idx = find(sign(th(1:end-1)) ~= sign(th(2:end)), 2); > cross = ts(idx) + h*th(idx)./(th(idx) - th(idx+1)); > T = 2*(cross(2) - cross(1)); > end > ``` > *Look for:* measured $2.007022$ s and $2.367842$ s matching quadrature to $\sim10^{-10}$; the sweep ratios $1.000476$ through $2.439363$; the clock number $19.7$ s/day. > > **Step 5: render the period curve and the phase portraits.** Export with `exportgraphics`, not `saveas`. > ```matlab > degs = 1:175; > ratios = arrayfun(@(d) T_exact(deg2rad(d), w0p)/T0p, degs); > figure; plot(degs, ratios); hold on > plot(degs, 1 + deg2rad(degs).^2/16, '--') > sweep = [5, 30, 60, 90, 120, 150, 170]; > meas = arrayfun(@(d) measured_period(deg2rad(d), 0.001, w0p, f_pend)/T0p, sweep); > plot(sweep, meas, 'o') > xlabel('amplitude theta_0 (degrees)'); ylabel('period ratio T/T_0') > legend('exact (quadrature)', 'series 1 + theta_0^2/16', 'measured (RK4)') > figure; hold on % phase portraits > for deg = [5, 90, 150] > th0 = deg2rad(deg); > [~, up] = integrate_ode(@rk4_step, f_pend, [th0; 0], 1.02*T_exact(th0, w0p), 0.001); > plot(rad2deg(up(:, 1)), up(:, 2)) > end > thsep = linspace(-pi, pi, 400); > plot(rad2deg(thsep), 2*w0p*cos(thsep/2), 'k--') > plot(rad2deg(thsep), -2*w0p*cos(thsep/2), 'k--') > xlabel('angle theta (degrees)'); ylabel('angular velocity (rad/s)') > legend('5 deg', '90 deg', '150 deg', 'separatrix') > ``` > *Look for:* measured points on the exact curve, the dashed series peeling away past $60°$; the $150°$ orbit hugging the separatrix. > [!example]- R > > **Step 1: define the helper, then warm up on the hand numbers.** Base R only; no packages needed anywhere in this project. > ```r > same <- function(name, X, Y) cat(sprintf("%s: max difference %.1e\n", name, max(abs(X - Y)))) > m <- 1.0; k <- 4.0 # the worked spring > w0 <- sqrt(k/m); T_lin <- 2*pi/w0 # 2 rad/s, pi seconds > same("omega0 equals 2", w0, 2.0) > same("quarter-period state equals (0, -2)", > c(cos(w0*T_lin/4), -w0*sin(w0*T_lin/4)), c(0, -2)) > g <- 9.81; L <- 1.0 # the worked pendulum > w0p <- sqrt(g/L); T0p <- 2*pi/w0p > cat(sprintf("pendulum omega0: %.6f T0: %.6f\n", w0p, T0p)) # 3.132092 2.006067 > cat(sprintf("period dial at 90 deg: %.6f\n", 1 + (pi/2)^2/16)) # 1.154213 > ``` > *Look for:* both helper differences at roundoff ($\sim10^{-16}$); $\omega_0=3.132092$ rad/s, $T_0=2.006067$ s; the dial $1.154213$. > > **Step 2: the period integral by quadrature (the mirror anchor).** Base R has no elliptic-$K$ function, so the $n$ versus $2n$ self-consistency check stands in for the library cross-check the other languages run. > ```r > K_trap <- function(kmod, n = 400) { # complete elliptic integral K(k) > phi <- seq(0, pi/2, length.out = n + 1) > y <- 1/sqrt(1 - (kmod*sin(phi))^2) > h <- (pi/2)/n # trapezoid superconverges here > h*(0.5*y[1] + sum(y[2:n]) + 0.5*y[n + 1]) > } > T_num <- (4/w0)*K_trap(0) # linear limit: K(0) = pi/2 exactly > same("linear period at machine precision", T_num, pi) > k90 <- sin(pi/4) # modulus at theta0 = 90 deg > cat(sprintf("K(sin 45 deg): %.6f\n", K_trap(k90))) # 1.854075 > same("quadrature self-consistency (n=400 vs 800)", K_trap(k90), K_trap(k90, 800)) > T_exact <- function(theta0, w) (4/w)*K_trap(sin(theta0/2)) > T90 <- T_exact(pi/2, w0p) > cat(sprintf("T(90 deg): %.6f ratio T/T0: %.6f\n", T90, T90/T0p)) # 2.367842 1.180341 > ``` > *Look for:* the linear-period difference at $\sim10^{-16}$ (machine precision); $K(\sin45°)=1.854075$; $T(90°)=2.367842$ s, ratio $1.180341$. > > **Step 3: build Euler and RK4 from scratch, then make RK4 prove its order.** > ```r > euler_step <- function(f, t, u, h) u + h*f(t, u) > rk4_step <- function(f, t, u, h) { > k1 <- f(t, u); k2 <- f(t + h/2, u + h/2*k1) > k3 <- f(t + h/2, u + h/2*k2); k4 <- f(t + h, u + h*k3) > u + h/6*(k1 + 2*k2 + 2*k3 + k4) > } > integrate_ode <- function(step, f, u0, t_end, h) { > n <- round(t_end/h); ts <- (0:n)*h > us <- matrix(0, n + 1, 2); us[1, ] <- u0; u <- u0 > for (i in 1:n) { u <- step(f, ts[i], u, h); us[i + 1, ] <- u } > list(ts = ts, us = us) > } > f_lin <- function(t, u) c(u[2], -(k/m)*u[1]) # the spring, u = (y, v) > errs <- sapply(c(0.01, 0.005), function(h) { > sol <- integrate_ode(rk4_step, f_lin, c(1, 0), T_lin, h) > max(abs(sol$us[, 1] - cos(w0*sol$ts))) > }) > print(errs) # 6.441e-09 4.019e-10 > cat("ratio:", errs[1]/errs[2], "\n") # 16.03 > ``` > *Look for:* errors $6.44\times10^{-9}$ and $4.02\times10^{-10}$, ratio $16.03$ (fourth order). > > **Step 4: simulate the pendulum, measure the period, sweep the amplitude.** > ```r > f_pend <- function(t, u) c(u[2], -(g/L)*sin(u[1])) # u = (theta, omega) > measured_period <- function(theta0, h = 0.001) { > sol <- integrate_ode(rk4_step, f_pend, c(theta0, 0), 0.9*T_exact(theta0, w0p), h) > th <- sol$us[, 1] > idx <- which(sign(th[-length(th)]) != sign(th[-1]))[1:2] > cross <- sol$ts[idx] + h*th[idx]/(th[idx] - th[idx + 1]) > 2*(cross[2] - cross[1]) > } > for (deg in c(5, 90)) { > th0 <- deg*pi/180 > same(sprintf("measured period matches the exact value (%d deg)", deg), > measured_period(th0), T_exact(th0, w0p)) > } > for (deg in c(5, 30, 60, 90, 120, 150, 170)) # the amplitude sweep > cat(deg, T_exact(deg*pi/180, w0p)/T0p, 1 + (deg*pi/180)^2/16, "\n") > r2 <- T_exact(2*pi/180, w0p)/T0p; r4 <- T_exact(4*pi/180, w0p)/T0p > cat("clock loses", (r4 - r2)*86400, "s/day\n") # 19.7 > ``` > *Look for:* measured $2.007022$ s and $2.367842$ s matching quadrature to $\sim10^{-10}$; the sweep ratios $1.000476$ through $2.439363$; the clock number $19.7$ s/day. > > **Step 5: render the period curve and the phase portraits.** Under headless `Rscript`, wrap each plot in a `draw()` closure and use `png(); draw(); dev.off()` to save (the COURSE pattern). > ```r > degs <- 1:175 > ratios <- sapply(degs, function(d) T_exact(d*pi/180, w0p)/T0p) > plot(degs, ratios, type = "l", xlab = "amplitude theta_0 (degrees)", > ylab = "period ratio T/T_0") > lines(degs, 1 + (degs*pi/180)^2/16, lty = 2) > sweep <- c(5, 30, 60, 90, 120, 150, 170) > points(sweep, sapply(sweep, function(d) measured_period(d*pi/180)/T0p), pch = 19) > legend("topleft", c("exact (quadrature)", "series 1 + theta_0^2/16", "measured (RK4)"), > lty = c(1, 2, NA), pch = c(NA, NA, 19)) > plot(NULL, xlim = c(-185, 185), ylim = c(-6.6, 6.6), # phase portraits > xlab = "angle theta (degrees)", ylab = "angular velocity (rad/s)") > for (deg in c(5, 90, 150)) { > th0 <- deg*pi/180 > sol <- integrate_ode(rk4_step, f_pend, c(th0, 0), 1.02*T_exact(th0, w0p), 0.001) > lines(sol$us[, 1]*180/pi, sol$us[, 2]) > } > thsep <- seq(-pi, pi, length.out = 400) > lines(thsep*180/pi, 2*w0p*cos(thsep/2), lty = 2) > lines(thsep*180/pi, -2*w0p*cos(thsep/2), lty = 2) > ``` > *Look for:* measured points on the exact curve, the dashed series peeling away past $60°$; the $150°$ orbit hugging the separatrix. > [!example]- Mathematica (numerical only) > > **Step 1: define the helper, then warm up on the hand numbers.** `E`, `N`, `D`, `C`, `K` are protected or reserved names, so use your own (`mm`, `kk`, `kTrap`, and friends); keep everything in machine reals. > ```wolfram > same[name_, x_, y_] := Print[name, ": max difference ", Max@Abs[Flatten[{x}] - Flatten[{y}]]]; > mm = 1.0; kk = 4.0; (* the worked spring *) > w0 = Sqrt[kk/mm]; tLin = 2. Pi/w0; (* 2 rad/s, pi seconds *) > same["omega0 equals 2", w0, 2.] > same["quarter-period state equals (0, -2)", > {Cos[w0 tLin/4], -w0 Sin[w0 tLin/4]}, {0., -2.}] > grav = 9.81; len = 1.0; (* the worked pendulum *) > w0p = Sqrt[grav/len]; t0p = 2. Pi/w0p; > Print["pendulum omega0: ", w0p, " T0: ", t0p] (* 3.132092 2.006067 *) > Print["period dial at 90 deg: ", 1. + (Pi/2.)^2/16] (* 1.154213 *) > ``` > *Look for:* both helper differences at roundoff ($\sim10^{-16}$); $\omega_0=3.132092$ rad/s, $T_0=2.006067$ s; the dial $1.154213$. > > **Step 2: the period integral by quadrature (the mirror anchor).** `EllipticK[m]` takes the parameter $m=k^2$ (evaluated numerically only, as a cross-check; the pipeline runs on our own trapezoid). > ```wolfram > kTrap[kmod_, n_: 400] := Module[{phi, y, h}, (* complete elliptic integral K(k) *) > phi = Subdivide[0., N[Pi/2], n]; > y = 1./Sqrt[1. - (kmod Sin[phi])^2]; > h = (Pi/2.)/n; (* trapezoid superconverges here *) > h (0.5 First[y] + Total[y[[2 ;; -2]]] + 0.5 Last[y])]; > tNum = (4./w0) kTrap[0.]; (* linear limit: K(0) = pi/2 exactly *) > same["linear period at machine precision", tNum, N[Pi]] > k90 = Sin[N[Pi/4]]; (* modulus at theta0 = 90 deg *) > Print["K(sin 45 deg): ", kTrap[k90]] (* 1.854075 *) > same["quadrature self-consistency (n=400 vs 800)", kTrap[k90], kTrap[k90, 800]] > same["quadrature K matches EllipticK", kTrap[k90], N[EllipticK[k90^2]]] > tExact[theta0_, w_] := (4./w) kTrap[Sin[theta0/2.]]; > t90 = tExact[N[Pi/2], w0p]; > Print["T(90 deg): ", t90, " ratio T/T0: ", t90/t0p] (* 2.367842 1.180341 *) > ``` > *Look for:* the linear-period difference at $\sim10^{-16}$ (machine precision); $K(\sin45°)=1.854075$; $T(90°)=2.367842$ s, ratio $1.180341$. > > **Step 3: build Euler and RK4 from scratch, then make RK4 prove its order.** `FoldList` carries the state through the steps; times enter as the fold's second argument. > ```wolfram > eulerStep[f_, t_, u_, h_] := u + h f[t, u]; > rk4Step[f_, t_, u_, h_] := Module[{c1, c2, c3, c4}, > c1 = f[t, u]; c2 = f[t + h/2, u + h/2 c1]; > c3 = f[t + h/2, u + h/2 c2]; c4 = f[t + h, u + h c3]; > u + h/6 (c1 + 2 c2 + 2 c3 + c4)]; > integrateODE[step_, f_, u0_, tEnd_, h_] := Module[{n = Round[tEnd/h], ts}, > ts = Range[0, n] h; > {ts, FoldList[step[f, #2, #1, h] &, N[u0], Most[ts]]}]; > fLin[t_, u_] := {u[[2]], -(kk/mm) u[[1]]}; (* the spring, u = (y, v) *) > errs = Table[Module[{ts, us}, > {ts, us} = integrateODE[rk4Step, fLin, {1., 0.}, tLin, h]; > Max[Abs[us[[All, 1]] - Cos[w0 ts]]]], {h, {0.01, 0.005}}]; > Print[errs] (* 6.441e-9, 4.019e-10 *) > Print["ratio: ", errs[[1]]/errs[[2]]] (* 16.03 *) > ``` > *Look for:* errors $6.44\times10^{-9}$ and $4.02\times10^{-10}$, ratio $16.03$ (fourth order). > > **Step 4: simulate the pendulum, measure the period, sweep the amplitude.** > ```wolfram > fPend[t_, u_] := {u[[2]], -(grav/len) Sin[u[[1]]]}; (* u = (theta, omega) *) > measuredPeriod[theta0_, h_: 0.001] := Module[{ts, us, th, idx, cross}, > {ts, us} = integrateODE[rk4Step, fPend, {theta0, 0.}, 0.9 tExact[theta0, w0p], h]; > th = us[[All, 1]]; > idx = Take[Flatten[Position[Most[th] Rest[th], _?(# < 0 &)]], 2]; > cross = ts[[#]] + h th[[#]]/(th[[#]] - th[[# + 1]]) & /@ idx; > 2 (cross[[2]] - cross[[1]])]; > Do[Module[{th0 = deg Pi/180.}, > same["measured period matches the exact value (" <> ToString[deg] <> " deg)", > measuredPeriod[th0], tExact[th0, w0p]]], {deg, {5, 90}}]; > Do[Module[{th0 = deg Pi/180.}, (* the amplitude sweep *) > Print[deg, " ", tExact[th0, w0p]/t0p, " ", 1. + th0^2/16]], > {deg, {5, 30, 60, 90, 120, 150, 170}}]; > r2 = tExact[2. Pi/180, w0p]/t0p; r4 = tExact[4. Pi/180, w0p]/t0p; > Print["clock loses ", (r4 - r2) 86400., " s/day"] (* 19.7 *) > ``` > *Look for:* measured $2.007022$ s and $2.367842$ s matching quadrature to $\sim10^{-10}$; the sweep ratios $1.000476$ through $2.439363$; the clock number $19.7$ s/day. > > **Step 5: render the period curve and the phase portraits.** `Export` saves but does not display; evaluate the plot on its own line to see it. > ```wolfram > degs = Range[1, 175]; > ratios = tExact[# Pi/180., w0p]/t0p & /@ degs; > Show[ListLinePlot[{Transpose[{degs, ratios}], > Transpose[{degs, 1. + (degs Pi/180.)^2/16}]}, PlotStyle -> {Automatic, Dashed}, > AxesLabel -> {"amplitude theta0 (degrees)", "period ratio T/T0"}], > ListPlot[Table[{deg, measuredPeriod[deg Pi/180.]/t0p}, > {deg, {5, 30, 60, 90, 120, 150, 170}}], PlotStyle -> PointSize[.015]]] > orbits = Table[Module[{th0 = deg Pi/180., ts, us}, (* phase portraits *) > {ts, us} = integrateODE[rk4Step, fPend, {th0, 0.}, 1.02 tExact[th0, w0p], 0.001]; > Transpose[{us[[All, 1]] 180./Pi, us[[All, 2]]}]], {deg, {5, 90, 150}}]; > thSep = Subdivide[-N[Pi], N[Pi], 400]; > Show[ListLinePlot[orbits, AxesLabel -> {"theta (degrees)", "angular velocity (rad/s)"}], > ListLinePlot[{Transpose[{thSep 180./Pi, 2 w0p Cos[thSep/2]}], > Transpose[{thSep 180./Pi, -2 w0p Cos[thSep/2]}]}, > PlotStyle -> {{Black, Dashed}, {Black, Dashed}}]] > ``` > *Look for:* measured points on the exact curve, the dashed series peeling away past $60°$; the $150°$ orbit hugging the separatrix. What the guided run should produce: ![Period ratio versus amplitude: the exact elliptic-integral curve, the small-angle series that peels away past 60 degrees, and the RK4-measured points sitting exactly on the exact curve (ratio 1.1803 at 90 degrees, 2.4394 at 170 degrees)](Media/pendulum_period.png) > [!warning] Verify against ground truth (required) > Record these checked numbers; "it ran" is not verification. > - Warm-up: $\omega_0=2$ rad/s, $T=\pi=3.141593$ s, $E=2$ J; the quarter-period state $(0,-2)$ through the helper at $\sim10^{-16}$; pendulum $\omega_0=3.132092$ rad/s, $T_0=2.006067$ s; the dial at $90°$ reads $1.154213$. > - The mirror anchor: $(4/\omega_0)\,K(0) = \pi$ with helper difference $\sim4\times10^{-16}$ (machine precision, not merely small); $K(\sin45°)=1.854075$; the $n=400$ versus $n=800$ self-consistency difference is $0.0$; the library cross-check agrees to $\sim2\times10^{-16}$. > - RK4 order: max errors against $\cos(2t)$ over one period are $6.441\times10^{-9}$ ($h=0.01$) and $4.019\times10^{-10}$ ($h=0.005$); the ratio $16.03$ is the fourth order. > - Measured periods: $T(5°)=2.007022$ s and $T(90°)=2.367842$ s, each matching its quadrature value to $\sim10^{-10}$ s (helper differences $6.1\times10^{-10}$ and $1.5\times10^{-10}$). > - The sweep, $T(\theta_0)/T_0$ exact: $1.000476,\ 1.017409,\ 1.073182,\ 1.180341,\ 1.372881,\ 1.762204,\ 2.439363$ at $\theta_0=5°, 30°, 60°, 90°, 120°, 150°, 170°$; the series errs by $-0.03\%$ at $30°$, $-2.21\%$ at $90°$, $-18.94\%$ at $150°$, $-36.45\%$ at $170°$. > - The clock number: ratios $1.00007616$ at $2°$ and $1.00030470$ at $4°$; a clock regulated at $2°$ but swinging at $4°$ loses $19.7$ s/day. ### Generalization: energy drift, Euler versus RK4 (Step 6) The period study trusted RK4 because it proved its order on an exact solution. Now stress the trust: energy conservation says $E(t)/E(0)=1$ forever, so any drift in a computed $E(t)$ is pure integrator error, visible without knowing the exact trajectory at all. This makes energy the sharpest cheap diagnostic for long simulations, and it is exactly the "open spiral" question from the Background made quantitative. **Prediction (before running).** Euler's global error is first order, so its energy drift after a fixed time should roughly halve when $h$ does; RK4's trajectory error is fourth order, so its drift should fall by about $16$. Euler adds energy (each straight-line step overshoots the curving orbit outward, so the phase spiral opens); RK4's tiny error turns out to be dissipative here (the spiral closes inward, very slowly). > [!example]- Step 6 code, Python (reference) > > **Step 6: integrate ten periods at $90°$ with both methods and watch $E(t)$.** Energy per unit $mL^2$ is $E=\tfrac12\dot\theta^2+\omega_0^2(1-\cos\theta)$. > ```python > th0, t_end = np.pi/2, 10*T90 > E_start = (g/L)*(1 - np.cos(th0)) # 9.81 * (1 - cos 90 deg) = 9.81 > energy = lambda us: 0.5*us[:, 1]**2 + (g/L)*(1 - np.cos(us[:, 0])) > for step, hs in ((euler_step, (0.001, 0.0005)), (rk4_step, (0.01, 0.005))): > for h in hs: > ts, us = integrate(step, f_pend, [th0, 0.0], t_end, h) > print(step.__name__, h, abs(energy(us)[-1]/E_start - 1)) > ts_e, us_e = integrate(euler_step, f_pend, [th0, 0.0], t_end, 0.001) > ts_r, us_r = integrate(rk4_step, f_pend, [th0, 0.0], t_end, 0.001) > plt.semilogy(ts_e, np.abs(energy(us_e)/E_start - 1), label="Euler, h = 0.001") > plt.semilogy(ts_r, np.maximum(np.abs(energy(us_r)/E_start - 1), 1e-17), label="RK4, h = 0.001") > plt.xlabel("time t (s)"); plt.ylabel("relative energy drift |E(t)/E_0 - 1|") > plt.legend(); plt.show() > ``` > *Look for:* Euler drift $1.602\times10^{-1}$ at $h=0.001$ and $7.892\times10^{-2}$ at $h=0.0005$ (ratio $2.03$, first order, and it *grows*); RK4 drift $1.774\times10^{-8}$ at $h=0.01$ and $5.544\times10^{-10}$ at $h=0.005$ (ratio $31.99$, and it *decays*); at the same $h=0.001$, Euler $1.602\times10^{-1}$ versus RK4 $1.615\times10^{-13}$, twelve orders of magnitude apart. > [!example]- Step 6 code, MATLAB > > **Step 6: integrate ten periods at $90°$ with both methods and watch $E(t)$.** > ```matlab > th0 = pi/2; t_end = 10*T90; > E_start = (g/L)*(1 - cos(th0)); % 9.81 > methods = {@euler_step, [0.001, 0.0005]; @rk4_step, [0.01, 0.005]}; > for j = 1:2 > for h = methods{j, 2} > [~, us] = integrate_ode(methods{j, 1}, f_pend, [th0; 0], t_end, h); > E = 0.5*us(:, 2).^2 + (g/L)*(1 - cos(us(:, 1))); > fprintf('h=%g: drift %.3e\n', h, abs(E(end)/E_start - 1)); > end > end > [ts_e, us_e] = integrate_ode(@euler_step, f_pend, [th0; 0], t_end, 0.001); > [ts_r, us_r] = integrate_ode(@rk4_step, f_pend, [th0; 0], t_end, 0.001); > E_e = 0.5*us_e(:, 2).^2 + (g/L)*(1 - cos(us_e(:, 1))); > E_r = 0.5*us_r(:, 2).^2 + (g/L)*(1 - cos(us_r(:, 1))); > figure; semilogy(ts_e, max(abs(E_e/E_start - 1), 1e-17)); hold on > semilogy(ts_r, max(abs(E_r/E_start - 1), 1e-17)) > xlabel('time t (s)'); ylabel('relative energy drift'); legend('Euler', 'RK4') > ``` > *Look for:* the same drift numbers as the Python reference (Euler ratio $2.03$, RK4 ratio $31.99$; twelve orders apart at the shared $h$). > [!example]- Step 6 code, R > > **Step 6: integrate ten periods at $90°$ with both methods and watch $E(t)$.** > ```r > th0 <- pi/2; t_end <- 10*T90 > E_start <- (g/L)*(1 - cos(th0)) # 9.81 > energy <- function(us) 0.5*us[, 2]^2 + (g/L)*(1 - cos(us[, 1])) > for (meth in list(list(euler_step, c(0.001, 0.0005)), list(rk4_step, c(0.01, 0.005)))) > for (h in meth[[2]]) { > sol <- integrate_ode(meth[[1]], f_pend, c(th0, 0), t_end, h) > E <- energy(sol$us) > cat(sprintf("h=%g: drift %.3e\n", h, abs(E[length(E)]/E_start - 1))) > } > sol_e <- integrate_ode(euler_step, f_pend, c(th0, 0), t_end, 0.001) > sol_r <- integrate_ode(rk4_step, f_pend, c(th0, 0), t_end, 0.001) > plot(sol_e$ts, pmax(abs(energy(sol_e$us)/E_start - 1), 1e-17), type = "l", log = "y", > xlab = "time t (s)", ylab = "relative energy drift", ylim = c(1e-17, 1)) > lines(sol_r$ts, pmax(abs(energy(sol_r$us)/E_start - 1), 1e-17), col = "darkorange") > legend("right", c("Euler", "RK4"), col = c("black", "darkorange"), lwd = 2) > ``` > *Look for:* the same drift numbers as the Python reference (Euler ratio $2.03$, RK4 ratio $31.99$; twelve orders apart at the shared $h$). > [!example]- Step 6 code, Mathematica (numerical only) > > **Step 6: integrate ten periods at $90°$ with both methods and watch $E(t)$.** > ```wolfram > th0 = N[Pi/2]; tEnd = 10 t90; > eStart = (grav/len) (1 - Cos[th0]); (* 9.81 *) > energyPend[us_] := 0.5 us[[All, 2]]^2 + (grav/len) (1 - Cos[us[[All, 1]]]); > Do[Do[Module[{ts, us, en}, > {ts, us} = integrateODE[meth[[1]], fPend, {th0, 0.}, tEnd, h]; > en = energyPend[us]; > Print[meth[[2]], " h=", h, ": drift ", Abs[Last[en]/eStart - 1]]], > {h, meth[[3]]}], > {meth, {{eulerStep, "Euler", {0.001, 0.0005}}, {rk4Step, "RK4", {0.01, 0.005}}}}]; > {tsE, usE} = integrateODE[eulerStep, fPend, {th0, 0.}, tEnd, 0.001]; > {tsR, usR} = integrateODE[rk4Step, fPend, {th0, 0.}, tEnd, 0.001]; > ListLogPlot[{Transpose[{tsE, Clip[Abs[energyPend[usE]/eStart - 1], {10.^-17, Infinity}]}], > Transpose[{tsR, Clip[Abs[energyPend[usR]/eStart - 1], {10.^-17, Infinity}]}]}, > Joined -> True, PlotLegends -> {"Euler", "RK4"}, > AxesLabel -> {"time t (s)", "relative energy drift"}] > ``` > *Look for:* the same drift numbers as the Python reference (Euler ratio $2.03$, RK4 ratio $31.99$; twelve orders apart at the shared $h$). What the generalization should produce: ![Relative energy drift over ten periods at 90 degrees: Euler climbs steadily to 16 percent while RK4 at the same step size holds near 1e-13, a twelve-order-of-magnitude gap](Media/pendulum_drift.png) > [!warning] Verify the generalization (required) > - Setup: $\theta_0=90°$, $E_0/mL^2 = 9.81$, run time $10\,T = 23.678$ s. > - Euler: relative drift $1.602\times10^{-1}$ at $h=0.001$ and $7.892\times10^{-2}$ at $h=0.0005$, both *growing*; ratio per halving $2.03$ (first order). > - RK4: relative drift $1.774\times10^{-8}$ at $h=0.01$ and $5.544\times10^{-10}$ at $h=0.005$, both *decaying*; ratio per halving $31.99$ (see the note below). > - Head to head at $h=0.001$: Euler $1.602\times10^{-1}$, RK4 $1.615\times10^{-13}$; RK4's drift is about $10^{12}$ times smaller at identical cost per step count. > [!note] A surprise worth seeing: RK4's energy ratio is 32, not 16 > Halving $h$ cut RK4's energy drift by $31.99$, a fifth-order signature, even though RK4's trajectory error is honestly fourth order (your Step 3 ratio was $16.03$). No contradiction: the leading term of RK4's error on an oscillation is *dissipative*, and it damps the energy one order faster than it bends the trajectory; you saw it as the slow inward spiral. Euler has no such luck: its leading error pumps energy in, first order, every step. The clean way to say it: order describes the trajectory; what an integrator does to *conserved quantities* is extra structure, and it is why long-run simulations (see Real-World Context) choose their integrators by energy behavior, not order alone. ### Your oscillators: three nonlinearity personalities Now run the study on three oscillators you choose yourself, picked in advance to have three different nonlinearity personalities. The order of operations matters every time: write down the parameters, read the dial $\theta_0^2/16$ (where it applies), predict the period, and only then compute. 1. **Oscillator A, linear.** Your own mass-spring system with $m$ and $k$ taken from something real (a car suspension corner, a bungee cord, a lab spring; cite where the numbers came from). Prediction to test: $T=2\pi\sqrt{m/k}$ exactly, and the period is *amplitude-independent*: simulate at two amplitudes that differ by a factor of ten and the two measured periods must agree with each other and with the formula through the helper (differences at the $10^{-9}$ level of the measurement, not percent level). 2. **Oscillator B, weakly nonlinear.** A pendulum with your own length $L$ (a playground swing, a chandelier, a wrecking ball) at one amplitude of your choice between $10°$ and $45°$. Prediction to test: read the dial first and state the expected period increase in percent *before* running; then measure and compare all three numbers (dial prediction, quadrature exact, RK4 measurement). Our sweep found the series good to $-0.03\%$ at $30°$; yours should behave comparably. 3. **Oscillator C, strongly nonlinear.** The same pendulum released from $150°$ or higher (state your choice; stay below $180°$). Prediction to test: the dial *fails* (at $150°$ it underpredicts by $18.94\%$, at $170°$ by $36.45\%$); the quadrature and the simulation must still agree to $\sim10^{-8}$ s; the time series should show flat-topped swings (the pendulum lingering near the turning points) and the phase orbit should hug the separatrix. For each oscillator, report: the parameters and their source; the predicted period (formula or dial) made before computing; the measured period and its exact value through the helper; and one energy-drift check (RK4, ten periods) showing your simulation is trustworthy. Close with a cross-oscillator comparison: rank the three by nonlinearity using the measured $T/T_0$ and the dial's error, and defend the ranking from those numbers, not from the pictures. #### Student task loop for implementation, analysis and reflection 1. **Predict** the period (and the dial's accuracy) from the parameters, before running. 2. **Implement** the workflow (Steps 2 to 5) on the oscillator. 3. **Compare** the measured period against the exact value and against your prediction; check the energy drift. 4. **Interpret** what the agreement (or the dial's failure) says about where linear thinking stops working for this system. ## Reflection Framework Address these (2-3 focused questions per category, no more): ### Linearity and its limits - The spring's period is amplitude-independent; the pendulum's is not. Point to the exact line in the derivation of $T=2\pi/\omega_0$ where linearity gets used, and explain what $\sin\theta\ne\theta$ does to it. - The series $1+\theta_0^2/16$ erred by $-0.03\%$ at $30°$ and $-36.45\%$ at $170°$. What, in Taylor-series terms, decides where a truncated expansion stops being trustworthy? ### Integrators and trust - Your measured periods matched the quadrature ground truth to $\sim10^{-10}$ s even though the integrator knows no pendulum theory. List what each verification layer (the machine-precision linear anchor, the RK4 order check, the crossing measurement) actually certified, and what would have slipped through if you had skipped one. - Euler's phase orbit spirals outward while RK4's stays closed for ten periods. Why is energy drift a sharper diagnostic than eyeballing $\theta(t)$, and what does it cost to compute? ### Mathematical insights - **[Separatrix](https://en.wikipedia.org/wiki/Separatrix_(mathematics)):** the curve $\dot\theta=\pm2\omega_0\cos(\theta/2)$ divides back-and-forth swinging from over-the-top rotation, and $T(\theta_0)\to\infty$ as $\theta_0\to180°$. What is the pendulum doing physically in that limit, and how did your $170°$ run ($T/T_0=2.44$) foreshadow it? - **[Euler-Maclaurin superconvergence](https://en.wikipedia.org/wiki/Euler%E2%80%93Maclaurin_formula):** the plain trapezoid rule delivered $K(k)$ to machine precision with a few hundred points. Which property of the integrand (not of the rule) explains that, and where did the same effect appear in the Orbital Mechanics project? - **Order versus structure:** RK4's energy drift fell by $32$ per halving while its trajectory error fell by $16$. What does that say about the character (dissipative versus dispersive) of its leading error term? ## (Optional) Mathematical Extensions These go beyond the project's stated level; attempt them if interested. - **Symplectic Euler in one line:** change explicit Euler to use the *updated* velocity in the position update (a [semi-implicit Euler](https://en.wikipedia.org/wiki/Semi-implicit_Euler_method) step). Rerun the drift study: the energy error stops growing and stays bounded forever. This is the doorway to [symplectic integrators](https://en.wikipedia.org/wiki/Symplectic_integrator), the standard for long-run orbital and molecular simulations. - **$K(k)$ by the [arithmetic-geometric mean](https://en.wikipedia.org/wiki/Arithmetic%E2%80%93geometric_mean):** $K = \pi/(2\,\mathrm{agm}(1, \sqrt{1-k^2}))$, converging quadratically (digits double per iteration). Implement it and race your trapezoid. - **The separatrix asymptotic:** as $k\to1$, $K(k)\sim\ln\bigl(4/\sqrt{1-k^2}\bigr)$, so the period diverges logarithmically. Check the rate against your quadrature at $178°$, $179°$, $179.5°$. - **Damping and driving:** add $-\gamma\dot\theta + F\cos(\Omega t)$ to the pendulum and watch [resonance](https://en.wikipedia.org/wiki/Resonance) and, for large driving, the onset of chaos; the sensitivity toolkit for that regime is the Lorenz project's. - **The exact trajectory:** the pendulum's $\theta(t)$ has a closed form in [Jacobi elliptic functions](https://en.wikipedia.org/wiki/Jacobi_elliptic_functions); compare it pointwise against your RK4 trajectory. ## (Optional) Real-World Context ### Applications - **[Pendulum clocks](https://en.wikipedia.org/wiki/Pendulum_clock):** [Christiaan Huygens](https://en.wikipedia.org/wiki/Christiaan_Huygens) built the first one in 1656 precisely because the period is *nearly* amplitude-independent at small swing; three centuries of horology then fought the $\theta_0^2/16$ term you measured (escapement design is amplitude regulation). - **[MEMS resonators](https://en.wikipedia.org/wiki/Microelectromechanical_systems):** micro-scale oscillators in phones and watches show the same amplitude-frequency pull at operating drive levels; designers budget it exactly the way your dial does. - **[Molecular dynamics](https://en.wikipedia.org/wiki/Molecular_dynamics):** simulations of proteins run billions of steps, where Euler-style energy pumping would boil the molecule; the field runs on energy-respecting (symplectic) integrators for exactly the reason your drift plot shows. - **[Structural dynamics](https://en.wikipedia.org/wiki/Tacoma_Narrows_Bridge_(1940)):** large-amplitude motion pushes bridges and towers out of the linear regime, where resonance calculations based on a single natural frequency stop being safe. ### Technical challenges - **[Long-time integration](https://en.wikipedia.org/wiki/Symplectic_integrator):** order alone does not guarantee faithful long-run behavior; conserving structure (energy, area in phase space) is a separate design goal, as your Euler-versus-RK4 comparison demonstrated. - **[Period measurement from data](https://en.wikipedia.org/wiki/Zero_crossing):** real signals carry noise, and zero-crossing timing degrades; your clean $10^{-10}$ s agreement is the noiseless best case. - **[Near-separatrix stiffness](https://en.wikipedia.org/wiki/Stiff_equation):** close to $180°$ the dynamics mixes a slow crawl with fast swings, and fixed-step integrators waste effort; adaptive stepping is the standard cure. ### Why it matters A grandfather clock regulated at a $2°$ swing but running at $4°$ loses $19.7$ seconds per day (your Step 4 number: the period ratio moves from $1.00007616$ to $1.00030470$). That error, one part in $4000$ from a barely visible change in swing, is why every precision pendulum clock ever built carries an amplitude-control mechanism, and why "the" natural frequency is a linear-world fiction that engineering has to budget for.