matlab仿真自由落体

通信仿真技术与实践上机作业一(实例1.1)试对空气中在重力作用下不同质量物体的下落过程进行建模和仿真。已知重力加速度g=9.8m/st,在初始时刻t0=0s时物体由静止开始坠落。考虑空气阻力的影响。(1)建立数学模型质量为m的物体在自由坠落过程中受到竖直向下的恒定重力和向上的空气阻力f的作用，由牛顿第二定律，我们知道，重力mg,加速度a以及物体质量m之间的关系是：mg-f=maf=k*(F2) 卜=空气阻力系数，为一恒定值a=g-k(v^2/m)(2)数学模型的解析分析v(t)=ats(t)=21at^2(3)根据数学模型建立计算机仿真模型(编程)将方程转换为一种在自变量(时间)上的“递推”表达式v(t+dt)=v(t)+dv=v(t)+adts(t+dt)=s(t)+ds=s(t)+v(t)dt(4)执行仿真和结果分析%自由落体.m%模拟受到空气阻力的小球g=9.8;%重力加速度a=g;m=10;k=0.5;%空气阻力系数v=0; %设定初始速度条件s=0; %设定初始位移条件t=0; %设定起始时间dt=0.1;%设置计算步长仿真时间等于N与仿真时间等于N与dt的乘积%计算新时刻的速度%空气阻力f=k*(v"2)/m%新位移%时间更新forf=1:Nv=v+a*dt;a=a-k.*(/2)./m;s(f+1)=s(f)+v*dt;t(f+1)=t(f)+dt;end%作图：受空气阻力落体结果与自由落体结果对比t=0:dt:N*dt;subplot(1,2,2)plot(t,s,，o，);xlabel('时间t');ylabel('位移s,);legendC受空气阻力的落体,)；

运行得到的结果:(5)仿真程序的功能扩展―-以动态方式来观察物体坠落的过程受到空气阻力落体动画.m(实例1.2〕对乒乓球的弹跳过程进行仿真。忽略空气对球的影响，乒乓球垂直下落，落点为光滑的水平面，乒乓球接触落点立即反弹。如果不考虑弹跳中的能量损耗，则反弹前后的瞬时速率不变，但方向相反。如果考虑撞击损耗，则反弹速率有所降低。我们希望通过仿真得出乒乓球位移随时间变化的关系曲线，并进行弹跳过程的“实时”动画显示。(1)数学模型首先对乒乓球弹跳过程进行一些理想化假设。设球是刚性的，质量为m，垂直下落。碰击面为水平光滑平面。在理想情况下碰击无能量损耗。如果考虑碰击面损耗，则碰击前后速度方向相反，大小按比例系0<=K<=1下降。在t时刻的速度设为v=v(t)，位移设为y=y(t)，并以碰击点为坐标原点，水平方向为坐标横轴建立直角坐标系。球体的速度以竖直向上方向为正方向。重力加速度为g=9.8m/s2初始条件假设：设初始时刻t0=0球体的初始速度为v0=v(t0),初始位移为y0=y(t0)。受力分析：在空中时小球受重力F=mg作用，其中，g=-(dt/dv)。则在t+dt时刻小球的速度为v(t+dt)=v(t)-gdt在t+dt时刻小球的位移为y(t+dt)=y(t)+v(t)dt在小球撞击水平面的瞬间，即y(t)=0的时刻，它的速度方向改变，大小按比例K衰减。当K=1时，就是无损耗弹跳情况。因此，小球反弹瞬间(t+dt时刻)的速度为v(t+dt)=-Kv(t)-gdt;0<=K<=1反弹瞬间的位移为y(t+dt)=y(t)-Kv(t)dt=-Kv(t)dt(2)仿真模型设计(程序)从数学模型中可见，小球在空中自由运动时刻与撞击时刻的动力方程不同。通过小球所处位置(位移)是否为零可判定小球处于何种状态。程序文件代码如下。乒乓球弹跳曲线.m%乒乓球弹跳曲线.mg=9.8;%重力加速度v0=0;%初始速度y0=1;%初始位置m=1;%小球质量t0=0;%起始时间K=1; %弹跳的损耗系数N=5000;%仿真的总步进数dt=0.001;%仿真步长v=v0;%初状态y=y0;fork=1:Nif(y>0)|(v>0) %小球在空中的(含刚刚弹起瞬间)动力方程计算v=v-g*dt;y=y+v*dt;else %碰击瞬间的计算y=y-K.*v*dt;v=-K.*v-g*dt;ends(k)=y; %将当前位移记录到s数组中以便作图endt=t0:dt:dt*(N-1);%仿真时间长度subplot(1,2,1);plot(t,s);title('无损耗时的弹跳')；xlabel('时间t');ylabel('位移y(t)');axis([0501.1]);g=9.8;%重力加速度v0=0;%初始速度y0=1;%初始位置m=1;%小球质量t0=0;%起始时间K=0.85;%弹跳的损耗系数N=5000;%仿真的总步进数dt=0.001;%仿真步长v=v0;%初状态y=y0;fork=1:Nif(y>0)|(v>0) %小球在空中的(含刚刚弹起瞬间)动力方程计算v=v-g*dt;y=y+v*dt;else %碰击瞬间的计算y=y-K.*v*dt;v=-K.*v-g*dt;ends(k)=y; %将当前位移记录到s数组中以便作图endt=t0:dt:dt*(N-1);%仿真时间长度subplot(1,2,2);plot(t,s);title('有损耗的弹跳曲线')；xlabel('时间t');ylabel('位移y(t)');axis([0501.1]);

%乒乓球仿真动画.m实例1.3）试用蒙特卡罗方法求出半径为1的圆的面积。并与理论值对比。（1）数学模型设有两个相互独立的随机变量x,y,服从［0,2］上的均匀分布。那么由它们所确定的坐标点（x;y）是均匀分布于边长为2的一个正方形区域中。如下图所示。+>,该正方形的内接圆的半径为1。显然，坐标点区丫)落入圆中的概率p等于该圆面积Sc与正方形面积S之比，即Sc=pS因此，只要通过随机试验统计出落入圆中点的频度，即可计算出圆的近似面积来。当随机试验的次数充分大的时候，计算结果就趋近于理论真值。(2)仿真试验求圆周率1.m%ch1example3prg1.msita=0:0.01:2*pi;x=sin(sita);y=cos(sita);%计算半径为1的圆周上的点，以便作出圆周观察m=0; %在圆内在落点计数器x1=2*rand(1000,1)-1;%产生均匀分布于[-1,+1]直接的两个独立随机数x1,y1y1=2*rand(1000,1)-1;N=1000; %设置试验次数forn=1:N %循环进行重复试验并统计p1=x1(1:n);q1=y1(1:n);if(x1(n)*x1(n)+y1(n)*y1(n))<1%计算落点到坐标原点的距离，判别落点是否在圆内m=m+1; %如果落入圆中，计数器加1endplot(p1,q1,'.',x,y,'-k',[T-111-1],[-111-1T],，-k，);axisequal; %坐标纵横比例相同axis([-22-22]);%固定坐标范围text(-1,-1.2,[，试验总次数n=，,num2str(n)]);%显示试验结果text(-1,-1.4,[，落入圆中数m=，,num2str(m)]);text(-1,-1.6,[，近似圆面积S_c=，,num2str(m/n*4)]);set(gcf,，DoubleBuffer，,，on，); %双缓冲避免作图闪烁drawnow; %显示结果end程序执行中，将动态显示随机落点情况和当前的统计计算结果。随着试验次数增加，计算结果将趋近于半径为1的圆面积的真值pi。蒙特卡罗法求圆周率2.m%蒙特卡罗法求圆周率2.m%ch1example3prg2.mticn=10000;fork=1:1000x1=2*rand(n,1)-1;y1=2*rand(n,1)-1;%计时器启动%每次随机落点10000个%重复试验1000次%随机落点产生m(k)=sum((x1.*x1+y1.*y1)<1);%求落入圆中的点数和endS_c=mean(m).*4./nendS_c=mean(m).*4./n%计算并显示结果time=tocS_c=3.1414time=0.6886这种算法计算速度更快，仿真的结果更准确(实例1.4)实际物理试验中，当我们让一个乒乓球垂直下落到一个完全水平的玻璃板上后，乒乓球不断弹跳，直到能量耗尽。我们要建立更加接近真实的物理环境的弹跳模型，就必须考虑被忽视的微小的扰动因素，根据大数定理，在数学上我们就可以将水平作用力建模为一个高斯随机变量。(1)分析在三维空间内仿真小球下落的过程：ax(t)=Fx(t)/mdvx(t)=ax(t)dtdsx(t)=vx(t)dtY、2平面也如此

（2）仿真代码三维弹跳动画.m%三维弹跳动画g=9.8; %重力加速度v0=0; %初始速度y0=1; %初始位置m=0.1; %小球质量t0=0; %起始时间K=0.8; %弹跳的损耗系数dt=0.005;%仿真步长初状态N=5000; %仿真的总步进数dt=0.005;%仿真步长初状态v=v0y=y0vx=0vz=0sx=0sz=0fork=1:Nify>0fork=1:Nify>0%小球在空中的动力方程计算v=v-g*dt;y=y+v*dt;else y=y+v*dt;else %碰击瞬间的计算y=-K.*v*dt;v=-K.*v-g*dt;endFx=randn; % x水平方向的随机力，方差为1ax=Fx./m; % Fx导致的x水平方向的加速度vx=vx+ax*dt; %小球在x水平方向的瞬时速度sx=sx+vx*dt; %小球在*水平方向的位移Fz=randn; % z水平方向的随机力，方差为1az=Fz./m; % Fz导致的z水平方向的加速度vz=vz+az