单摆振动运规律的傅立叶分析

单摆振动运动规律的傅立叶分析
院系 XX 专业XX 姓名 XX
[问题]单摆振动运动规律的傅立叶分析 对单摆振动规律进行傅立叶分析。

[数学模型]
如A5.2图所示，设摆锤质量为m ，角位置为θ，摆锤的运动方程为
22d sin d ml mg t
θ
θ=-， 即 22d sin d g
t l
θθ=-， (5.4.1) 在小角度的情况下，sin θ ≈ θ，可得
22
02
d 0d t
θωθ+=， (5.4.2)
其中0ω=
，ω0为圆频率。

可知：单摆在小角度时作简谐振动，
小角度周期为
00
2π
2T ω=
=。

(5.4.3) 可见：在小角振动的情况下，单摆的周期与角振幅无关，这称为单摆的等时性。

摆锤的角速度为ω = d θ/d t ，因此
22d d d d d d d d d d t t t θωωθω
ωθθ
===， 由(5.4.1)式可得
d sin d g
l
ωωθθ=-，
积分得
21cos 2g
C l
ωθ=+， 当t = 0时，ω = 0，θ = θm ，可得C = -g cos θm /l 。

因此角速度大小为
d d t θω==。

(5.4.4) 注意：角速度是单位时间内角度的变化率d θ/d t ，圆频率是简谐运动中2π时间内周期性运动
m g
A5.2图
的次数2π/T，它们常用字母ω表示，单位也相同，但意义不同。

单摆的周期为
m m
00
π
T T
==⎰。

(5.4.5)
对于任何角振幅θm，通过数值积分和符号积都能计算周期。

利用半角公式可得
m
1
π
T T
θ
=⎰
设
m
sin
2
k
θ
=，(5.4.6) 并设k sin x = sin(θ/2)，因此
1
cos d cos d
22
k x x
θ
θ
=，可得
π/2π/2
00
00
12
ππ
T T T
==
⎰⎰，
即
π/2
2
π
T T
=⎰(5.4.7) 这是椭圆积分。

第一类完全椭圆积分定义为
π/2
K()k=⎰(5.4.8) 周期为
2
K()
π
T T k
=。

(5.4.9)
[算法]对于任何一个角振幅θm，利用(5.4.5)式，通过MA TLAB数值积分指令quadl和符号积分指令int都可计算单摆的周期。

利用MATLAB的完全椭圆积分指令ellipke也能计算单摆的周期。

不过，MA TLAB定义的一类完全椭圆积分定义为
π/2
K()m=⎰，(5.4.8*) 周期可表示为
2
K()
π
T T m
=。

(5.4.9*) 其中
2m
sin
2
m
θ
=。

(5.4.6*)
取t * = ω0t 为量纲时间，(5.4.1)式可化为
2*2
d sin d t θ
θ=-， (5.4.1*) 取θ(1) = θ，θ(2) = d θ/d t *，可得
*d (1)(2)d t θθ=，*
d (2)
sin (1)d t θθ=-。

(5.4.1**) 初始条件是θ(1) = θm ，θ(2) = 0。

由于
*0
d (1)d (2)d d t t θθω
θωω=
==， 因此θ(2)表示以ω0为单位的角速度。

由于
00002π
t t
t T ωωω=
=， 因此，无量纲时间除以2π就是以小角周期T 0为单位的时间。

取l 为坐标单位，则摆锤的坐标可表示为
x * = x /l = sin θ，y * = y /l = -cos θ。

(5.4.9*)
简谐振动的角位移可表示为
m m 0*0012π1
cos(
)cos()/t t T T T T
θθθω==。

(5.4.10*)
其中T * = T /T 0是约化周期。

简谐振动的角速度可表示为
m m
00**
00002π112π1sin()sin()//t t T T T T T T T T
θωθωω=-=-。

(5.4.11*) [程序]zhou13_1ode.m 如下。

%单摆振动规律的傅立叶分析
clear %清除变量 thetam=input('请输入单摆角振幅的度数:');%键盘输入角振幅 thm=thetam*pi/180; %化为弧度
T=ellipke(sin(thm/2).^2)*2/pi; %用椭圆积分计算精确周期 n=2^8; %数据个数 wt=linspace(0,2*pi*T,n); %无量纲时间向量 options.RelTol=1e-6; %相对精确
[wt,TH]=ode45('zhou13_1ode_fun',wt,[thm,0],options);%计算角度和角速度 theta=TH(:,1)*180/pi; %角位移向量
fg=figure; %建立图形窗口并取句柄 plot(wt/2/pi,theta,'linewidth',2) %画角位移 grid on %加网格 fs=16; %字体大小 xlabel('\itt\rm/\itT\rm_0','fontsize',fs)%x 标签 ylabel('\it\theta\rm/\circ','fontsize',fs)%y 标签
title('单摆角位移的傅立叶分析','fontsize',fs)%标题
text(0,-thetam,['\itT/T\rm_0=',num2str(T)],'fontsize',fs)%
pause %暂停
hold on%保持图像
h=plot(wt/2/pi,zeros(size(wt)),'r.','erasemode','xor');%画线并取句柄Cn=fft(theta); %离散序列的傅立叶变换
cn=Cn/n; %求复系数
a0=cn(1); %平均值
%text(0,a0,num2str(a0),'fontsize',fs) %显示平均值文本
n=1; %求和下标
x=a0; %常量
An=[]; %频谱变量置空
while 1 %无限循环
an=2*real(cn(n+1)); %取实部
bn=-2*imag(cn(n+1)); %取虚部
An=[An,2*abs(cn(n+1))]; %连接频谱
xn=an*cos(n*wt/T)+bn*sin(n*wt/T); %求分量
x=x+xn; %累加分量
set(h,'ydata',x) %设合成线
plot(wt/2/pi,xn,'color',rand(3,1)) %画分支线
yn=input('还要吗(y/n)?','s'); %键盘输入
if lower(yn)=='n',break,end%判断退出
n=n+1; %波的个数加1
end%结束循环
figure %创建图形窗口
stem(1:n,An) %画杆图
grid on%加网格
xlabel('\it\omega','fontsize',fs) %x标签
ylabel('\itA_n/A','fontsize',fs) %y标签
title('单摆角位移的频谱图','fontsize',fs)%标题
text(1:n,An,num2str(An',4)) %显示频谱
theta=ifft(Cn); %求傅立叶逆变换
figure %创建图形窗口
plot(wt/2/pi,theta) %画角位移曲线
grid on%加网格
xlabel('\itt\rm/\itT\rm_0','fontsize',fs)%x标签
ylabel('\it\theta\rm/\circ','fontsize',fs)%y标签
title('单摆的角位移','fontsize',fs) %标题
程序要调用函数zhou13_1ode_fun.m。

%单摆运动规律的函数
function f=Fun(t,th)
f=[th(2); %角速度
-sin(th(1))]; %角加速度
[说明]程序执行时从键盘输入角振幅。

P5_4_2a图P5_4_2b图。