常微分MATLAB质点振动

质点振动振动是日常生活和工程技术中常见的一种运动形式。利用常系数线性微分方程的理论来讨论有关自由振动和强迫振动的相关问题。利用MATLAB数学软件大致可分四类情况：（1）无阻尼自由振动情况；（2）有阻尼自由振动；（3）无阻尼强迫振动；（4）有阻尼强迫振动求其数值解和解析解；1、MATLAB软件求解微分方程解析解的命令“dsolve()”(1)求通解的命令格式：（’微分方程’,’自变量’）注：微分方程在输入时，一阶导数y’应输入Dy，y’’应输入D2y等，D应大写。（I）无阻尼自由振动情况常见的数学摆的无阻尼微小振动方程代码如下：t=0:pi/50:2*pi;y=2*sin(3*t+2);plot(t,y,'b')01234567-2-1.5-1-0.500.511.52（III）无阻尼强迫振动用ＭＡＴＬＡＢ操作的代码如下：t=0:0.1:10;y=3*exp(-0.1*t).*sin(5*t+4);plot(t,y,'k-')2，有阻尼自由振动由无阻尼振动的通解可以看出，无阻尼振动是按照正弦规律运动的，摆动似乎可以无限期的进行下去，但事实上，空气从在阻力，在运动时，我们必须把空气阻力考虑在内，所以我们得到有阻尼摆动方程为：记u/m=2n,g/l=w^2,这里n,w是正常数，所以：y=dsolve('D2y+2*n*Dy+w^2*y=0','t');（4.43）解得:y=C3*exp(-t*(n+((n+w)*(n-w))^(1/2)))+C2*exp(-t*(n-((n+w)*(n-w))^(1/2)))（1）小阻尼情形:nw时，方程（4.43）的通解为：y=exp(-n*t)*(c1*cos(w1*t)+c2*sin(w1*t))和前面无阻尼的情形一样，可以把上式的通解改写为一下形式：y=A*exp(-n*t)*sin(w1*t+Q),(4.45)这里的A，Q为任意常数。用matlab操作得到：t=0:0.1:10;y=3*exp(-0.1*t).*sin(5*t+4);plot(t,y,'k-')如图：012345678910-3-2-10123由（4.45）可见，摆动的运动不是周期的，振动的幅度随着时间的增加而不断减小。（2）大阻尼情形：你w时；r2r10;方程（4.43）的通解：Y=c1*exp(r1*t)+c2*exp(r2*t),这里的c1,c2为任意常数；（3）临界的情形：即，n=w的情形，方程（4.43）的通解解为y=exp(-n*t)*(c1+c2*t),这里的，c1,c2为任意常数；由MAYLAB绘制图像得：t=0:0.1:100;yy=exp(-0.2*t).*(0.5-0.2*t);plot(t,yy,'r-')0102030405060708090100-0.3-0.2-0.100.10.20.30.40.53有阻尼自由振动无阻尼自由振动和有阻尼自由振动都属于自由振动，它对应于一个二阶常数齐次线性微分方程。当一个振动系统还经常受到一个外力作用时，这种振动称为强迫振动。＝Ａｓｉｎ（ｔ＋）＋Ｈ／（２－ｐ２）ｓｉｎ（ｐｔ）取Ａ＝２；＝５；ｐ＝３；t=0:pi/50:2*pi;y=2*sin(5*t+2)+1/9*sin(3*t);plot(t,y,'k')01234567-2.5-2-1.5-1-0.500.511.522.54.有阻尼强迫振动这时摆动的运动方程（1.11）变为：D2y+2*n*Dy+w.^2*y=H*sin(p*t),（4.52）根据实际情况，我们只讨论小阻尼的情况，即nw的情形，这时的（4.52）对应的齐次线性方程的通解为：Y=A*exp(-n*t)*sin(w1*t+Q),(4.45)这里的A，Q为任意常数,w1=sqrt(w^2-n^2)现在求（4.52）的一个特解，这时可以寻求如：Y=M*cos(p*t)+N*sin(p*t)的特解，这里M，N是待定系数，代入（4,53），（4,52）化简得：Y=H*sin(Q)*cos(p*t)+H*cos(Q)*sin(P*t)=H*sin(p*t+Q),因此，（4.52）的通解为：Y=A*exp(-n*t)*sin(w1*t+Q)+H/sqrt((w^2-p^2)^2+4*n^2*p^2)*sin(p*t+Q);Matlab操作如下：A=3;n=0.2;w1=2;Q=0.5;H=5;w=3;p=2;t=0:0.1:50;Y=A*exp(-n*t).*sin(w1*t+Q)+H/sqrt((w.^2-p^2)^2+4.*n.^2.*p^2).*sin(p.*t+Q);plot(t,Y,'k-')05101520253035404550-3-2-101234