%% CAS_F0027_VERIFY.m — Simple harmonic oscillator (mẍ + kx = 0; ω² = k/m) % Pillar : Mechanics % CalRef : Math Appendix §4B–4C, Mechanics Calibration §2A % Chain : quadratic potential → force → Newton → EOM → ω² identification % Hardening: isAlways(..., 'Unknown', 'false') on every symbolic assertion % symvar-based presence checks; solver and has-function avoided clear; clc; fprintf('\n=== CAS F0027 : Simple harmonic oscillator ===\n\n'); %% ── Symbols ────────────────────────────────────────────────── syms x t real syms m k real positive syms omega real positive syms A_amp phi_0 real % amplitude and phase % Position as symfun of t syms x_t(t) %% ── A. Inputs ──────────────────────────────────────────────── % Potential: V(x) = V₀ + (1/2)kx² (V₀ drops out of force) syms V0 real V_x = V0 + k * x^2 / 2; %% ── Step 1: Force from potential ──────────────────────────── % F = −dV/dx = −kx dVdx = diff(V_x, x); F_x = -dVdx; res1 = simplify(F_x + k*x); % F + kx should be 0 assert(isAlways(res1 == 0, 'Unknown', 'false'), ... 'FAIL Step 1: F ≠ −kx'); fprintf('Step 1 PASS — F = −dV/dx = −kx\n'); %% ── Step 2: V₀ drops out of force ────────────────────────── % dV/dx should not contain V₀ assert(~any(symvar(dVdx) == V0), ... 'FAIL Step 2: V₀ appears in dV/dx'); fprintf('Step 2 PASS — V₀ drops out of force (constant offset irrelevant)\n'); %% ── Step 3: Newton's law → EOM ───────────────────────────── % mẍ = F = −kx → mẍ + kx = 0 % Verify: substituting F into mẍ = F gives mẍ + kx = 0 % EOM residual: m·d²x/dt² + k·x = 0 % For general solution x(t) = A cos(ωt + φ) with ω² = k/m: x_sol = A_amp * cos(omega * t + phi_0); d2x_sol = diff(x_sol, t, 2); % = −Aω² cos(ωt + φ) eom_residual = m * d2x_sol + k * x_sol; % = m(−Aω²cos) + k(Acos) = A cos(ωt+φ)(k − mω²) % Substitute ω² = k/m: eom_sub = subs(eom_residual, omega^2, k/m); res3 = simplify(eom_sub); assert(isAlways(res3 == 0, 'Unknown', 'false'), ... 'FAIL Step 3: x=Acos(ωt+φ) with ω²=k/m does not satisfy EOM'); fprintf('Step 3 PASS — x=Acos(ωt+φ), ω²=k/m satisfies mẍ+kx=0\n'); %% ── Step 4: ω² identification ────────────────────────────── % Divide EOM by m: ẍ + (k/m)x = 0 % Compare with ẍ + ω²x = 0 → ω² = k/m % Verify algebraically: k/m − ω² = 0 when ω² = k/m omega_sq = k / m; res4 = simplify(omega_sq - omega^2); res4_sub = subs(res4, omega^2, k/m); assert(isAlways(res4_sub == 0, 'Unknown', 'false'), ... 'FAIL Step 4: ω² ≠ k/m'); fprintf('Step 4 PASS — ω² = k/m\n'); %% ── Step 5: Period T = 2π/ω = 2π√(m/k) ──────────────────── T_period = 2*sym(pi) / omega; T_explicit = 2*sym(pi) * sqrt(m/k); % Substitute ω = √(k/m): T_sub = subs(T_period, omega, sqrt(k/m)); res5 = simplify(T_sub - T_explicit); assert(isAlways(res5 == 0, 'Unknown', 'false'), ... 'FAIL Step 5: T ≠ 2π√(m/k)'); fprintf('Step 5 PASS — T = 2π/ω = 2π√(m/k)\n'); %% ── Step 6: Energy conservation ───────────────────────────── % T = (1/2)mv², V = (1/2)kx² % E = T + V = (1/2)m(dx/dt)² + (1/2)kx² % For x = A cos(ωt+φ): v = −Aω sin(ωt+φ) % E = (1/2)mA²ω²sin² + (1/2)kA²cos² = (1/2)A²(mω²sin² + kcos²) % With mω² = k: E = (1/2)kA²(sin² + cos²) = (1/2)kA² v_sol = diff(x_sol, t); T_energy = m * v_sol^2 / 2; V_energy = k * x_sol^2 / 2; E_total = T_energy + V_energy; E_sub = subs(E_total, omega^2, k/m); E_expected = k * A_amp^2 / 2; res6 = simplify(E_sub - E_expected); assert(isAlways(res6 == 0, 'Unknown', 'false'), ... 'FAIL Step 6: E ≠ (1/2)kA²'); fprintf('Step 6 PASS — E = (1/2)kA² (energy conserved, time-independent)\n'); %% ── Step 7: Virial cross-check ⟨T⟩ = ⟨V⟩ ────────────────── % For SHO (n=2 potential): virial gives ⟨T⟩ = ⟨V⟩ = E/2 % ⟨T⟩ = (1/4)mA²ω² = (1/4)kA² = E/2 % ⟨V⟩ = (1/4)kA² = E/2 T_avg_sho = m * A_amp^2 * omega^2 / 4; V_avg_sho = k * A_amp^2 / 4; T_avg_sub = subs(T_avg_sho, omega^2, k/m); res7 = simplify(T_avg_sub - V_avg_sho); assert(isAlways(res7 == 0, 'Unknown', 'false'), ... 'FAIL Step 7: ⟨T⟩ ≠ ⟨V⟩'); fprintf('Step 7 PASS — ⟨T⟩ = ⟨V⟩ = (1/4)kA² (virial, consistent with F0026)\n'); %% ── Step 8: Numerical — 1 kg mass, k = 4 N/m ────────────── % ω = √(4/1) = 2 rad/s, T = 2π/2 = π ≈ 3.1416 s m_val = 1.0; k_val = 4.0; omega_val = sqrt(k_val / m_val); T_val = 2*pi / omega_val; assert(abs(omega_val - 2.0) < 1e-10, 'FAIL Step 8a: ω mismatch'); assert(abs(T_val - pi) < 1e-10, 'FAIL Step 8b: T mismatch'); fprintf('Step 8 PASS — m=1kg, k=4N/m → ω=%.1f rad/s, T=%.4f s\n', omega_val, T_val); %% ── Step 9: Unit consistency ──────────────────────────────── % [k/m] = (N/m)/kg = (kg/s²)/kg = 1/s² = [ω²] % Verify via dimensional proxy: if k = m·C for C [1/s²], then ω² = C syms C_dim real positive k_proxy = m * C_dim; omega_sq_proxy = k_proxy / m; res9 = simplify(omega_sq_proxy - C_dim); assert(isAlways(res9 == 0, 'Unknown', 'false'), ... 'FAIL Step 9: unit consistency'); fprintf('Step 9 PASS — [k/m] = [1/s²] = [ω²] (unit consistency)\n'); %% ── Step 10: Self-test — wrong frequency ──────────────────── % Wrong: ω² = k/(2m) (factor of 2 error) x_wrong = A_amp * cos(sqrt(k/(2*m)) * t + phi_0); d2x_wrong = diff(x_wrong, t, 2); eom_wrong = simplify(m * d2x_wrong + k * x_wrong); assert(~isAlways(eom_wrong == 0, 'Unknown', 'false'), ... 'FAIL Step 10a: wrong frequency not detected'); fprintf('Step 10a PASS — Wrong ω²=k/(2m) does not satisfy EOM\n'); % Quantify: residual = (k − k/2)·x_wrong = (k/2)·x_wrong eom_wrong_expected = k * x_wrong / 2; res10b = simplify(eom_wrong - eom_wrong_expected); assert(isAlways(res10b == 0, 'Unknown', 'false'), ... 'FAIL Step 10b: wrong residual ≠ (k/2)·x'); fprintf('Step 10b PASS — Wrong residual = (k/2)·x (quantified)\n'); %% ── Summary ───────────────────────────────────────────────── fprintf('\n✓ ALL 11 ASSERTIONS PASSED — F0027 audit complete.\n'); fprintf(' Chain: V=(1/2)kx² → F=−kx → mẍ+kx=0 → ω²=k/m\n'); fprintf(' Cross-ref: F0026 (virial ⟨T⟩=⟨V⟩ for n=2 potential)\n');