MATLAB - 仿真单摆的周期性摆动

系列文章目录

前言

本例演示如何使用 Symbolic Math Toolbox™ 模拟单摆的运动。推导摆的运动方程，然后对小角度进行分析求解，对任意角度进行数值求解。

一、步骤 1：推导运动方程

摆是一个遵循微分方程的简单机械系统。摆最初静止在垂直位置。当摆移动一个角度 θ 并释放时，重力将其拉回静止位置。它的动量会使它过冲并到达 -θ 角（如果没有摩擦力），以此类推。由于重力的作用，钟摆运动的恢复力为 -mgsinθ。因此，根据牛顿第二定律，质量乘以加速度必须等于 -mgsinθ。

syms m a g theta(t)
eqn = m*a == -m*g*sin(theta)

eqn(t) = a m=−g m sin(θ(t))

对于长度为 r 的摆锤，摆锤的加速度等于角加速度乘以 r。

$a=r\frac{d^2\theta}{dt^2}.$

用子项代替 a。

syms r
eqn = subs(eqn,a,r*diff(theta,2))

$\begin{aligned}eqn(t)&=\\mr\frac{\partial^2}{\partial t^2}&\theta(t)=-gm\sin(\theta(t))\end{aligned}$

使用 isolate 隔离公式中的角加速度。

eqn = isolate(eqn,diff(theta,2))

$\begin{aligned}\text{eqn =} \frac{\partial^2}{\partial t^2}\left.\theta(t)=-\frac{g\sin(\theta(t))}r\right.\end{aligned}$

将常数 g 和 r 合并为一个参数，也称为固有频率。

$\omega_0=\sqrt{\dfrac{g}{r}}.$

syms omega_0
eqn = subs(eqn,g/r,omega_0^2)

$\begin{aligned}\text{eqn =}\frac{\partial^2}{\partial t^2}\theta(t)=-\omega_0^2\sin(\theta(t))\end{aligned}$

二、步骤 2：运动方程线性化

运动方程是非线性的，因此难以用解析法求解。假设角度很小，利用 sinθ 的泰勒展开式将方程线性化。

syms x
approx = taylor(sin(x),x,'Order',2);
approx = subs(approx,x,theta(t))

$approx = \theta(t)$

运动方程变成了线性方程。

eqnLinear = subs(eqn,sin(theta(t)),approx)

$\begin{aligned}\text{eqnLinear}&=\\\frac{\partial^2}{\partial t^2}\theta(t)&=-\omega_0^2\theta(t)\end{aligned}$

三、步骤 3：分析求解运动方程

使用 dsolve 求解方程 eqnLinear。将初始条件指定为第二个参数。使用 assume 假设 ω0 为实数，简化解法。

syms theta_0 theta_t0
theta_t = diff(theta);
cond = [theta(0) == theta_0, theta_t(0) == theta_t0];
assume(omega_0,'real')
thetaSol(t) = dsolve(eqnLinear,cond)

$\begin{aligned}\text{thetasol(t)}&=\theta_0\cos(\omega_0t)+\frac{\theta_{t0}\sin(\omega_0t)}{\omega_0}\end{aligned}$

四、步骤 4：ω0 的物理意义

项 ω 0 t 称为相位。余弦函数和正弦函数每 2π 重复一次。改变 ω 0 t 变化 2π 所需的时间称为时间周期。

$T=\dfrac{2\pi}{\omega0}=2\pi\sqrt{\dfrac{r}{g}}.$

时间周期 T 与摆长的平方根成正比，与质量无关。对于线性运动方程，时间周期与初始条件无关。

五、步骤 5：绘制摆的运动图

绘制小角度近似的摆运动图。

定义物理参数：

重力加速度 $g=9.81\mathrm{m/s}^{2}$
摆长 $\text{r=1 m}$

gValue = 9.81;
rValue = 1;
omega_0Value = sqrt(gValue/rValue);
T = 2*pi/omega_0Value;

设置初始条件。

theta_0Value  = 0.1*pi; % Solution only valid for small angles.
theta_t0Value = 0;      % Initially at rest.

将物理参数和初始条件代入一般解法。

vars   = [omega_0      theta_0      theta_t0];
values = [omega_0Value theta_0Value theta_t0Value];
thetaSolPlot = subs(thetaSol,vars,values);

绘制谐摆运动图。

fplot(thetaSolPlot(t*T)/pi, [0 5]);
grid on;
title('Harmonic Pendulum Motion');
xlabel('t/T');
ylabel('\theta/\pi');

求出 θ(t) 的解后，想象一下摆的运动。

x_pos = sin(thetaSolPlot);
y_pos = -cos(thetaSolPlot);
fanimator(@fplot,x_pos,y_pos,'ko','MarkerFaceColor','k','AnimationRange',[0 5*T]);
hold on;
fanimator(@(t) plot([0 x_pos(t)],[0 y_pos(t)],'k-'),'AnimationRange',[0 5*T]);
fanimator(@(t) text(-0.3,0.3,"Timer: "+num2str(t,2)+" s"),'AnimationRange',[0 5*T]);

输入 playAnimation 命令播放钟摆运动的动画。

六、步骤 6：使用恒定能量路径确定非线性摆运动

为了理解摆的非线性运动，请使用总能量方程来直观显示摆的运动轨迹。总能量是守恒的。

$E=\dfrac12mr^2{\left(\dfrac{d\theta}{dt}\right)}^2+mgr(1-\cos\theta)$

使用三角函数特性 $1-\cos\theta=2\sin^2(\theta/2)$ 和关系式 $\omega_0=\sqrt{g/r}$ 重写比例能量。

$\dfrac E{mr^2}=\dfrac12\left[\left(\dfrac{d\theta}{dt}\right)^2+\left(2\omega_0\sin\dfrac\theta2\right)^2\right]$

由于能量守恒，摆的运动可以用相空间中的恒定能量路径来描述。相空间是一个抽象空间，坐标为 θ 和 dθ/dt。使用 fcontour 将这些路径可视化。

syms theta theta_t omega_0
E(theta, theta_t, omega_0) = (1/2)*(theta_t^2+(2*omega_0*sin(theta/2))^2);
Eplot(theta, theta_t) = subs(E,omega_0,omega_0Value);

figure;
fc = fcontour(Eplot(pi*theta, 2*omega_0Value*theta_t), 2*[-1 1 -1 1], ...
    'LineWidth', 2, 'LevelList', 0:5:50, 'MeshDensity', 1+2^8);
grid on;
title('Constant Energy Contours in Phase Space ( \theta vs. \theta_t )');
xlabel('\theta/\pi');
ylabel('\theta_t/2\omega_0');

恒定能量等值线围绕 θ 轴和 dθ/dt 轴对称，沿 θ 轴呈周期性分布。图中显示了两个行为截然不同的区域。

等值线图的较低能量相互靠近。摆锤在两个最大角度和速度之间来回摆动。摆锤的动能不足以克服重力能，使摆锤绕一圈。

等值线图中的高能量不会自行闭合。摆锤始终沿着一个角度方向运动。钟摆的动能足以克服重力能，使钟摆能够绕一圈。

七、步骤 7：求解非线性运动方程

非线性运动方程是二阶微分方程。使用 ode45 求解器对这些方程进行数值求解。由于 ode45 只接受一阶系统，因此请将系统简化为一阶系统。然后生成函数句柄，作为 ode45 的输入。

将二阶 ODE 重写为一阶 ODE 系统。

syms theta(t) theta_t(t) omega_0
eqs = [diff(theta)   == theta_t;
       diff(theta_t) == -omega_0^2*sin(theta)]

$\left.\text{eqs}(t)=\\\left(\begin{aligned}\frac{\partial}{\partial t}&\theta(t)=\theta_t(t)\\\\\frac{\partial}{\partial t}&\theta_t(t)=-\omega_0^2\sin(\theta(t))\end{aligned}\right.\right)$

eqs  = subs(eqs,omega_0,omega_0Value);
vars = [theta, theta_t];

求出系统的质量矩阵 M 和方程 F 的右边。

[M,F] = massMatrixForm(eqs,vars)

$\begin{array}{rcl}\text{M}&=&\\&&\begin{bmatrix}1&0\\0&1\end{bmatrix}\end{array}$

$\begin{aligned}{F}&=\\&\begin{bmatrix}\theta_t(t)\\-\dfrac{981\sin(\theta(t))}{100}\end{bmatrix}\end{aligned}$

M 和 F 指的就是这种形式。

$M(t,x(t))\dfrac{dx}{dt}=F(t,x(t)).$

为简化进一步计算，可将系统改写为 $dx/dt=f(t,x(t)).$ 的形式。

f = M\F

$\begin{aligned}\text{f}&=\\&\begin{bmatrix}\theta_t(t)\\-\dfrac{981\sin(\theta(t))}{100}\end{bmatrix}\end{aligned}$

使用 odeFunction 将 f 转换为 MATLAB 函数句柄。生成的函数句柄是 MATLAB ODE 求解器 ode45 的输入。

f = odeFunction(f, vars)

f = function_handle with value:
    @(t,in2)[in2(2,:);sin(in2(1,:)).*(-9.81e+2./1.0e+2)]

八、步骤 8：求解封闭能量等值线的运动方程

使用 ode45 求解封闭能量等值线的 ODE。

从能量等值线图来看，封闭等值线满足条件 $\theta_0=0,\theta_{t0}/2\omega_0\leq1$ . 将 θ 和 dθ/dt 的初始条件存储在变量 x0 中。

x0 = [0; 1.99*omega_0Value];

指定一个从 0 秒到 10 秒的时间间隔，用于求解。这个时间间隔允许摆锤经历两个完整的周期。

tInit  = 0;
tFinal = 10;

求解 ODE。

sols = ode45(f,[tInit tFinal],x0)

sols = struct with fields:
     solver: 'ode45'
    extdata: [1x1 struct]
          x: [0 3.2241e-05 1.9344e-04 9.9946e-04 0.0050 0.0252 0.1259 0.3449 0.6020 0.8591 1.1161 1.3597 1.5996 1.8995 2.2274 2.4651 2.7028 2.9567 3.2138 3.4709 3.7150 3.9511 4.2483 4.5759 4.8239 5.0719 5.3182 5.5764 5.8346 6.0803 ... ] (1x45 double)
          y: [2x45 double]
      stats: [1x1 struct]
      idata: [1x1 struct]

sols.y(1,:) 表示角位移 θ，sols.y(2,:) 表示角速度 dθ/dt。

绘制闭合路径解。

figure;

yyaxis left;
plot(sols.x, sols.y(1,:), '-o');
ylabel('\theta (rad)');

yyaxis right;
plot(sols.x, sols.y(2,:), '-o');
ylabel('\theta_t (rad/s)');

grid on;
title('Closed Path in Phase Space');
xlabel('t (s)');

可视化钟摆的运动。

x_pos = @(t) sin(deval(sols,t,1));
y_pos = @(t) -cos(deval(sols,t,1));
figure;
fanimator(@(t) plot(x_pos(t),y_pos(t),'ko','MarkerFaceColor','k'));
hold on;
fanimator(@(t) plot([0 x_pos(t)],[0 y_pos(t)],'k-'));
fanimator(@(t) text(-0.3,1.5,"Timer: "+num2str(t,2)+" s"));

输入 playAnimation 命令播放钟摆运动的动画。

九、步骤 9：开放式能量等值线的求解

使用 ode45 求解开放式能量等值线的 ODE。从能量等值线图来看，开放式等值线满足条件 $\theta_0=0\text{,}\theta_{t0}/2\omega_0>1.$ .

x0 = [0; 2.01*omega_0Value];
sols = ode45(f, [tInit, tFinal], x0);

绘制开放式能量等值线的解。

figure;

yyaxis left;
plot(sols.x, sols.y(1,:), '-o');
ylabel('\theta (rad)');

yyaxis right;
plot(sols.x, sols.y(2,:), '-o');
ylabel('\theta_t (rad/s)');

grid on;
title('Open Path in Phase Space');
xlabel('t (s)');

可视化钟摆的运动。

x_pos = @(t) sin(deval(sols,t,1));
y_pos = @(t) -cos(deval(sols,t,1));
figure;
fanimator(@(t) plot(x_pos(t),y_pos(t),'ko','MarkerFaceColor','k'));
hold on;
fanimator(@(t) plot([0 x_pos(t)],[0 y_pos(t)],'k-'));
fanimator(@(t) text(-0.3,1.5,"Timer: "+num2str(t,2)+" s"));