钟摆模型主要物理参量如下图所示:
syms m a g theta(t)
eqn = m*a == -m*g*sin(theta)
对于有长度的摆锤,摆锤的加速度:
syms r
eqn = subs(eqn,a,r*diff(theta,2))
使用隔离,隔离方程eqn中的角加速度。得到:
eqn = isolate(eqn,diff(theta,2))
将常数和频率收集到一个参数中,该参数也称为固有频率:
syms omega_0
eqn = subs(eqn,g/r,omega_0^2)
运动方程是非线性的,所以很难解析求解。假设角度很小,用泰勒展开式将方程线性化:
syms x
approx = taylor(sin(x),x,'Order',2);
approx = subs(approx,x,theta(t))
eqnLinear = subs(eqn,sin(theta(t)),approx)
用dsolve解方程组。将初始条件指定为第二个参数。通过使用假设来简化解决方案:
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)
钟摆周期:
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');
钟摆动画:
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]);
如若转载,请注明出处:https://www.dasum.com/28549.html