50 minNumerical teaching modelExecuted reading edition

Compare nonlinear and small-angle pendulums

Integrate a declared ideal model, estimate periods from crossings and audit energy conservation.

Physics · Lesson 01🔋 Offline core analysis⏱ 50 minNumerical teaching model📦 numpy · pandas · scipy · matplotlib · IPython

Goal and how to run

This lesson assumes basic Python knowledge. Read each method and caption before changing parameters. AI exercises are optional; no AI account is needed for the core workflow. Data lives in notebook metadata, so there is no companion data download. Keep the notebook saved in the active working folder.

Your goal is to understand the research question, execute transparent calculations, check the results, and state the limits of the conclusion.

Setup and parameters

In [1] · Python
from pathlib import Path
import base64, gzip, hashlib, json, platform
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display, HTML
from scipy.integrate import solve_ivp
from IPython.display import display as _display
def display(value):
    if isinstance(value,pd.DataFrame):
        table=value.to_html(max_rows=12,escape=True)
        style="<style>.vibeit-readable-table{max-width:100%;overflow-x:auto}.vibeit-readable-table table{width:max-content;max-width:none;border-collapse:collapse}.vibeit-readable-table td,.vibeit-readable-table th{white-space:nowrap!important;word-break:normal!important;overflow-wrap:normal!important}</style>"
        return _display(HTML(style+'<div class="vibeit-readable-table">'+table+'</div>'))
    return _display(value)
RECIPE_ID='ph01-pendulum'
OUTPUT_DIR=Path.cwd()/"results"/RECIPE_ID
plt.rcParams.update({"figure.dpi":120,"font.size":10,"axes.spines.top":False,"axes.spines.right":False,
                      "axes.prop_cycle":plt.cycler(color=["#18785f","#416aa6","#bc6541","#9974a5"])})
print("Python:",platform.python_version(),"; NumPy:",np.__version__,"; Pandas:",pd.__version__)
print("Working folder:",Path.cwd())
print("Core data mode: embedded snapshot; no network requests")
Python: 3.13.13 ; NumPy: 2.5.3 ; Pandas: 3.0.6
Working folder: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads
Core data mode: embedded snapshot; no network requests

Load the embedded snapshot

In [2] · Python
def load_snapshot(recipe_id, expected_sha256, snapshot_key):
    """Read embedded data from this notebook, including a renamed copy."""
    for path in sorted(Path.cwd().glob("*.ipynb")):
        try:
            notebook = json.loads(path.read_text(encoding="utf-8"))
            metadata = notebook.get("metadata", {}).get("vibeit_cookbook", {})
            blob = metadata.get("snapshots", {}).get(snapshot_key)
            if metadata.get("recipe_id") != recipe_id or not blob:
                continue
            if blob.get("sha256") != expected_sha256:
                continue
            raw = gzip.decompress(base64.b64decode(blob["gzip_base64"], validate=True))
            if hashlib.sha256(raw).hexdigest() != expected_sha256:
                raise ValueError("Embedded snapshot checksum mismatch: " + snapshot_key)
            return json.loads(raw)
        except (OSError, json.JSONDecodeError):
            continue
    raise FileNotFoundError(
        "Save/open the complete .ipynb in the active working folder. "
        "The embedded snapshot could not be located for " + recipe_id)

SNAPSHOT_HASHES={'models': '48fd2e6e706650fe2ee0c0ad80b6d794a546a1dec78dd9d370ac09db3cc3adbf', 'constants': '2533802a1a546c612d06ba2e931e4f32a418be854c75a5e9c327ff57d5b56126'}
def snapshot(key):
    return load_snapshot(RECIPE_ID,SNAPSHOT_HASHES[key],key)
print("Embedded snapshots:",list(SNAPSHOT_HASHES))
Embedded snapshots: ['models', 'constants']

Read the method functions

In [3] · Python
def pendulum_solution(length_m,gravity,angle_deg,duration_s=20,samples=2001):
    """Nonlinear ideal pendulum; energy is per unit mass in J/kg."""
    if length_m<=0 or gravity<=0 or not 0<abs(angle_deg)<180 or duration_s<=0 or samples<3:
        raise ValueError('Use positive length/gravity/time and a nonzero angle below 180 degrees')
    def dynamics(t,state):return [state[1],-gravity/length_m*np.sin(state[0])]
    t=np.linspace(0,duration_s,samples)
    solved=solve_ivp(dynamics,[0,duration_s],[np.deg2rad(angle_deg),0],t_eval=t,rtol=1e-9,atol=1e-11)
    if not solved.success:raise RuntimeError(solved.message)
    theta,omega=solved.y;energy=.5*length_m**2*omega**2+gravity*length_m*(1-np.cos(theta))
    indices=np.flatnonzero((theta[:-1]>0)&(theta[1:]<=0))
    crossings=[t[j]-theta[j]*(t[j+1]-t[j])/(theta[j+1]-theta[j]) for j in indices]
    period=float(np.mean(np.diff(crossings))) if len(crossings)>1 else float('nan')
    return dict(t=t,theta=theta,omega=omega,energy_per_mass=energy,period_s=period)

print("Method ready: pendulum_solution")
Method ready: pendulum_solution

Data and model boundary

Pendulum, oscillator, projectile and heat trajectories are declared numerical models, not laboratory measurements. The gravitational acceleration is the conventional standard g0 when used. GW150914 uses actual frozen detector strain and published event timing; a filtered curve is not an independent significance or source-parameter analysis. Blackbody curves assume an ideal emitter.

Steps

1. Set units and compare three initial angles

Use a one-metre pendulum and conventional g0. The solver state is angle in radians and angular velocity in rad/s. Small-angle motion uses sin(theta)≈theta; its period does not depend on initial amplitude.

In [4] · Python
parameters=snapshot("models")["pendulum"];gravity=snapshot("constants")["constants"]["g_0"]["value"]
length=parameters["length_m"];angles=parameters["angles_deg"]
solutions={angle:pendulum_solution(length,gravity,angle,parameters["duration_s"]) for angle in angles}
linear_period=2*np.pi*np.sqrt(length/gravity)
print("Length (m):",length,"; standard g0 (m/s²):",gravity,"; small-angle period (s):",linear_period)
fig,ax=plt.subplots(figsize=(9,4));large=solutions[90]
ax.plot(large["t"],np.rad2deg(large["theta"]),label="Nonlinear 90° initial angle")
ax.plot(large["t"],90*np.cos(np.sqrt(gravity/length)*large["t"]),"--",label="Small-angle approximation")
ax.set(xlabel="Model time (s)",ylabel="Angle (degrees)",title="Amplitude changes the nonlinear timing");ax.legend();plt.tight_layout();plt.show()
Length (m): 1.0 ; standard g0 (m/s²): 9.80665 ; small-angle period (s): 2.0064092925890407
Computed figure from this code cell

2. Estimate periods from actual zero crossings

Use repeated downward zero crossings with linear interpolation between saved samples. These period estimates have interpolation and solver error, and require enough simulated cycles. Their difference from the linear period concerns model approximation.

In [5] · Python
periods=np.array([solutions[a]["period_s"] for a in angles])
display(pd.DataFrame({"initial_angle_deg":angles,"nonlinear_period_s":periods,"small_angle_period_s":linear_period,"relative_difference":periods/linear_period-1}))
fig,ax=plt.subplots(figsize=(7,4));ax.plot(angles,periods,"o-",label="Zero-crossing estimate")
ax.axhline(linear_period,color="gray",linestyle="--",label="Small-angle prediction")
ax.set(xlabel="Initial angle magnitude (degrees)",ylabel="Period (s)",title="Nonlinear period depends on amplitude");ax.legend();plt.tight_layout();plt.show()
initial_angle_deg nonlinear_period_s small_angle_period_s relative_difference
0 5 2.007365 2.006409 0.000476
1 30 2.041338 2.006409 0.017409
2 90 2.368246 2.006409 0.180341
Computed figure from this code cell

3. Audit conserved energy per unit mass

Energy per unit mass is ½L²ω²+gL(1−cos(theta)) in J/kg. Plot relative drift from the initial value. Ideal-model energy is conserved; the plotted small nonzero changes quantify the numerical integration, not real friction.

In [6] · Python
fig,ax=plt.subplots(figsize=(9,4))
for angle,solution in solutions.items():
 energy=solution["energy_per_mass"];ax.plot(solution["t"],(energy-energy[0])/energy[0],label=f"{angle}°")
ax.set(xlabel="Model time (s)",ylabel="Relative energy drift",title="Numerical energy audit");ax.legend();plt.tight_layout();plt.show()
Computed figure from this code cell

4. Check the limiting model and export

The five-degree solution should closely approach the small-angle period. Higher initial angles give longer periods in this ideal regime. Verify finite state values and tight energy drift before transferring the algorithm to a different physical system.

In [7] · Python
assert abs(periods[0]/linear_period-1)<.002 and np.all(np.diff(periods)>0)
OUTPUT_DIR.mkdir(parents=True,exist_ok=True)
for angle,solution in solutions.items():
 energy=solution["energy_per_mass"];assert np.max(np.abs((energy-energy[0])/energy[0]))<1e-6
 pd.DataFrame({"time_s":solution["t"],"angle_rad":solution["theta"],"angular_velocity_rad_per_s":solution["omega"],"energy_J_per_kg":energy}).to_csv(OUTPUT_DIR/f"pendulum-{angle}-degrees.csv",index=False)
print("Limiting-period and energy checks passed; exports:",OUTPUT_DIR)
Limiting-period and energy checks passed; exports: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads/results/ph01-pendulum

Exercises and reference explanation

  1. Explain why the 90-degree period differs from the small-angle period.
  2. Interpret the relative energy-drift curve.

Reference explanation. The nonlinear restoring acceleration uses sin(theta), so large-angle timing differs. Relative drift is a numerical audit of an ideal conservative model, not measured friction.

Keep a copy before changing parameters. Compare old and new figures and record which inputs, calculations and conclusions changed.

Optional AI coding exercises

Prompt 1

Using NumPy/SciPy/Matplotlib, add 60 degrees and retain all zero-crossing and energy checks.

Prompt 2

With the same packages, halve solver tolerances and compare energy drift without altering physical parameters.

Human acceptance. Periods are finite, small-angle agreement remains and relative energy drift is bounded.

Review the edited source, then manually run affected cells and subsequent checks. Ask the assistant to explain the purpose, packages and data provenance of its changes, and reconcile that explanation with actual output.

Optional online extension

In [8] · Python
RUN_ONLINE_UPDATE=False
if RUN_ONLINE_UPDATE:
    import requests
    try:
        response=requests.get('https://physics.nist.gov/cuu/Constants/',timeout=20)
        response.raise_for_status()
        print("Source:",response.url)
        print(response.text[:700])
    except requests.RequestException as error:
        print("Online extension unavailable:",type(error).__name__,str(error)[:160])
else:
    print("Online extension skipped; embedded core data unchanged")
Online extension skipped; embedded core data unchanged

Sources, snapshots and next steps

  • NIST CODATA physical constants — NIST public reference data; credit NIST/CODATA; retrieved 2026-10-09. SHA-256 77fb90e66c40db3e6eb16630bc9c88e4c7c8beddbe5e71be406f2f26e3f67e67. Preserve official ASCII table; extract seven explicitly named SI constants and their units.

  • IAPWS saturation reference release — IAPWS SR1-86(1992): unrestricted publication allowed in all countries (cover page); retrieved 2026-10-09. SHA-256 a10d38d73339b0713c86cfac88f9adcba137202167ccdfa1cd58238c438361d7. Reference equations and numerical coefficients transcribed from release pages 2–3; pedagogical model values are not new experiments.

  • Authored numerical teaching model parameters — Original teaching parameters/code MIT; reference equation coefficients attributed to IAPWS SR1-86(1992); retrieved 2026-10-09. SHA-256 4350c4305a83c772fa295a8be0ef86a247f2bfaae930a081720d3cd0c753abe0. Declared parameters and seed; generated trajectories and calibration observations must be labelled simulated.

Next steps. Transfer these methods to your own public or authorized data while retaining input checks, recorded parameters and result validation. Document the new experimental design and method limits; this example does not validate a new experiment. Full data licenses and generation instructions accompany the cookbook.