单摆相空间与能量 notebook

这是手机端静态预览：保留说明和代码；图表请在电脑端 Jupyter 中运行 notebook 查看。

????????????

Audience:

• ?????????????????????????-???????????

Prerequisites:

• ???????????
• ????????? L = T - V?
• ????????? theta???? omega?

Learning goals:

• ????????????
• ??? Runge-Kutta ?????????
• ?????????????
• ????????????????

Outline

1. ???????????
2. ??????????
3. ??????????
4. ????????????
5. ?????????????????????

from __future__ import annotations

from dataclasses import dataclass

import matplotlib.pyplot as plt
import numpy as np

plt.rcParams["figure.figsize"] = (9, 4.8)
plt.rcParams["axes.grid"] = True
plt.rcParams["grid.alpha"] = 0.25

1. ??

?????????

??-????????

???????????????theta_dot = omega?omega_dot = -(g/l)sin(theta)?

@dataclass
class PendulumParams:
    length: float = 1.2
    gravity: float = 9.8
    mass: float = 1.0
    damping: float = 0.0


def derivatives(y: np.ndarray, params: PendulumParams) -> np.ndarray:
    theta, omega = y
    theta_dot = omega
    omega_dot = -(params.gravity / params.length) * np.sin(theta) - params.damping * omega
    return np.array([theta_dot, omega_dot], dtype=float)


def rk4_step(y: np.ndarray, dt: float, params: PendulumParams) -> np.ndarray:
    k1 = derivatives(y, params)
    k2 = derivatives(y + 0.5 * dt * k1, params)
    k3 = derivatives(y + 0.5 * dt * k2, params)
    k4 = derivatives(y + dt * k3, params)
    return y + (dt / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)


def energy(theta: np.ndarray, omega: np.ndarray, params: PendulumParams) -> np.ndarray:
    kinetic = 0.5 * params.mass * params.length**2 * omega**2
    potential = params.mass * params.gravity * params.length * (1 - np.cos(theta))
    return kinetic + potential


def lagrangian(theta: np.ndarray, omega: np.ndarray, params: PendulumParams) -> np.ndarray:
    kinetic = 0.5 * params.mass * params.length**2 * omega**2
    potential = params.mass * params.gravity * params.length * (1 - np.cos(theta))
    return kinetic - potential


def simulate(theta0: float, omega0: float, params: PendulumParams, t_max: float = 12.0, dt: float = 0.002) -> dict[str, np.ndarray]:
    steps = int(t_max / dt) + 1
    t = np.linspace(0.0, t_max, steps)
    y = np.zeros((steps, 2), dtype=float)
    y[0] = [theta0, omega0]

    for i in range(1, steps):
        y[i] = rk4_step(y[i - 1], dt, params)

    theta = y[:, 0]
    omega = y[:, 1]
    return {
        "t": t,
        "theta": theta,
        "omega": omega,
        "energy": energy(theta, omega, params),
        "lagrangian": lagrangian(theta, omega, params),
        "action": np.cumsum(lagrangian(theta, omega, params)) * dt,
    }

2. ?????????

???????????????????? 55???????? sin(theta)?

params = PendulumParams(length=1.2, gravity=9.8, damping=0.0)
result = simulate(theta0=np.deg2rad(55), omega0=0.0, params=params, t_max=12.0, dt=0.002)

summary = {
    "theta_min_rad": float(result["theta"].min()),
    "theta_max_rad": float(result["theta"].max()),
    "omega_max_abs_rad_s": float(np.abs(result["omega"]).max()),
    "energy_initial_J": float(result["energy"][0]),
    "energy_final_J": float(result["energy"][-1]),
    "action_final": float(result["action"][-1]),
}
summary

3. ???????

????????????????????????????????????????

fig, axes = plt.subplots(1, 2, figsize=(11, 4.4))

axes[0].plot(result["t"], result["theta"], color="#087c75", lw=2)
axes[0].set_title("Angle over time")
axes[0].set_xlabel("t (s)")
axes[0].set_ylabel("theta (rad)")

axes[1].plot(result["theta"], result["omega"], color="#c35a2b", lw=2)
axes[1].set_title("Phase portrait")
axes[1].set_xlabel("theta (rad)")
axes[1].set_ylabel("omega (rad/s)")

fig.tight_layout()
plt.show()

4. ??????

???????????????????????????????????????????????????

energy_drift = result["energy"] - result["energy"][0]
relative_drift = energy_drift / result["energy"][0]

fig, ax = plt.subplots(figsize=(9, 4))
ax.plot(result["t"], relative_drift, color="#7458b6", lw=2)
ax.set_title("Relative energy drift")
ax.set_xlabel("t (s)")
ax.set_ylabel("(E(t) - E(0)) / E(0)")
plt.show()

float(np.max(np.abs(relative_drift)))

5. ????

????????????????????? beta * omega ???????????????????????

damped_params = PendulumParams(length=1.2, gravity=9.8, damping=0.18)
damped = simulate(theta0=np.deg2rad(55), omega0=0.0, params=damped_params, t_max=16.0, dt=0.002)

fig, axes = plt.subplots(1, 2, figsize=(11, 4.4))

axes[0].plot(damped["t"], damped["energy"], color="#087c75", lw=2)
axes[0].set_title("Damped energy")
axes[0].set_xlabel("t (s)")
axes[0].set_ylabel("E (J)")

axes[1].plot(damped["theta"], damped["omega"], color="#c35a2b", lw=2)
axes[1].set_title("Damped phase portrait")
axes[1].set_xlabel("theta (rad)")
axes[1].set_ylabel("omega (rad/s)")

fig.tight_layout()
plt.show()

??

1. ?????? 5??????????????? 2*pi*sqrt(l/g)?
2. ? dt ? 0.002 ?? 0.02????????????
3. ? damping ?? 0.5?????????????????

def estimate_period_from_zero_crossings(t: np.ndarray, theta: np.ndarray) -> float:
    crossings = []
    for i in range(1, len(theta)):
        if theta[i - 1] > 0 and theta[i] <= 0:
            crossings.append(t[i])
    if len(crossings) < 2:
        return float("nan")
    return 2 * float(np.mean(np.diff(crossings)))

small_angle = simulate(
    theta0=np.deg2rad(5),
    omega0=0.0,
    params=PendulumParams(length=1.2, gravity=9.8, damping=0.0),
    t_max=12.0,
    dt=0.002,
)

numerical_period = estimate_period_from_zero_crossings(small_angle["t"], small_angle["theta"])
small_angle_period = 2 * np.pi * np.sqrt(params.length / params.gravity)

{
    "numerical_period_s": numerical_period,
    "small_angle_period_s": float(small_angle_period),
    "absolute_error_s": abs(numerical_period - small_angle_period),
}

????

• ???????? sin?Python ??????????????
• ???? dt ??????????????????????????????????????
• ????? S ???????????????????????????????????