个人网站建设主要功能,南京自助建站网站,wordpress小工具安装,公众号平台制作系列文章目录 前言
本例演示如何使用 Symbolic Math Toolbox™ 模拟单摆的运动。推导摆的运动方程#xff0c;然后对小角度进行分析求解#xff0c;对任意角度进行数值求解。 一、步骤 1#xff1a;推导运动方程
摆是一个遵循微分方程的简单机械系统。摆最初静止在垂直位置…系列文章目录 前言
本例演示如何使用 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。
syms r
eqn subs(eqn,a,r*diff(theta,2)) 使用 isolate 隔离公式中的角加速度。
eqn isolate(eqn,diff(theta,2)) 将常数 g 和 r 合并为一个参数也称为固有频率。 syms omega_0
eqn subs(eqn,g/r,omega_0^2) 二、步骤 2运动方程线性化
运动方程是非线性的因此难以用解析法求解。假设角度很小利用 sinθ 的泰勒展开式将方程线性化。
syms x
approx taylor(sin(x),x,Order,2);
approx subs(approx,x,theta(t)) 运动方程变成了线性方程。
eqnLinear subs(eqn,sin(theta(t)),approx) 三、步骤 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) 四、步骤 4ω0 的物理意义
项 ω 0 t 称为相位。余弦函数和正弦函数每 2π 重复一次。改变 ω 0 t 变化 2π 所需的时间称为时间周期。 时间周期 T 与摆长的平方根成正比与质量无关。对于线性运动方程时间周期与初始条件无关。
五、步骤 5绘制摆的运动图
绘制小角度近似的摆运动图。
定义物理参数
重力加速度 摆长 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使用恒定能量路径确定非线性摆运动
为了理解摆的非线性运动请使用总能量方程来直观显示摆的运动轨迹。总能量是守恒的。 使用三角函数特性 和关系式 重写比例能量。 由于能量守恒摆的运动可以用相空间中的恒定能量路径来描述。相空间是一个抽象空间坐标为 θ 和 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, 12^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)] eqs subs(eqs,omega_0,omega_0Value);
vars [theta, theta_t]; 求出系统的质量矩阵 M 和方程 F 的右边。
[M,F] massMatrixForm(eqs,vars) M 和 F 指的就是这种形式。 为简化进一步计算可将系统改写为 的形式。
f M\F 使用 odeFunction 将 f 转换为 MATLAB 函数句柄。生成的函数句柄是 MATLAB ODE 求解器 ode45 的输入。
f odeFunction(f, vars)
f function_handle with value:(t,in2)[in2(2,:);sin(in2(1,:)).*(-9.81e2./1.0e2)]
八、步骤 8求解封闭能量等值线的运动方程
使用 ode45 求解封闭能量等值线的 ODE。
从能量等值线图来看封闭等值线满足条件 . 将 θ 和 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: ode45extdata: [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。从能量等值线图来看开放式等值线满足条件 .
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)); 输入 playAnimation 命令播放钟摆运动的动画。
