基于Mathematica的阻尼拋體數(shù)值仿真實(shí)驗(yàn)

胡宗元

（首都經(jīng)濟(jì)貿(mào)易大學(xué) 安全與環(huán)境工程學(xué)院，北京 100070）

·計(jì)算機(jī)科學(xué)與實(shí)驗(yàn)·

基于Mathematica的阻尼拋體數(shù)值仿真實(shí)驗(yàn)

胡宗元

（首都經(jīng)濟(jì)貿(mào)易大學(xué) 安全與環(huán)境工程學(xué)院，北京 100070）

該文基于Mathematica設(shè)計(jì)了一個(gè)能夠進(jìn)行阻尼拋物運(yùn)動(dòng)仿真試驗(yàn)的程序。在任意給定的輸入條件下，該程序既可以輸出運(yùn)動(dòng)實(shí)況動(dòng)畫和頻閃圖等仿真結(jié)果，又可以輸出阻尼拋物運(yùn)動(dòng)的軌跡曲線。

數(shù)值仿真；阻尼拋體；Mathematica；實(shí)驗(yàn)

10.3969/j.issn.1672-4550.2016.05.016

運(yùn)用計(jì)算機(jī)來進(jìn)行探究和設(shè)計(jì)是大學(xué)生應(yīng)具備的一項(xiàng)基本功。我國(guó)許多非計(jì)算機(jī)專業(yè)的大學(xué)生的計(jì)算機(jī)應(yīng)用能力目前還處于較低水平。在教學(xué)中開展數(shù)值仿真實(shí)驗(yàn)，不僅可以幫助學(xué)生更好地掌握專業(yè)知識(shí)，同時(shí)可以培養(yǎng)學(xué)生的計(jì)算機(jī)應(yīng)用能力。

常用的科學(xué)計(jì)算軟件如Matlab、Maple和Mathematica都可以選作設(shè)計(jì)數(shù)值仿真實(shí)驗(yàn)的工具，它們各有特色［1］。Mathematica是美國(guó)Wolfram Research公司開發(fā)的數(shù)學(xué)軟件，主要用途是科學(xué)研究與工程技術(shù)中的計(jì)算［1］。它簡(jiǎn)單易學(xué)、高效多能，一句簡(jiǎn)單的Mathematica命令就能完成一些編程軟件需用幾十句才能完成的工作。該軟件除了能進(jìn)行任意精度的數(shù)值計(jì)算，還可以進(jìn)行微分、積分、冪級(jí)數(shù)展開等符號(hào)演算。此外，Mathematica還提供了多種方式作圖、發(fā)聲和動(dòng)畫的命令，使人容易直觀地把握住對(duì)象的特性［3］?；谶@些特點(diǎn)，本文選擇Mathematica作為開發(fā)仿真實(shí)驗(yàn)的工具。

下面以阻尼拋體運(yùn)動(dòng)實(shí)驗(yàn)為例，具體展示用Mathematica設(shè)計(jì)數(shù)值仿真實(shí)驗(yàn)的方法和結(jié)果。

1 物理模型與求解

為方便起見，假設(shè)拋體的質(zhì)量為m，可以看成質(zhì)點(diǎn)，在重力和空氣阻力作用下運(yùn)動(dòng)。拋體受到的重力為mg，空氣阻尼力為mb，b為阻尼系數(shù)，取拋出點(diǎn)的位置為原點(diǎn)建立坐標(biāo)系，x軸為水平方向，y軸為豎直向上。在初始時(shí)刻t=0，物體以速度v0，仰角θ拋出，則由牛頓第二定律［4］得到運(yùn)動(dòng)微分方程和初始條件:

由式（1）可解出速度方程:

和運(yùn)動(dòng)方程:

式中，v∞=-g/b為垂直方向的收尾速度。

由此可以進(jìn)一步求出拋體所能達(dá)到的最大高度為:

落地時(shí)間tmax滿足:

對(duì)應(yīng)落地點(diǎn)的水平位置為:

當(dāng)阻尼系數(shù)足夠小，即b＜＜1時(shí)，有近似結(jié)果tmax≈2v0sinθ/g。

2 仿真程序設(shè)計(jì)

2.1 需求分析

按照仿真要求，阻尼拋體運(yùn)動(dòng)仿真實(shí)驗(yàn)需要得出的效果為阻尼拋體運(yùn)動(dòng)過程實(shí)況、具體記錄運(yùn)動(dòng)過程的頻閃圖和由此繪出的運(yùn)動(dòng)軌跡圖等。其中，運(yùn)動(dòng)過程實(shí)況可以通過Mathematica的動(dòng)畫功能來實(shí)現(xiàn)；頻閃圖和軌跡圖可以直接調(diào)用Mathematica命令里的Graphics和parametricPlot。

從控制論的角度看，一個(gè)具有特定功能的系統(tǒng)，需要有適當(dāng)?shù)妮斎牒洼敵觥Ｗ鳛榉抡鎸?shí)驗(yàn)程序，輸入的是控制變量，輸出的是實(shí)驗(yàn)結(jié)果。按照前面的分析，需要輸入的控制變量應(yīng)該有初速度v0、仰角θ和阻尼系數(shù)b；需要輸出的實(shí)驗(yàn)結(jié)果為阻尼拋體運(yùn)動(dòng)過程的動(dòng)畫顯示、過程的頻閃圖和運(yùn)動(dòng)的軌跡圖等。由于存在多種輸出，因此還需要增加一個(gè)用以選擇輸出模式的輸入變量——模式變量。

2.2 程序的設(shè)計(jì)

仿真程序可以認(rèn)為是將輸入信息加工為輸出的一種變換。利用輸入的初速度、仰角和阻尼系數(shù)，完全確定了運(yùn)動(dòng)方程式（3）。這部分內(nèi)容構(gòu)成了程序的基本模塊:

其中，引入兩個(gè)內(nèi)部變量g=9.8，vf=-g/b。

在此基礎(chǔ)上，調(diào)用Mathematica命令里的Show、Graphics和parametricPlot可以分別實(shí)現(xiàn)運(yùn)動(dòng)實(shí)況動(dòng)畫、頻閃圖和軌跡圖的輸出。

物理分析計(jì)算提供了計(jì)算機(jī)模擬的基礎(chǔ)，但是要獲得良好的演示效果，還需要注意展示形式。

1）首先是拋體的形狀和大小。由于不計(jì)形狀的影響，可以將拋體設(shè)為小球。從演示的效果看，球的半徑不宜過小，不然看不清楚；也不宜過大，不然占據(jù)過多的畫面難以看清運(yùn)動(dòng)過程。因此，適當(dāng)?shù)剡x擇是規(guī)定一個(gè)半徑的相對(duì)大小，如占整個(gè)畫面尺度的5。

2）其次是整個(gè)畫面的大小。為了充分利用屏幕空間，應(yīng)當(dāng)使拋體的運(yùn)動(dòng)范圍盡可能地充滿屏幕。按照第1節(jié)中的結(jié)果，拋體中心的運(yùn)動(dòng)范圍為x∈［0，xmax］，y∈［0，ymax］，因此展示范圍應(yīng)該包含上述區(qū)域。如果考慮到拋體的大小，運(yùn)動(dòng)范圍還需要適當(dāng)加大。

還要注意閃光間隔大小的選擇，按照實(shí)物實(shí)驗(yàn)的情況，取閃光間隔dt=0.2比較適當(dāng)。動(dòng)畫的時(shí)間間隔原則上可以選擇任意小，為了節(jié)約計(jì)算機(jī)資源，選擇閃光間隔的二分之一或者三分之一即可。

3 仿真程序及其運(yùn)行

3.1 Mathematica源程序

由上面的設(shè)計(jì)，不難編寫出相應(yīng)的Mathematica模塊程序如下:

dampingthrowing［v0_，q_，b_，m_Integer］:=

Module［｛g，ymax，t0，tmax，xmax，dt，vf｝，

g=9.8；dt=0.2；vf=-g/b；

x［t_］:=v0*cos［q］*（1-exp［-b* t］）/b；

y［t_］:=vf*t+（v0*sin［q］-vf）（1-exp［-b*t］）/b；

ymax=v0*sin［q］/b+vf*log［1-v0* sin［q］/vf］/b；

t0=2*v0*sin［q］/g；

tmax=t/.FindRoot［y［t］==0，｛t，t0｝］；

xmax=x［tmax］；

Switch［m，0， Table［Show［Graphics［｛PointSize［0.03］，Point［｛x［t］，y［t］｝］｝，PlotRange?｛｛-0.05xmax， 1.05xmax｝，｛-0.05ymax，1.05ymax｝｝］］，｛t，0，tmax，dt/3｝］，

1，Show［Table［Graphics［｛PointSize［0.03］，Point［｛x［t］，y［t］｝］｝，Axes? True，PlotRange? ｛-0.05ymax，1.05ymax｝］，｛t，0，tmax，dt｝］］，

2，ParametricPlot［｛x［t］，y［t］｝，｛t，0，tmax｝，PlotRange? ｛｛0，xmax｝，｛0，ymax｝｝］

］

］

3.2 程序中的參數(shù)說明

程序中各個(gè)參數(shù)說明如下。

輸入?yún)?shù):初速度v0，仰角θ，阻尼系數(shù)b，模式參數(shù)m。

輸出模式:m=0，動(dòng)畫；m=1，閃光照；m=2，運(yùn)動(dòng)軌跡圖。

內(nèi)部變量:重力加速度g=9.8；收尾速度vf；閃光間隔dt=0.2；最大高度ymax；最大時(shí)間tmax；最遠(yuǎn)距離xmax。

3.3 仿真程序的使用

在應(yīng)用時(shí)調(diào)用模塊dampingthrowing，輸入拋體的初速度v0；仰角θ；阻尼系數(shù)b；模式參數(shù)m等參數(shù)后，便可運(yùn)行。

例設(shè)拋體的初速度v0=20、仰角θ=π/4、阻尼參數(shù)b=2。

如果要得到動(dòng)畫過程，輸入命令為dampingthrowing［20，π/4，2，0］，輸出結(jié)果為阻尼拋體運(yùn)動(dòng)實(shí)況的仿真動(dòng)畫；由于動(dòng)畫效果在文章中無法展示，請(qǐng)讀者自行測(cè)試。

如果要得到頻閃圖，輸入命令為dampingthrowing［20，π/4，2，1］，輸出結(jié)果如圖1所示。

圖1 阻尼拋體運(yùn)動(dòng)的頻閃圖

如果要得到軌跡圖，輸入命令為dampingthrowing［20，π/4，2，2］，輸出結(jié)果如圖2所示。

圖2 阻尼拋體運(yùn)動(dòng)的軌跡圖

4 結(jié)束語

本文介紹了一個(gè)阻尼拋物運(yùn)動(dòng)仿真模塊的設(shè)計(jì)方案和運(yùn)行效果，為大學(xué)生在專業(yè)學(xué)習(xí)中使用計(jì)算機(jī)提供了一個(gè)可以模仿的案例。該程序主要有以下4個(gè)方面的優(yōu)點(diǎn):

1）輸入?yún)?shù)任意可調(diào)，可以按需要進(jìn)行各種模擬實(shí)驗(yàn)；

2）可以定量地演示阻尼拋物運(yùn)動(dòng)現(xiàn)象，仿真程度高；

3）既能夠輸出運(yùn)動(dòng)實(shí)況，又能夠輸出頻閃照片，演示效果好；

4）源程序簡(jiǎn)單，具有開放性，在該模塊的基礎(chǔ)上略加修改，就可以演示高處斜拋或者變阻尼的拋物運(yùn)動(dòng)等現(xiàn)象，還可以繪制相空間中的運(yùn)動(dòng)軌跡。

［1］彭芳麟，梁穎，忻蓓.計(jì)算軟件在計(jì)算物理課程中的地位和作用［J］.大學(xué)物理，2013（8）:7.

［2］洪維恩.數(shù)學(xué)運(yùn)算大師Mathematica4［M］.北京:人民郵電出版社，2002:1-16.

［3］倪致祥.科研的有力工具——Mathematica簡(jiǎn)介［J］.阜陽師范學(xué)院學(xué)報(bào)，2005（2）:7.

［4］周衍柏.理論力學(xué)教程［M］.2版.北京:高等教育出版社，1986:24-41.

Numerical Simulation Experiment of Damped Projectile Based on Mathematica

HU Zongyuan
（College of Safety and Environment Engineering，Capital University of Economics and Business，Beijing 100070，China）

We design a program to simulate the damping projectile experiment based on Mathematica.The program can output many simulation results such as animation，stroboscopic map and the track curve in coordinate space for given inputs.

numerical simulation；damped projectile；Mathematica；experiment

TP 391.6

A

2015-07-12；修改日期：2015-08-12

北京市自然科學(xué)基金（3144027）。

胡宗元（1984-），男，博士生，實(shí)驗(yàn)師，主要從事微波、天線與電磁環(huán)境方面的教學(xué)與研究。