matlab仿真报告.docx
《matlab仿真报告.docx》由会员分享，可在线阅读，更多相关《matlab仿真报告.docx（16页珍藏版）》请在冰点文库上搜索。
matlab仿真报告
下面我們通過對自由落體的仿真試驗來說明計算機仿真的過程。
〔實例1.1〕試對空氣中在重力作用下不同質量物體的下落過程進行建模和仿真。
已知重
力加速度g=9:
8m/s^2，在初始時刻t0=0s時物體由靜止開始墜落。
空氣對落體的影響可以
忽略不計。
（1）建立數學模型
首先，我們根據物理原理建立自由落體的數學模型。
在空氣阻力可忽略不計時，質量為m
的物體在自由墜落過程中受到豎直向下的恆定重力的作用，由牛頓第二定律，我們知道，重
力F，加速度a以及物體質量m之間的關係是：
F=ma（1.1）
其中加速度就是重力加速度，即a=g。
根據題設，初始時刻為t0=0，物體的初始速度
為v（t0）=0，並設物體下落的瞬時速度為v（t）。
設物體在t時刻的位移為s（t），並設初始位移為零，即s（t0）=0。
根據加速度、速度、位移三者之間的微積分關係，我們得到一組數學方程
a=dv/dt
v=ds/dt
F=ma
以及初始條件（也稱為方程的邊界條件）
v（t0）=0（1.5）
s（t0）=0m（1.6
我們需要得出不同時刻物體的運動狀態，即物體的瞬時速度和瞬時位移。
至此，我們得到
了自由落體的數學描述。
（2）數學模型的解析分析
數學模型建立之後，可以嘗試對其進行解析求解。
解析結果可以幫助我們驗證仿真數值
結果。
對於這個數學模型，其求解十分簡單。
只要對加速度方程、速度方程進行積分並代入初始條件，就能得到我們熟知的結果：
（3）根據數學模型建立計算機仿真模型（編程）
計算機仿真就是對數學模型的數值求解。
我們將微分方程進行形式上的變換以便於數值
求解。
由（1.2）和（1.3）式得
v（t+dt）=v（t）+dv=v（t）+adt（1.9）
和
s（t+dt）=s（t）+ds=s（t）+v（t）dt（1.10）
注意，這種變形僅僅是將方程轉換為一種在自變量（時間）上的「遞推」表達式，並沒有進行
解析求解。
利用（1.9）和（1.10），在已知當前時刻t的瞬時位移，瞬時速度和加速度的情況下，我們就可以推知下一個無限鄰近的時刻t+dt上物體新的瞬時位移，瞬時速度和加速度，這也就是微分方程數值求解的基本思想。
在數值求解中，無窮小量dt需要用一個很小的數量¢t來近似，¢t稱為微分方程的數值求解步長，通常也稱為仿真步進。
顯然，這種微分方程的遞推求解總是近似的，求解精度與步長有關。
下面我們就用程序來實現這個求解過程。
仿真時間範圍設置為0到2秒，為了使得仿真計算的誤差明顯一些，我們故意採用了較大的仿真步長¢t=0:
1。
讀者可以自己修改仿真步長來觀察計算精度的變化情況。
我們將瞬時位移作為仿真輸出變量，並同時通過（1.8）式計算出解析結果以相互對照。
（4）執行仿真和結果分析
g=9.8;%重力加速度
v=0;%設定初始速度條件
s=0;%設定初始位移條件
t=0;%設定起始時間
dt=0.1;%設置計算步長
N=20;%設置仿真遞推次數.仿真時間等於N與dt的乘積
fork=1:
N
v=v+g*dt;%計算新時刻的速度
s（k+1）=s（k）+v*dt;%新位移
t（k+1）=t（k）+dt;%時間更新
end
%理論計算,以便與仿真結果對照
t_theory=0:
0.01:
N*dt;%設置解析計算的時間點
v_theory=g*t_theory;%解析計算的瞬時速度
s_theory=1/2*g*t_theory.^2;%解析計算的瞬時位移
%作圖:
仿真結果與解析結果對比
t=0:
dt:
N*dt;
plot（t,s,'o',t_theory,s_theory,'-'）;
xlabel（'時間t'）;ylabel（'位移s'）;
legend（'仿真結果','理論結果'）;
仿真程序編寫完成並調試正確之後，運行得到的結果如圖1.1所示。
從圖中可見，仿真得
出的位移與理論結果之間存在差別，這種差別是由於微分方程數值求解的算法和採用步長較
大而引起的。
事實上，我們這裡採用的數值求解算法是最簡單的矩形積分的方法，精度不高，
僅僅是為了說明仿真過程而已。
現代仿真技術和數值計算方法中已經開發出許多更好的微分
方程求解算法，可供直接利用。
（5）仿真程序的功能擴展
從ch1example1prg1.m的程序代碼中可見，我們採用了循環語句來實現對微分方程的遞
推求解，每次循環就將計算時刻向前推進一個步長。
全部循環執行完畢後，就得到了一系列時刻上物體的瞬時速度和瞬時位移值，最後通過繪圖語句將結果數據用曲線表達出來。
如果我們希望以動態方式來觀察物體墜落的過程，可以這樣設計仿真程序：
使得在數值求
解的過程中能將求解結果以圖形方式輸出出來。
這樣，在數值求解不斷更新的過程中輸出圖
形也隨之同步更新，形成一種「動畫」的效果。
這種一邊計算一邊輸出可視化結果的方式更
加形象直觀，更便於理解物理系統的工作過程，同時也方便演示，教學講解和學術交流。
我們可以將作圖語句放在遞推計算循環內，並設置即時作圖刷新方式，從而得到這種「動
畫」仿真的效果。
當然，這樣是以犧牲計算速度為代價的。
修改後的程序文件代碼如下。
〔程序代碼〕ch1example1prg2.m
%ch1example1prg2.m
g=9.8;%重力加速度
forL=1:
5%仿真重複5次以便於觀察
v=0;%初始速度
s=0;%初始位置
t=0;
dt=0.01;%計算步長
fork=1:
200
v=v+g*dt;%速度
s=s+v*dt;%位移
t=t+dt;%時間
plot（0,-s,'o'）;
axis（[-22-200]）;%坐標範圍固定
text（0.5,-1,['当前时间：
t=',num2str（t）]）;
text（0.5,-2,['当前速度：
v=',num2str（v）]）;
text（0.5,-3,['当前位置：
s=',num2str（s）]）;
set（gcf,'DoubleBuffer','on'）;%雙緩衝避免作圖閃爍
drawnow;%立即作圖
end
end
在程序中，計算步進重新設置為0:
01，並且重複仿真多次以便於演示。
在作圖語句之後設
置了固定的顯示坐標範圍，並通過text語句在圖上動態顯示當前計算時刻的結果值。
語句
set（gcf,'DoubleBuffer','on'）;
將顯示設備雙緩衝開啟，以避免屏幕刷新引起的閃爍。
而
drawnow
是立即執行繪圖指令。
在仿真中途可按鍵Ctrl+C來終止程序執行。
程序執行過程中將顯示出物體墜落的動畫效果，其中的一個幀如圖1.2所示。
實例1.2〕對乒乓球的彈跳過程進行仿真。
忽略空氣對球的影響，乒乓球垂直下落，落點
為光滑的水平面，乒乓球接觸落點立即反彈。
如果不考慮彈跳中的能量損耗，則反彈前後的瞬時速率不變，但方向相反。
如果考慮撞擊損耗，則反彈速率有所降低。
我們希望通過仿真得出乒乓球位移隨時間變化的關係曲線，並進行彈跳過程的「實時」動畫顯示。
（1）數學模型
首先對乒乓球彈跳過程進行一些理想化假設。
設球是剛性的，質量為m，垂直下落。
碰
擊面為水平光滑平面。
在理想情況下碰擊無能量損耗。
如果考慮碰擊面損耗，則碰擊前後速
度方向相反，大小按比例係數K;0在t時刻的速度設為v=v（t），位移設
為y=y（t），並以碰擊點為坐標原點，水平方向為坐標橫軸建立直角坐標系。
球體的速度以豎
直向上方向為正方向。
重力加速度為g=9:
8m=s2。
初始條件假設：
設初始時刻t0=0球體的初始速度為v0=v（t0），初始位移為y0=y（t0）。
受力分析：
在空中時小球受重力F=mg作用，其中，g=-dv/dt。
則在t+dt時刻小球的速度為（注意其中負號是考慮了速度的方向）
v（t+dt）=v（t）¡gdt（1.11）
在t+dt時刻小球的位移為
y（t+dt）=y（t）+v（t）dt（1.12）
在小球撞擊水平面的瞬間，即y（t）=0的時刻，它的速度方向改變，大小按比例K衰減。
當K=1時，就是無損耗彈跳情況。
因此，小球反彈瞬間（t+dt時刻）的速度為
v（t+dt）=-Kv（t）-gdt;0反彈瞬間的位移為
y（t+dt）=y（t）-Kv（t）dt=-Kv（t）dt
（2）仿真模型設計（程序）
從數學模型中可見，小球在空中自由運動時刻與撞擊時刻的動力方程不同。
通過小球所
處位置（位移）是否為零可判定小球處於何種狀態。
程序中採用if語句來作出判斷，以決定使
用（1.11）還是使用（1.13）來計算。
程序文件代碼如下。
\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:
N
if（y>0）|（v>0）%小球在空中的（含剛剛彈起瞬間）動力方程計算
v=v-g*dt;
y=y+v*dt;
else%碰擊瞬間的計算
y=y-K.*v*dt;
v=-K.*v-g*dt;
end
s（k）=y;%將當前位移記錄到s數組中以便作圖
end
t=t0:
dt:
dt*（N-1）;%仿真時間長度
plot（t,s）;
xlabel（'时间t'）;
ylabel（'位移y（t）'）;
axis（[0501.1]）;
holdon;
......
legend（'K=0.85','K=1'）;
圖1.3中分別作出了是碰擊衰減係數K=1和K=0:
85兩種情況下的小球彈跳位移曲
線。
對程序稍加修改就可以得到顯示小球彈跳過程的動畫。
修改後的程序文件代碼如下，讀者可運行該程序觀察不同的碰擊衰減係數下的小球彈跳過程。
〔程序代碼〕ch1example2prg2.m
%ch1example2prg2.m
g=9.8;%重力加速度
v0=0;%初始速度
y0=1;%初始位置
m=1;%小球質量
t0=0;%起始時間
K=0.85;%彈跳的損耗係數
N=5000;%仿真的總步進數
dt=0.005;%仿真步長
v=v0;%初狀態
y=y0;
fork=1:
N
ify>0%小球在空中的動力方程計算
v=v-g*dt;
y=y+v*dt;
else%碰擊瞬間的計算
y=-K.*v*dt;
v=-K.*v-g*dt;
end
plot（0,y,'o'）;
axis（[-2201]）;%坐標範圍固定
set（gcf,'DoubleBuffer','on'）;%雙緩衝避免作圖閃爍
drawnow;
end
用Matlab常微分方程求解器重新仿真實例1.2的乒乓球彈跳模型，要求也能
夠在仿真過程中動態顯示球的墜落和反彈過程，在仿真結束後輸出小球位移和速度的變化曲
線。
設小球的初始位置距離水平面1米，由靜止狀態開始自由墜落，並設反彈瞬間的速度衰減
係數為K=0:
85。
數學模型與實例1.2完全相同，為了利用標準求解器進行求解，首先將模型中的方程改寫
為標準的狀態方程形式，小球在空中時，有運動方程
functionxdot=ch2example6statefun（t,x,flag）
%乒乓球彈跳模型的標準狀態方程
%x
（1）為小球速度，x
（2）為小球位移
xdot=zeros（2,1）;%狀態變量矩陣初始化
xdot
（1）=-9.8;%速度加速度方程
xdot
（2）=x
（1）;%位移速度方程
並注意小球碰撞水平面瞬間（在反彈之前）的條件是y60且v60。
碰撞瞬間速度發生反向
並以係數K衰減。
據此編寫程序，程序中以x1（t）表示速度狀態變量v（t），以x2（t）表示位移
狀態變量y（t），並以矩陣形式描述，代碼如下。
v0=0;y0=1;%球的初始狀態
x_state=[v0,y0];%將初始狀態賦值到狀態變量中
dt=0.01;%仿真步進
t=0:
dt:
5;%仿真時間序列
K=0.85;%碰撞衰減係數
fork=1:
length（t）%仿真開始，每次循環向前推進一個仿真步進dt
x（k,:
）=x_state;%記錄並保存當前狀態的計算結果
[t_out,x_out]=ode45（'ch2example6statefun',[t（k）,t（k）+dt],x_state）;
%計算下一個時刻的新狀態
%可換用ode23、ode113、ode23t、ode15s、ode23s以及ode23tb求解器
x_state=x_out（length（x_out）,:
）;%更新狀態
if（x_state
（2）<=0）&（x_state
（1）<0）%當速度為負（球向下運動）且已經接觸碰撞面
x_state
（1）=-K*x_state
（1）;%處理碰撞瞬間情況：
速度反向並衰減K
end
%動畫作圖：
顯示小球彈跳過程
y=x_state
（2）;
subplot（2,1,1）;
legend（'小球运动过程'）;
axis（[-2201]）;%坐標範圍固定
set（gcf,'DoubleBuffer','on'）;%雙緩衝避免作圖閃爍
drawnow;%立即顯示作圖
end
%仿真結束。
最後輸出計算時間序列上的結果：
隨時間變化的速度和位移曲線
subplot（2,1,2）;
[AX,H1,H2]=plotyy（t,x（:
1）,t,x（:
2）,'plot'）;
set（AX
（1）,'Ylim',[-5,5]）;
set（AX
（2）,'Ylim',[0,1]）;
set（AX
（1）,'yTick',[-5:
5:
5]）;
set（AX
（2）,'yTick',[0:
0.5:
1]）
legend（'小球瞬时速度反弹瞬间跃变','小球瞬时位移'）;
xlabel（'时间（秒）'）;
set（get（AX
（1）,'ylabel'）,'string','速度'）;
set（get（AX
（2）,'ylabel'）,'string','位移'）;
〔實例1.4〕實際物理試驗中，當我們讓一個乒乓球垂直下落到一個完全水平的玻璃板上
後，乒乓球不斷彈跳，直到能量耗盡。
假定空氣是靜止的，沒有風，但是我們仍然會觀察到彈跳中的乒乓球在玻璃板上的落點不會是同一點。
這說明在乒乓球運動過程中受到微弱的水平面方向力的作用，產生了水平面方向上的漂移。
這些水平力在實例1.2中被我們忽略了，所以那裡仿真的結果中小球落點總是在坐標原點處。
如果我們要建立更加接近真實的物理環境的彈跳模型，就必須考慮這些被忽視的微小的擾動因素。
通過物理實驗觀察，我們可以做這樣的合理假設：
水平面方向上對乒乓球的微弱作用力可能來自多種因素的綜合，其中各因素對合力的貢獻甚小。
根據大數定理，在數學上我們就可以將水平作用力建模為一個高斯隨機變量。
為簡單起見，這裡仍然忽略了空氣對小球的其他作用因素，如球運動中的阻力，空氣的浮力等等。
同時，我們將實例1.2推廣到三維空間中的情況。
大数定律又称大数法则、大数率。
在一个随机事件中，随着试验次数的增加，事件发生的频率趋于一个稳定值；同时，在对物理量的测量实践中，大量测定值的算术平均也具有稳定性。
在数理统计中，一般有三个定理，贝努利定理和辛钦定理，如：
反映算术平均值和频率的稳定性。
当n很大时，算术平均值接近数学期望；频率以概率收敛于事件的概率
設水平面為xz坐標平面，y軸指向為垂直方向。
小球在x方向上的受力Fx（t）是一個零均值獨立高斯隨機過程。
小球在z方向上的受力Fz（t）與Fx（t）具有相同的分佈，但兩者相互獨立。
即
x、z方向相應的加速度、速度和位移分別用
表示，小球的質量為m。
由牛頓第二運動定律，可得出以下運動方程。
z方向的運動方程類似。
據此編寫仿真程序代碼如下。
M文件
g=9.8;%重力加速度
v0=0;%初始速度
y0=1;%初始位置
m=0.1;%小球質量
t0=0;%起始時間
K=0.8;%彈跳的損耗係數
N=5000;%仿真的總步進數
dt=0.005;%仿真步長
v=v0;%初狀態
y=y0;
vx=0;
vz=0;
sx=0;
sz=0;
fork=1:
N
ify>0%小球在空中的動力方程計算
v=v-g*dt;
y=y+v*dt;
else%碰擊瞬間的計算
y=-K.*v*dt;
v=-K.*v-g*dt;
End
Fx=randn;%x水平方向的隨機力,方差為1
ax=Fx./m;%Fx導致的x水平方向的加速度
vx=vx+ax*dt;%小球在x水平方向的瞬時速度
sx=sx+vx*dt;%小球在x水平方向的位移
Fz=randn;%z水平方向的隨機力,方差為1
az=Fz./m;%Fz導致的z水平方向的加速度
vz=vz+az*dt;%小球在z水平方向的瞬時速度
sz=sz+vz*dt;%小球在z水平方向的位移
plot3（sx,sz,y,'.'）;
gridon;
holdon;
axis（[-22-2201]）;%坐標範圍固定
set（gcf,'DoubleBuffer','on'）;%雙緩衝避免作圖閃爍
xlabel（'水平方向x'）;
ylabel（'水平方向z'）;
zlabel（'垂直方向y'）;
drawnow;%作圖
End
仿真以動畫方式進行，以便於觀察。
圖1.6是程序運行的結果。
圖中顯示了小球的運動軌
跡，彈跳的落點是隨機的。
修改小球的質量，將看到彈跳落點的概率特性也發生變化。
質量大的球落點相對集中。
讀者也可將空氣阻力考慮到數學模型中，從而仿真出根據真實的彈跳過程。
球类的自由落体运动是物理学中一个简单而又经典的模型，对于这个模型的充分研究，可以让我们更加深刻地体会物体运动的规律，了解物体运动的本质。
本文从简单的一维直线落体讲起，扩展到球体在二维、三维空间的运动模型，借助matlab对该运动系统进行仿真。
系统介绍
不受任何阻力，只在重力作用下而降落的物体，叫“自由落体”。
如在地球引力作用下由静止状态开始下落的物体。
地球表面附近的上空可看作是恒定的重力场。
如不考虑大气阻力，在该区域内的自由落体运动是匀加速直线运动。
其加速度恒等于重力加速度g。
自由落体运动的特点，体现在“自由”二字上，其含意为：
物体开始下落时是静止的即v0＝0。
如果给物体一个初速度竖直下落，不能算自由落体。
物体在下落过程中，除受重力作用外，不再受其他任何作用力（包括空气阻力）。
本文利用matlab对小球作自由落体运动的细节作出了仔细的研究分析。
结论
在这篇文章里，我们利用matlab软件实现了对于球体自由落体运动系统仿真分析。
首先，我们从一维直线运动开始（模型1），对物体下落的理想状态进行了建模，并由仿真图象分析出步长对于仿真结果的影响：
步长越短，仿真越逼真；其次，我们研究了二维平面中小球掉落以及其发生弹性和非弹性碰撞的不同情形（模型2），由仿真图象分析出在实际情况下（K），小球在多次反弹后最终将因为能量耗尽而停止运动，并且知道在碰撞瞬间小球速度大小不变，方向相反；最后，在考虑了水平扰动力作用下，我们研究了小球在三维空间的运动情况，利用随机函数对小球的运动进行了预测，并得出了小球质量越大，落点相对越集中地结论。
升级会员 