MATLAB计算方法和技巧9_1单摆的傅立叶分析.doc

单摆振动规律的傅立叶分析问题单摆振动规律的傅立叶分析对单摆振动规律进行傅立叶分析。数学模型如A9.1图所示，设摆锤质量为m，角位置为，摆锤的运动方程为lmgOftTA9.1图，即 ， (9.1.1)在小角度的情况下，sin ，可得 ， (9.1.2)其中，0为小角振幅的圆频率。可知：单摆在小角度时作简谐振动，小角度周期为。 (9.1.3)可见：在小角振动的情况下，单摆的周期与角振幅无关，这称为单摆的等时性。摆锤的角速度为 = d/dt，因此，由(9.1.1)式可得，积分得，当t = 0时， = 0， = m，可得C = -gcosm/l。因此角速度大小为。 (9.1.4)单摆的周期为。 (9.1.5)对于角振幅m，通过数值积分和符号积都能计算周期。利用半角公式可得。设， (9.1.6)并设ksinx = sin(/2)，因此，可得，即 。 (9.1.7)这是椭圆积分。第一类完全椭圆积分定义为， (9.1.8)周期为。 (9.1.9)算法对于角振幅m，利用(9.1.5)式，通过MATLAB数值积分指令quadl和符号积分指令int都可计算单摆的周期。利用MATLAB的完全椭圆积分指令ellipke也能计算单摆的周期。不过，MATLAB定义的一类完全椭圆积分定义为， (9.1.8*)周期可表示为。 (9.1.9*)其中。 (9.1.6*)取t* = 0t为无量纲时间，(9.1.1)式可化为， (9.1.1*)取(1) = ，(2) = d/dt*，可得，。初始条件是(1) = m，(2) = 0。由于，因此(2)表示以0为单位的角速度。由于，因此，无量纲时间除以2就是以小角周期T0为单位的时间。程序zqy9_1fourier.m如下。%单摆振动规律的傅立叶分析clear %清除变量thetam=input(请输入单摆角振幅的度数:);%键盘输入角振幅thm=thetam*pi/180; %化为弧度T=ellipke(sin(thm/2).2)*2/pi; %用椭圆积分计算精确周期n=28; %数据个数wt=linspace(0,2*pi*T,n); %无量纲时间向量options.RelTol=1e-6; %相对精确wt,TH=ode45(zqy9_1fun,wt,thm,0,options);%计算角度和角速度theta=TH(:,1)*180/pi; %角位移向量fg=figure; %建立图形窗口并取句柄plot(wt/2/pi,theta,.-,LineWidth,2) %画角位移grid on %加网格fs=16; %字体大小xlabel(ittrm/itTrm_0,FontSize,fs)%x标签ylabel(itthetarm/circ,FontSize,fs)%y标签 title(单摆角位移的傅立叶分析,FontSize,fs)%标题text(0,thetam,num2str(thetam),FontSize,fs) %显示平均值文本text(0,-thetam,itT/Trm_0=,num2str(T),FontSize,fs)%pause %暂停hold on %保持图像h=plot(wt/2/pi,zeros(size(wt),r,LineWidth,2);%画线并取句柄Cn=fft(theta); %离散序列的傅立叶变换cn=Cn/n; %求复系数a0=cn(1); %平均值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),LineWidth,2)%画分量线 yn=input(还要吗(y/n)?,s); %键盘输入 if lower(yn)=n,break,end %判断退出 n=n+1; %波的个数加1end %结束循环figure %创建图形窗口stem(1:n,An) %画杆图grid on %加网格xlabel(itomega/omegarm_0,FontSize,fs)%x标签ylabel(itA_n/A,FontSize,fs) %y标签 title(单摆角位移的频谱图(itthetarm_m=,. num2str(thetam),circ),FontSize,fs)%标题text(1:n,An,num2str(An,4),FontSize,fs)%显示频谱 theta=ifft(Cn); %求傅立叶逆变换figure(fg) %重开图形窗口plot(wt/2/pi,theta,o) %画角位移曲线程序