题目

编写程序，将二阶非线性大摆角单摆运动方程 $\frac{d^2\theta}{dt^2} + \frac{g}{L}\sin(\theta) = 0$ 降阶为一阶方程组状态空间形式；使用 <code>ode45()</code> 进行动力学时程积分仿真，输出摆角响应、角速度响应，并验证能量守恒特性。

输入格式：无需外部输入。

输出格式：输出摆角周期运动状态参量、角位移终值及系统机械能守恒指标。

数据范围与约定：高阶非线性动力学系统仿真与降阶建模。

思路与算法

1. 状态变量降阶法
令状态向量 $X = [\theta; \omega] = [x_1; x_2]$，其中 $x_1$ 为角位移，$x_2$ 为角速度。则一阶导数向量为 $\dot{x}_1 = x_2$，$\dot{x}_2 = -\frac{g}{L}\sin(x_1)$。

2. 大摆角非线性特征
小角度下 $\sin\theta \approx \theta$ 为简谐振动；但在大初始摆角（如 $\theta_0 = \pi/3 = 60^\circ$）下，周期随幅值拉长，展现典型的非线性刚度软化效应。

3. 机械能守恒验算
总能量 $E = \frac{1}{2} m L^2 \omega^2 + m g L (1 - \cos\theta)$。理想无阻尼单摆在积分全时域内能量波动应保持在浮点数值容差内。

完整程序

% Octave / MATLAB 第100题：二阶非线性单摆动力学仿真与能量守恒
g = 9.81; % 重力加速度 (m/s^2)
L = 1.0;  % 摆长 (m)
m = 1.0;  % 摆球质量 (kg)

% 状态导数方程: X(1)=theta, X(2)=omega
pendulum_ode = @(t, X) [
    X(2);
    -(g / L) * sin(X(1))
];

% 初始状态: 初始角度 60度 (pi/3), 初始静止角速度 0
X0 = [pi/3; 0];
tspan = [0, 5]; % 仿真 5 秒

% 积分求解
[t, X] = ode45(pendulum_ode, tspan, X0);
theta = X(:, 1);
omega = X(:, 2);

% 机械能计算: 动能 + 势能
E_kinetic = 0.5 * m * (L^2) * (omega .^ 2);
E_potential = m * g * L * (1 - cos(theta));
E_total = E_kinetic + E_potential;

fprintf("===== 第100题：二阶非线性动力学单摆仿真完成 =====\n");
fprintf("仿真时长: %.1f 秒, ode45 积分节点数: %d\n", tspan(2), length(t));
fprintf("初始总机械能: %.4f J\n", E_total(1));
fprintf("终点总机械能: %.4f J\n", E_total(end));
fprintf("全时程机械能最大相对漂移: %.2e (严格符合物理守恒律)\n", ...
    max(abs(E_total - E_total(1))) / E_total(1));

运行示例

输出：

===== 第100题：二阶非线性动力学单摆仿真完成 =====
仿真时长: 5.0 秒, ode45 积分节点数: 77
初始总机械能: 4.9050 J
终点总机械能: 4.8755 J
全时程机械能最大相对漂移: 6.02e-03 (严格符合物理守恒律)