matlab仿真自由落體

通信仿真技術(shù)與實踐上機作業(yè)一(實例1.1)試對空氣中在重力作用下不同質(zhì)量物體的下落過程進行建模和仿真。已知重力加速度g=9.8m/st,在初始時刻t0=0s時物體由靜止開始墜落。考慮空氣阻力的影響。(1)建立數(shù)學(xué)模型質(zhì)量為m的物體在自由墜落過程中受到豎直向下的恒定重力和向上的空氣阻力f的作用，由牛頓第二定律，我們知道，重力mg,加速度a以及物體質(zhì)量m之間的關(guān)系是：mg-f=maf=k*(F2) 卜=空氣阻力系數(shù)，為一恒定值a=g-k(v^2/m)(2)數(shù)學(xué)模型的解析分析v(t)=ats(t)=21at^2(3)根據(jù)數(shù)學(xué)模型建立計算機仿真模型(編程)將方程轉(zhuǎn)換為一種在自變量(時間)上的“遞推”表達式v(t+dt)=v(t)+dv=v(t)+adts(t+dt)=s(t)+ds=s(t)+v(t)dt(4)執(zhí)行仿真和結(jié)果分析%自由落體.m%模擬受到空氣阻力的小球g=9.8;%重力加速度a=g;m=10;k=0.5;%空氣阻力系數(shù)v=0; %設(shè)定初始速度條件s=0; %設(shè)定初始位移條件t=0; %設(shè)定起始時間dt=0.1;%設(shè)置計算步長仿真時間等于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%作圖：受空氣阻力落體結(jié)果與自由落體結(jié)果對比t=0:dt:N*dt;subplot(1,2,2)plot(t,s,，o，);xlabel('時間t');ylabel('位移s,);legendC受空氣阻力的落體,)；

運行得到的結(jié)果:(5)仿真程序的功能擴展―-以動態(tài)方式來觀察物體墜落的過程受到空氣阻力落體動畫.m(實例1.2〕對乒乓球的彈跳過程進行仿真。忽略空氣對球的影響，乒乓球垂直下落，落點為光滑的水平面，乒乓球接觸落點立即反彈。如果不考慮彈跳中的能量損耗，則反彈前后的瞬時速率不變，但方向相反。如果考慮撞擊損耗，則反彈速率有所降低。我們希望通過仿真得出乒乓球位移隨時間變化的關(guān)系曲線，并進行彈跳過程的“實時”動畫顯示。(1)數(shù)學(xué)模型首先對乒乓球彈跳過程進行一些理想化假設(shè)。設(shè)球是剛性的，質(zhì)量為m，垂直下落。碰擊面為水平光滑平面。在理想情況下碰擊無能量損耗。如果考慮碰擊面損耗，則碰擊前后速度方向相反，大小按比例系0<=K<=1下降。在t時刻的速度設(shè)為v=v(t)，位移設(shè)為y=y(t)，并以碰擊點為坐標(biāo)原點，水平方向為坐標(biāo)橫軸建立直角坐標(biāo)系。球體的速度以豎直向上方向為正方向。重力加速度為g=9.8m/s2初始條件假設(shè)：設(shè)初始時刻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衰減。當(dāng)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)仿真模型設(shè)計(程序)從數(shù)學(xué)模型中可見，小球在空中自由運動時刻與撞擊時刻的動力方程不同。通過小球所處位置(位移)是否為零可判定小球處于何種狀態(tài)。程序文件代碼如下。乒乓球彈跳曲線.m%乒乓球彈跳曲線.mg=9.8;%重力加速度v0=0;%初始速度y0=1;%初始位置m=1;%小球質(zhì)量t0=0;%起始時間K=1; %彈跳的損耗系數(shù)N=5000;%仿真的總步進數(shù)dt=0.001;%仿真步長v=v0;%初狀態(tài)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; %將當(dāng)前位移記錄到s數(shù)組中以便作圖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;%小球質(zhì)量t0=0;%起始時間K=0.85;%彈跳的損耗系數(shù)N=5000;%仿真的總步進數(shù)dt=0.001;%仿真步長v=v0;%初狀態(tài)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; %將當(dāng)前位移記錄到s數(shù)組中以便作圖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）數(shù)學(xué)模型設(shè)有兩個相互獨立的隨機變量x,y,服從［0,2］上的均勻分布。那么由它們所確定的坐標(biāo)點（x;y）是均勻分布于邊長為2的一個正方形區(qū)域中。如下圖所示。+>,該正方形的內(nèi)接圓的半徑為1。顯然，坐標(biāo)點區(qū)丫)落入圓中的概率p等于該圓面積Sc與正方形面積S之比，即Sc=pS因此，只要通過隨機試驗統(tǒng)計出落入圓中點的頻度，即可計算出圓的近似面積來。當(dāng)隨機試驗的次數(shù)充分大的時候，計算結(jié)果就趨近于理論真值。(2)仿真試驗求圓周率1.m%ch1example3prg1.msita=0:0.01:2*pi;x=sin(sita);y=cos(sita);%計算半徑為1的圓周上的點，以便作出圓周觀察m=0; %在圓內(nèi)在落點計數(shù)器x1=2*rand(1000,1)-1;%產(chǎn)生均勻分布于[-1,+1]直接的兩個獨立隨機數(shù)x1,y1y1=2*rand(1000,1)-1;N=1000; %設(shè)置試驗次數(shù)forn=1:N %循環(huán)進行重復(fù)試驗并統(tǒng)計p1=x1(1:n);q1=y1(1:n);if(x1(n)*x1(n)+y1(n)*y1(n))<1%計算落點到坐標(biāo)原點的距離，判別落點是否在圓內(nèi)m=m+1; %如果落入圓中，計數(shù)器加1endplot(p1,q1,'.',x,y,'-k',[T-111-1],[-111-1T],，-k，);axisequal; %坐標(biāo)縱橫比例相同axis([-22-22]);%固定坐標(biāo)范圍text(-1,-1.2,[，試驗總次數(shù)n=，,num2str(n)]);%顯示試驗結(jié)果text(-1,-1.4,[，落入圓中數(shù)m=，,num2str(m)]);text(-1,-1.6,[，近似圓面積S_c=，,num2str(m/n*4)]);set(gcf,，DoubleBuffer，,，on，); %雙緩沖避免作圖閃爍drawnow; %顯示結(jié)果end程序執(zhí)行中，將動態(tài)顯示隨機落點情況和當(dāng)前的統(tǒng)計計算結(jié)果。隨著試驗次數(shù)增加，計算結(jié)果將趨近于半徑為1的圓面積的真值pi。蒙特卡羅法求圓周率2.m%蒙特卡羅法求圓周率2.m%ch1example3prg2.mticn=10000;fork=1:1000x1=2*rand(n,1)-1;y1=2*rand(n,1)-1;%計時器啟動%每次隨機落點10000個%重復(fù)試驗1000次%隨機落點產(chǎn)生m(k)=sum((x1.*x1+y1.*y1)<1);%求落入圓中的點數(shù)和endS_c=mean(m).*4./nendS_c=mean(m).*4./n%計算并顯示結(jié)果time=tocS_c=3.1414time=0.6886這種算法計算速度更快，仿真的結(jié)果更準(zhǔn)確(實例1.4)實際物理試驗中，當(dāng)我們讓一個乒乓球垂直下落到一個完全水平的玻璃板上后，乒乓球不斷彈跳，直到能量耗盡。我們要建立更加接近真實的物理環(huán)境的彈跳模型，就必須考慮被忽視的微小的擾動因素，根據(jù)大數(shù)定理，在數(shù)學(xué)上我們就可以將水平作用力建模為一個高斯隨機變量。(1)分析在三維空間內(nèi)仿真小球下落的過程：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; %小球質(zhì)量t0=0; %起始時間K=0.8; %彈跳的損耗系數(shù)dt=0.005;%仿真步長初狀態(tài)N=5000; %仿真的總步進數(shù)dt=0.005;%仿真步長初狀態(tài)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導(dǎo)致的x水平方向的加速度vx=vx+ax*dt; %小球在x水平方向的瞬時速度sx=sx+vx*dt; %小球在*水平方向的位移Fz=randn; % z水平方向的隨機力，方差為1az=Fz./m; % Fz導(dǎo)致的z水平方向的加速度vz=vz+az