数学实验：7-数值仿真

1、实验六 微分方程数值仿真一、试验目的与要求1动态仿真单自由度阻尼系统，将阻尼系数的影响可视化；2直接仿真弹簧振动，观察阻尼系数对振动情形的影响；3熟练掌握ode指令的使用方法；4学会将方程框图化，初步了解 Simulink仿真的思想设计；5学习实时动画编程的原理。二、问题描述1振动曲线的可视化：考虑无强迫力的情形 ，固定，初始位置，初始速度，计算终点时间，分析单自由度阻尼系统的阻尼系统对其固有模态的影响。2模拟弹簧振动，讨论最简单的情形，一弹簧系统做横向运动，其位移由振动规律所决定，仿真弹簧振动。3方程的框图表示和 Simulink仿真三、问题分析对问题1，可以用ode45求解微分方程的解，然

2、后绘制弹簧振动曲线，这要求把该方程先转化位一阶微分方程组：令，则方程化为其中，u(t)为强迫项，在观察程序中选取初始项对问题2，运行观察程序exp2_3.m，图形窗口将显示一弹簧块的动画振动（图 是其中的一帧）。对问题3，考虑微分方程 ，其中u(t)是一个输入信号函数。四、知识背景1单自由度阻尼系统研究汽车系统减振问题时，我们对汽车在道路上行驶时的垂直位移x(t)感兴趣。对于所关心的问题，可以假定整个车体的质量集中在理想的质量m上，并用一个理想的弹簧k和减震器a分别表示汽车的弹性和冲击阻尼，这样一个质量弹簧减震器(m-k-a)如图1所示。图2-1 质量弹簧减震器假定该机械系统的3种元件都是线性

3、的，则可用下面的二阶线性微分方程来描述该系统： (1)其中称为二阶线性系统的阻尼系数，称为二阶系统的稳态增益，为二阶系统的自然频率。实验指出，当物体以不太大的速度在介质中运动时，介质对物体的阻力与物体运动速度成正比，方向总和速度方向相反。在图1的系统中，系统的输入量为外力f(t),输出量为质量m的位移x(t)。根据牛顿第二定律，在外力作用下质量为m的物体将产生加速度，于是由虎克定律，将弹簧的恢复力，阻尼力代入上式得，这样就得到二阶常系数微分方程在方程两边同除以m得，令，则可将上式化为（1）的形式。 可以用解析方法得到方程（1）的解，从而由显示解的表达式直接分析这种影响，这是最好的方法。 但是我

4、们感兴趣的是仿真方法，即用曲线将阻尼系数的影响可视化，或者直接仿真弹簧的振动，观察阻尼系数对振动情形的影响。2实时动画编程实时动画是通过设置绘图指令Plot的擦除属性Erasemode来实现的。（1）属性Erasemode的选项：实时动画是计算新对象并予以显示，同时擦除旧对象。Matlab通过改变画图命令的擦除方式的属性来达到这一目的。例如命令：Plot(x,y,Erasemode,background)Erasemode=background 其含义是画图时将旧对象变为背景色，同时画出新对象，达到动画显示的目的。Normal 重画整个显示区None 不作任何擦除Xor 只画和屏幕颜色不一致的

5、新对象点，只擦除和屏幕颜色不一致的旧对象点（2）句柄图形、set 和get 指令 ：在设计好弹簧的位置给出相应的数据后，可通过使用句柄图形命令获得图形：；但是现在并没有画图令。如果要画图应使用drawnow命令，同时用将当前轴对象gca的UserData属性设置为画图句柄hndl，这样设置好以后，当需要画图时，只需要用get指令就获得这一句柄，用drawnow就可以画图了。可用set指令改变句柄图形中的数据，例如指令是将新的横坐标数据赋予当前图形句柄hndl的绘图向量x中。重复使用set和drawnow就可以得到动画效果。3Simulink仿真工具简介框图是方程的等价物，既使是普通的微分方程，

6、也可以用框图来描述。在工程设计和计算中，人们往往用框图来表示复杂的系统。因为框图可形象的把各个环节区分开来，便于研究各个环节的参数对系统的影响，simulink仿真工具箱提供了这样的手段。Simulink是基于框图化思想实现建模和仿真的，主要用于微分和差分方程描写的非线性动态系统。五、实验过程1阻尼系数对振动曲线的影响观察1）选取阻尼系数c的一系列值，对这些阻尼系数用数值方法求解微分方程得到相应的数值解，将这些解曲线绘制在同一图形上，标注相应阻尼系数的值，观察这些曲线的变化和差异，比较不同阻尼系数对振幅曲线的影响。运行exp2_1.m 图显示结果。图2-2 不同阻尼系数对应的振动曲线2）另一种

7、方法是时间阻尼系数振动曲线构成的三维空间中表示振动模态和阻尼系数的关系。使用彗星指令comet3。运行exp2_2.m, 结果显示在图中， 图中的彗星正在运动。图2-3 在三维空间中逐一观察振动曲线观察程序以及说明如下：exp2_1.m （图2-2） clear;clf; global c wx0=1;0;w=10;tspan=linspace(0,4,100);cc=0.1,0.4,0.7,1;hold on, xlabel('t'), ylabel('x');for i=1:length(cc); c=cc(i); t,x=ode45('exp2_1

8、f',tspan,x0); %求微分方程的数值解。 text(t(10),x(10,1),'leftarrow c=',num2str(c),'fontsize',15) plot(t,x(:,1),endexp2_1f.m od45指令中的微分方程组函数子程序 function dx=exp2_1f(t,x) %建立微分方程global c wu=0;dx=0;0;dx(1)=x(2);dx(2)=u-c*w*x(2)-w2*x(1); exp2_2.m （图2-3）x0=1;0; w=10; tspan=linspace(0,4,1500); cc=

9、1:-1/10:0;xx=;for i=1:length(cc); c=cc(i);t,x=ode45('exp2_1f',tspan,x0); xx=x(:,1),xx;endaniminit('one'); %逐条观察振动图形。 for i=1:length(cc) c=cc(i); plot3(t,c*ones(length(t),1),xx(:,i),'r:'),hold on, view(30,60); %适当选择观察角度 comet3(t,c*ones(length(t),1),xx(:,i);axis( 0, 4, -0.2, 1.

10、2, -1.1, 1.1)end2振动弹簧的实时动画1）确定弹簧的位置为简单起见，用一条直线代替弹簧。在弹簧的平衡点u=2处，最左端垂线、直线（弹簧）和矩形块的x-y坐标依次是xy=0, 0, 0, 0, u, u, u+1, u+1, u, u ; -1.2, 0, 1.2, 0, 0, 1.2, 1.2, -1.2, -1.2, 0 ； 其中x=xy(1,:)、y=xy(2,:)分别为这些点的横坐标和纵坐标。用画上一条直线代表地板。当弹簧振动时纵坐标不变，弹簧的横坐标随时间按u=2+cos(t)的规律运动，所以在每一时刻只需重新计算xy(1,:)中的u的新坐标即可。2）实时动画编程观察程序

11、exp2_3.m以及说明： animinit('onecart1 Animation') axis(-2 6 -10 10); hold on; u=2; xy= 0, 0, 0, 0, u, u, u+1, u+1, u, u; -1.2, 0, 1.2, 0, 0, 1.2, 1.2, -1.2, -1.2, 0; x=xy(1,:);y=xy(2,:); % Draw the floor under the sliding masses plot(-10 20,-1.4 -1.4,'b-','LineWidth',2); hndl=plot

12、(x,y,'b-','EraseMode','XOR','LineWidth',2); set(gca,'UserData',hndl); for t=1:0.025:100 u=2+cos(t); x=0, 0, 0, 0, u, u, u+1, u+1, u, u; hndl=get(gca,'UserData'); set(hndl,'XData',x); pause(0.005); drawnowend 运行观察程序exp2-3.m，图形窗口将显示一弹簧块的动画振动（图2-4

13、是其中一祯）。图24振动的弹簧3方程的框图表示及Simulink仿真图2-5是用Smulink做成的模拟框图 图2-5 框图实现Simulink创建上述模型的步骤是（1） 在命令窗口键入Simulink时，弹出Simulink Library Brower窗口。或者用鼠标在命令窗口的工具栏单击Simulink Library Brower图标（用鼠标在该图标上停留片刻，下方将显示Simulink Library Brower的说明，下面找图标的方法也是一样的）（2） 单击Simulink Library Brower工具条上的creat a new model 图标，弹出模型编辑窗口；（3）

14、在Simulink Library Brower窗口中的Simulink项单击右键，然后单击左键open the simulink library条，这时弹出Simulink Block Library窗口。（4） 双击左键continuous图标，弹出continuous block窗口。用鼠标按住积分器integrator图标将其拖动模型文件窗口后松开左键，这时integrator框便放入模型文件中，等待编辑。（5） 同样在Math下找到gain和sum图标，在sources下找到signal generator图标，在sink下找到scope图标，将它们分别移到模型窗口中。（6） 在模型

15、窗口中将这些基本图标按框图移到所需要的位置，然后直接用鼠标在窗口中画出它们之间的连线，这很容易做到（在连线需要转折的地方，只要单击一下鼠标左键然后继续连线即可）。在需要标注文字的位置单击左键，然后直接输入文字，十分方便。（7） 在Gain图标上单击右键弹出一菜单，在Parameter下可选择所需参数，在Format下可对该图标做旋转等操作，其它情形类似。画好框图并设定好参数后，就可以进行模拟了。点击filesave保存文件，然后点击Simulation Start 开始模拟。双击scope图标将输出模拟的图形，呈现示波器样，如图所示。图2-6 示波器的输出回到弹簧的振动问题，方程为可用框图表示

16、这个方程，开始仿真后，屏幕将出现振动的弹簧，双击Actual position 可观察波动图形。下图是弹簧振动的模拟框图。图2-7 弹簧振动的模拟框图六、结论与应用1图2-2 中每一条曲线对应着弹簧在某个阻尼系数下的振动波形，可以看到，随着阻尼系数的增加，曲线振动的衰减也越来越快，但是当曲线过多时这种方法不容易观察清楚。使用彗星指令comet3，可以逐条的观察不同阻尼系数所对应的振动曲线运行的情况。2通过动画仿真可以更直观地感受到阻尼系数对弹簧振动的影响。3Signal Generator（信号发生器）模块产生了输入信号u，sum模块将u和-2x加在一起得到-2x+u，由方程有，通过积分器integrator得到x，图中integrator 模块用1/s表示，这可由式解释。X一方面倒输出终端scope模块，显示输出结果，另一方面，通过Gain（增益）模块乘以增