MATLAB计算方法和技巧6_2阻尼振动.pdf

1 弹簧振子的阻尼振动[问题 ]一弹簧振子的质量为m，倔强系数为k振子还受到与速度大小成正比、方向相反的阻力， 比例系数为γ 当振子从静止开始运动时，初位移为A物体的运动规律是什么？不同的阻尼下的运动曲线和速度曲线有什么特点？[数学模型 ] 根据牛顿运动定律，物体运动的微分方程为22ddddxxmkx tt，(6.2.1) 取 k/m = ω02，γ /m = 2β ， ω 0就是无阻尼时物体的固有角频率，β是阻尼因子物体的运动方程可表示为2 202dd20 ddxxx tt6.2.2) 设微分方程的解为x = ert，代入上式可得特征方程r2- 2βr + ω02= 06.2.3) 特征方程的解为220r，(6.2.4) 设220，α可以是实数和零以及虚数，则r1 = - β + α ，r2 = - β–α ，r1和 r2可以是实数或复数微分方程的解为12 1212eee(ee)r tr ttttxCCCC，(6.2.5) 其中 C1和 C2是由初始条件决定的常数物体的速度为12 1 12212d eee[()e()e] dr tr ttttx vC rC rCC t6.2.6) 当 t = 0 时， x = A，v = 0，因此可得A = C1 + C2，0 = C1(-β + α ) + C2(-β - α )，(6.2.7) 如果 β≠ ω0，即 α≠ 0 ，解得两个常数分别为12CA，22CA。

因此物体的位移为e[()e()e] 2tttAx6.2.8) [讨论 ]①当 β > ω0时，即 α > 0，上式就是过阻尼的情况 ②当 β→ ω0时，即 α → 0，不论用罗必塔法则还是用公式eαt→ 1 + αt和 e-αt→ 1 - αt，都可得0 0(1)etxAt6.2.9) 这是临界阻尼的情况③当 β < ω0时，设220，则 α = iω ，i=-1为虚数单位，利用欧拉公式2 eiθ= cosθ + isinθ ，e-i θ= cosθ - isinθ ，可得i-i1cos(ee) 2，i-i1sin(ee) 2i，由(6.2.8) 式可得iie[(i)e(i)e] 2itttAx，即e( c o ss i n)txAtt，(6.2.10) 或0ec o s ()txAt，(6.2.11) 其中arctan这就是欠阻尼的情况，物体可作准周期性运动，ω是其角频率，振幅按指数规律衰减物体作阻尼运动的周期为2202π2πT，(6.2.12) 可见：阻尼因子越大，周期越长或者说：阻尼使振动变慢了④当 β = 0 时，则 α = iω0，可得00 0[ee]cos 2ititAxAt，可见：在不计阻尼的情况下，物体作简谐振动。

由于 MATLAB 能够进行复数运算，不论是过阻尼、临界阻尼，还是欠阻尼的情况，位移都可以用 (6.2.8)式或 (6.2.10) 式计算 (6.2.8)式可用双曲函数表示e(coshsinh)txAtt6.2.13) 利用欧拉公式可以证明：coshiθ = cosθ ，sinhiθ = isinθ ，因此， (6.2.10) 式和 (6.2.13)式是完全相同的由(6.2.10)式可得物体运动的速度为20esintvAt 6.2.14) [算法 ]方法一：用解析解取A 为坐标的单位取ω0的倒数为时间单位，则约化时间为 t*= ω0t；取约化阻尼因子为β* = β /ω 0，则约化阻尼角频率为220**2001，物体的坐标可表示为* ********exp()(cossin)xttt，(6.2.10*) 其中， x*= x/A可见：阻尼因子β*= β /ω0决定了物体的位移曲线取v0 = ω0A 为速度单位，物体的速度可表示为3 ******01exp() sinvvtt v6.2.14*) 方法二：用微分方程的数值解利用约化物理量，(6.2.2)式可化为2*****2*dd20 ddxxx tt。

6.2.2*) 设 x(1) = x*，x(2) = dx*/dt*，可得*d(1)(2) dxx t，**d(2)2(2)(1) dxxx t6.2.2**) 由于**00d1d(2) ddxxvx tAtA，因此 x(2) 就是约化速度初始条件为x(1) = 1 ，x(2) = 0 方法三：用微分方程的符号解将微分方程(6.2.2*) 式可化为符号形式，利用dsolve 指令可求符号解[程序 ]zqy6_2ode.m 如下 阻尼运动的类型clear % 清除变量b=input(' 请输入相对阻尼因子beta/w0:'); % 键盘输入阻尼因子t=0:0.1:20; % 固有角频率与时间的乘积w0t 向量 ( 约化时间向量 )if b==1,b=b+eps;end% 阻尼因子为 1时则增加一小量w=sqrt(1-b^2); % 准角频率向量x=exp(-b*t).*(cos(w*t)+b/w*sin(w*t)); % 位移函数v=-exp(-b*t).*sin(w*t)/w; % 速度函数figure % 创建图形窗口subplot(2,1,1) % 选子图plot(t,x) % 画位移曲线grid on% 加网格fs=16; % 字体大小xlabel('\it\omega\rm_0\itt', 'FontSize',fs)% 标记横坐标ylabel('\itx/A', 'FontSize',fs) % 标记纵坐标title(' 质点阻尼运动的位移曲线' , 'FontSize',fs)% 标题subplot(2,1,2) % 选子图plot(t,v) % 画速度曲线grid on% 加网格xlabel('\it\omega\rm_0\itt', 'FontSize',fs)% 标记横坐标ylabel('\itv/\omega\rm_0\itA', 'FontSize',fs)% 标记纵坐标title(' 质点阻尼运动的速度曲线' , 'FontSize',fs)% 标题[tm,X]=ode45('zqy6_2fun',t,[1;0],[],b);% 计算位移和速度hold on% 保持图像plot(t,X(:,2),'r.') % 画速度曲线4 subplot(2,1,1) % 选子图hold on% 保持图像plot(t,X(:,1),'r.') % 画位移曲线s=[ 'D2x+' ,num2str(2*b),'*Dx+x']; % 微分方程字符串sx=dsolve(s,'x(0)=1', 'Dx(0)=0'); % 微分方程的符号解sv=diff(sx); % 求速度的符号解x=subs(sx,'t',t); % 位移v=subs(sv,'t',t); % 速度plot(t,x,'ko') % 画位移曲线legend( ' 解析解 ' , ' 数值解 ' , ' 符号解 ' ) % 图例txt=['\it\beta/\omega\rm_0:',num2str(b)];% 参数文本text(0,0,txt,'FontSize',fs) % 显示文本subplot(2,1,2) % 选子图plot(t,v,'ko') % 画速度曲线legend( ' 解析解 ' , ' 数值解 ' , ' 符号解 ' ) % 图例程序要调用函数zqy6_2fun.m 。

阻尼运动的二阶微分方程的函数function f=fun(t,x,flag,b)f=[x(2);-2*b*x(2)-x(1)]; % 速度和加速度[说明 ]用解析解画了曲线之后，可取图片，按任一键继续用微分方程的数值解和符号解所画的点和圈与公式法的曲线完全重合，说明三种方法都是正确的P6.2a 图P6.2b 图P6.2c图作业：用微分方程的数值解计算和绘制单摆阻尼振动的曲线湖南大学物电院周群益。