加農炮彈運動軌跡的研究

05-18

加農炮彈運動軌跡的研究

來自專欄計算物理學習筆記

摘要

加農炮在軍事上有著廣泛的應用，研究加農炮彈出膛後的運動軌跡具有重要意義。炮彈運動的實質是拋體運動。拋體運動是我們熟知的經典問題，在不計及空氣阻力的情況下，其運動方程十分簡單，可以解析求解，但在空氣阻力不可忽略的情況下，拋體運動方程將變得複雜，一般只能數值求解。本文利用Euler方法和Runge-Kutta方法對忽略空氣阻力和計及空氣阻力兩種情況下加農炮彈的運動軌跡進行了求解，研究了炮彈在不同發射角下的炮彈運動軌跡，考察了數值計算中時間步長對計算結果的影響。我們還研究了炮彈水平射程對發射角的依賴關係，並估計了達到最大水平射程時的發射角。進一步，我們研究了在考慮阻力情況下的中靶問題，給出了對空中目標實現精確打擊所需發射角的計算。

研究背景

拋體運動

加農炮彈出膛後的運動本質上是拋體運動，為便於處理，忽略炮彈在飛行過程中的轉動和形變，將其視為質點。質點的拋體運動是我們所熟知的，可以利用牛頓力學對其進行描述。在不計及空氣阻力的情況下，質點只受到重力作用，因而其運動方程由牛頓運動定律給出

$egin{equation} mfrac{mathrm{d}^{2}vec r}{mathrm{d} t^{2}}=vec F_{g} end{equation} ag{1}$

其中 $m$ 為炮彈的質量， $vec F_{g}$ 為重力。將上式寫成分量形式

$egin{equation} egin{cases} displaystylefrac{mathrm{d}^{2}x }{mathrm{d} t^{2}}=0\ displaystylefrac{mathrm{d}^{2}y}{mathrm{d} t^{2}}=-g end{cases} end{equation} ag{2}$

其中 $g$ 為重力加速度。這是一組容易解析求解的微分方程，若給定初始條件，炮彈的運動情況可以完全地確定。但在考慮空氣阻力的情況下，上述方程右側需要增加與阻力相關的項，這將使得炮彈的實際運動情況變得複雜。可以將炮彈的運動方程一般性地寫成如下形式

$egin{equation} egin{cases} displaystylefrac{mathrm{d}^{2}x }{mathrm{d} t^{2}}=f_{x}(x,y,v_{x},x_{y};t)\ displaystylefrac{mathrm{d}^{2}y}{mathrm{d} t^{2}}=-g+f_{y}(x,y,v_{x},x_{y};t) end{cases} end{equation} ag{3}$

其中 $f_{x}(x,y,v_{x},x_{y};t), f_{y}(x,y,v_{x},x_{y};t)$ 分別為單位質量的物體在 $x$ 方向和 $y$ 方向受到的除重力以外的力， $f_{x}$ 和 $f_{y}$ 的形式與具體情況有關。我們僅考察空氣阻力，並考慮簡單的阻力形式

$egin{equation} vec F_{mathrm{drag}}=-B_{1}vec v-B_{2}|v|vec v end{equation} ag{4}$

其中 $B_{1}, B_{2}$ 為阻力係數，是不顯著依賴於物體性質和狀態的常數。在許多情況下，僅考慮空氣阻力對速度二次方的依賴項，從而簡化模型，得到描述炮彈運動的方程

$egin{equation} egin{cases} displaystylefrac{mathrm{d}^{2}x }{mathrm{d} t^{2}}=-frac{B_{2}vv_{x}}{m}\ displaystylefrac{mathrm{d}^{2}y}{mathrm{d} t^{2}}=-g-frac{B_{2}vv_{y}}{m} end{cases} end{equation} ag{5}$

上述方程結合一定的初始條件可以完全確定炮彈的運動狀態，但方程 $(5)$ 的解析求解是困難的，因此我們將採用數值方法給出其數值解。

拋體運動方程的數值解法

Euler方法

我們利用Euler方法來求解上面的方程，將兩個二階微分方程降階，得到兩組一階微分方程組，按照Euler法的思想，將微分近似為差分，這裡採用向前差分，得到方程 $(5)$ 的顯式Euler法計算格式

$egin{equation} left { egin{aligned} x_{i+1}&=x_{i}+v_{x,i}Delta t\ v_{x,i+1}&=v_{x,i}-frac{B_{2}vv_{x,i}}{m}Delta t\ y_{i+1}&=y_{i}+v_{y,i}Delta t\ v_{y,i+1}&=v_{y,i}-g-frac{B_{2}vv_{y,i}}{m}Delta t end{aligned} ight. end{equation} ag{6}$

特別地，當 $B_{2}=0$ 時，上式退化為忽略空氣阻力的情形。由上述計算格式可以看出，只要給定初始條件 $x_{0},v_{x,0},y_{0}$ 和 $v_{y,0}$ ，即可逐步遞推後續各時刻拋體的位置、速度，也就是說，只要知道前一時刻 $t_{i}=iDelta t$ 拋體的位置、速度，即可求得下一時刻 $t_{i+1}=(i+1)Delta t$ 拋體的位置、速度，如此反覆迭代，最終可以得到一系列時刻拋體的運動情況。由於微分方程的數值解法是存在誤差的，上述顯式Euler法的整體截斷誤差與時間步長 $Delta t$ 相關，因此適當縮小步長可以提高計算精度，但這往往會以增加計算時間作為代價，因此需要權衡兩者的利弊，選定合適的時間步長進行數值求解，後面我們將會考察步長對計算結果的影響。

Runge-Kutta方法

作為對比，我們採用四階Runge-Kutta方法對方程 $(5)$ 進行求解，同樣地，先將二階微分方程降階，得到一階微分方程組

$egin{equation} egin{split} frac{mathrm{d} ilde{x}_{m}}{mathrm{d} t}&=f_{m}(x,y,v_{x},v_{y})qquad m=1,2,3,4\ ilde{x}_{1}&=x, ilde{x}_{2}=y, ilde{x}_{3}=v_{x}, ilde{x}_{4}=v_{y} end{split} end{equation} ag{7}$

為便於表述，上面並未寫出微分方程組的具體形式。按照Runge-Kutta法的基本思想，得到計算格式

$egin{equation} ilde{x}_{m,i+1}= ilde{x}_{m,i}+frac{Delta t}{6}(k_{m,1}+2k_{m,2}+2k_{m,3}+k_{m,4})\ end{equation} ag{8}$

$egin{equation} egin{cases} displaystyle k_{m,1}=f_{m}( ilde{x}_{1,i}, ilde{x}_{2,i}, ilde{x}_{3,i}, ilde{x}_{4,i})\ displaystyle k_{m,2}=f_{m}( ilde{x}_{1,i}+frac{Delta t}{2}cdot k_{1,1}, ilde{x}_{2,i}+frac{Delta t}{2}cdot k_{2,1}, ilde{x}_{3,i}+frac{Delta t}{2}cdot k_{3,1}, ilde{x}_{4,i}+frac{Delta t}{2}cdot k_{4,1})\ displaystyle k_{m,3}=f_{m}( ilde{x}_{1,i}+frac{Delta t}{2}cdot k_{1,2}, ilde{x}_{2,i}+frac{Delta t}{2}cdot k_{2,2}, ilde{x}_{3,i}+frac{Delta t}{2}cdot k_{3,2}, ilde{x}_{4,i}+frac{Delta t}{2}cdot k_{4,2})\ displaystyle k_{m,4}=f_{m}( ilde{x}_{1,i}+Delta tcdot k_{1,3}, ilde{x}_{2,i}+Delta tcdot k_{2,3}, ilde{x}_{3,i}+Delta tcdot k_{3,3}, ilde{x}_{4,i}+Delta tcdot k_{4,3}) end{cases} end{equation} ag{9}$

將 $f_{m}$ 的具體形式代入上式，即可得到具體的計算格式表達式，當 $B_{2}=0$ 時，上式退化為忽略空氣阻力的情形。同樣地，在給定初始條件 $x_{0},v_{x,0},y_{0}$ 和 $v_{y,0}$ 後，可以逐步遞推得到後續各時刻 $t_{i}=iDelta t$ 拋體的位置、速度，從而確定拋體的運動情況。四階Runge-Kutta法的整體截斷誤差與時間步長 $(Delta t)^{4}$ 相關，是比Euler法計算精度更高的一種方法，適當縮小時間步長同樣可以進一步提高該方法的精度，因此選擇合適的步長是很重要的，後面也將考察步長對計算結果的影響。

加農炮彈運動的研究

運動方程的解

在忽略空氣阻力的情況下，炮彈的運動方程如式 $(2)$ 所示。給定初始條件：炮彈的初始位置 $(x,y)=(0,0)$ ，初始速度 $v_{0}=500mathrm{m/s}$ ，通過分離變數法，可以得到解析解

$egin{equation} egin{cases} displaystyle x=v_{x0}t\ displaystyle y=v_{y0}t-frac{1}{2}gt^{2} end{cases} end{equation} ag{10}$

其中， $v_{x0},v_{y0}$ 是初始速度 $v_{0}$ 在 $x$ 方向和 $y$ 方向的分量。式 $(10)$ 表明炮彈的運動軌跡是一拋物線，這是我們所熟知的。對於方程 $(2)$ ，我們也可以採用數值方法進行求解，通過Matlab編程，得到計算結果如圖1所示。

在考慮空氣阻力的情況下，炮彈的運動狀態變得複雜，如前所述，解析地求解運動方程 $(5)$ 變得十分困難，這時適合採用數值計算的方法給出方程的數值解。這裡我們分別利用Euler法和Runge-Kutta法求出炮彈運動方程的數值解。選取發射角 $heta_{0}=45^{circ}$ ，繪製炮彈的運動軌跡如圖1所示。

圖1: 拋射角為45度時炮彈運動方程的解

從圖中可以看到，在忽略空氣阻力的情況下，數值解的各離散點較好地符合於式 $(10)$ 所給出的解析解，可以明顯地看出，Runge-Kutta法的精度要優於Euler法，前者的計算結果全部落在解析解的圖線上，而通過Euler法得到的解在自變數 $t$ 增大時對解析解存在小的偏離。總的來說，上述結果驗證了我們採用的數值方法的正確性和可行性，這為我們在考慮阻力情況下的數值計算提供了保障

現在我們對有無空氣阻力這兩種情況下炮彈運動狀態的差別作一簡要分析。在圖1中，我們選定了發射角為 $45^{circ}$ ，初始速度為 $500mathrm{m/s}$ 的情形進行考察。可以看到，在忽略阻力時，炮彈的拋體運動遵從式 $(10)$ 所給出的運動規律，其運動軌跡是一標準的拋物線，關於最高點確定的豎直線左右對稱；在考慮阻力時，炮彈的運動軌跡不再是標準的拋物線，沒有對稱的特性，並且炮彈的水平射程縮短，所能達到的最高點降低，體現出空氣阻力對拋體運動狀態的影響。關於這一問題，我們在後面還會詳細介紹。

如前所述，我們初步驗證了數值方法Euler法和Runge-Kutta法在求解微分方程時的可靠性，為了考察在不同發射角下炮彈的運動狀態以及進一步驗證上述兩種數值計算方法的正確性，我們選取了四個發射角（分別為 $35^{circ},45^{circ},55^{circ}$ 和 $65^{circ}$ ）對忽略空氣阻力情況下炮彈的運動軌跡進行了數值計算，並與該拋射角下的解析解進行了對比，通過Matlab編程，我們得到如圖2所示的結果。

圖2: 忽略空氣阻力時炮彈在不同拋射角下的運動狀態

由圖可知，通過Euler法和Runge-Kutta法得到的數值解都較好地符合於精確解，後者的精度更高。炮彈在忽略空氣時的運動軌跡呈現出標準的拋物線，且水平射程隨著拋射角的增大先逐漸增大，後逐漸減小，在拋射角 $heta_{0}=45^{circ}$ 時達到最大射程，這是我們熟知的結果。

值得說明的是，這裡我們選取 $35^{circ},45^{circ},55^{circ}$ 和 $65^{circ}$ 這四個角度僅作為代表，不具有任何特殊性。若希望考察其他拋射角度下炮彈的運動情況，可以運行Github上的程序，並手動輸入 $0^{circ}$ 到 $90^{circ}$ 之間的任意角度及其它參數，程序將給出該發射角下的運行結果。

數值解的穩定性與收斂性

我們在前面提到，利用Euler法或Runge-Kutta法求解微分方程得到的數值解是存在誤差的。根據兩種方法的計算格式，不難發現，解的誤差與計算時所採用的時間步長相關，因此步長會在一定程度上影響數值解的好壞。過長的步長可能導致解偏離真實情況，而過短的步長會顯著增加計算時間，因此選擇一個合適的步長對數值計算是極為重要的。

為了研究數值計算的結果與時間步長的關係，我們選取一個初始步長，進行計算，之後逐步減小步長取值，考察數值解對時間步長的依賴性，我們期待當步長減小到某一較小值時，計算結果將不顯著依賴於步長的選取。取步長分別為 $Delta t=20.0 mathrm{s},10.0 mathrm{s},5.0 mathrm{s},1.0 mathrm{s},0.5 mathrm{s},0.1 mathrm{s}$ ，並選定拋射角 $heta_{0}=45^{circ}$ 。利用Euler法和Runge-Kutt法求得考慮阻力情況下炮彈軌跡的數值解，如圖3所示。

圖3: 考慮空氣阻力情況下炮彈軌跡的數值解與時間步長的關係

從圖中可以發現，隨著時間步長的減小，炮彈運動軌跡的數值解發生形態上的變化，當步長取到 $Delta t=0.5 mathrm{s}$ 和 $0.1 mathrm{s}$ 時，在這兩個步長下計算得到的結果沒有顯著的差別。據此，我們做出如下論斷：

隨著時間步長的縮小，數值解的精度顯著提升，沒有出現振蕩的情況，因此可以相信解是穩定的；
隨著時間步長的縮小，Runge-Kutta法的計算結果較Euler法更快地收斂到高精度的數值解，體現出Runge-Kutta法的優越性；
當時間步長取到 $Delta t=0.5 mathrm{s}$ 時，繼續減小步長已經不能顯著提高解的精確性，但會增加計算耗費的時間。因此，我們斷定 $0.5 mathrm{s}$ 是一個合適的步長，在後面的計算中，我們將一律使用此步長。

空氣阻力對拋體運動的影響

在許多情況下，我們在研究拋體運動時認為物體只受到重力作用，這稱為理想拋體運動，而在實際情況下，物體還可能受到除重力以外的其他外力，情況要複雜得多。正如我們現在研究的加農炮彈出膛後的運動，炮彈在飛行過程中，除受到重力外，一般還受到空氣阻力的作用，其運動方程為

$egin{equation} mfrac{mathrm{d}^{2}vec r}{mathrm{d} t^{2}}=vec F_{g}+vec F_{mathrm{drag}} end{equation} ag{11}$

其中 $vec F_{mathrm{drag}}$ 為空氣阻力。在忽略阻力時，炮彈的運動規律由式 $(10)$ 描述，這是熟知的結果，但在考慮阻力時，炮彈的運動將遵從方程 $(11)$ ，其運動行為將由該方程的解給出。前面提到，方程 $(11)$ 的解析求解是困難的，其解也可能是複雜的，一般地寫為

$egin{equation} vec r=vec r(t) end{equation} ag{12}$

考慮到上述原因，我們利用數值方法對其進行研究。

為了對比忽略空氣阻力和考慮空氣阻力兩種情況下炮彈的運動狀態，我們將數值方法得到的解繪製於同一幅圖中，並選取了多個拋射角。通過Matlab編程，得到結果如圖4所示。

圖4: 有無空氣阻力時炮彈在不同拋射角下的運動軌跡

從圖中可以看到，空氣阻力的存在會顯著影響炮彈的運動狀態。當發射角取不同值時，炮彈的軌跡雖有所差別，但都具有一致的本質特徵，即相較於無阻力的情形，空氣阻力的存在使得炮彈的運動軌跡不再是標準的拋物線，水平射程和最高點都有大幅度的降低，直觀地看，空氣阻力將使炮彈的水平射程大約減小為無阻力時的一半，因此在實際應用時，要提高炮彈的射程，一個有效的辦法就是減小空氣阻力，這可以通過對炮彈形狀作出適當的設計而做到。此外，空氣阻力還將影響炮彈達到最大射程時的發射角，後面的計算將表明，當空氣阻力存在時，炮彈的最大射程並不是在發射角 $heta_{0}=45^{circ}$ 時取到，這與無阻力時的情形不同。這意味著空氣阻力顯著影響了炮彈射程對發射角的依賴關係。

拋體運動的最大射程

我們知道，拋體在同一發射速度情況下，所能達到的最大射程與拋射角密切相關。由於我們給定了炮彈的初始速度 $v_{0}=500 mathrm{m/s}$ ，若要得到最大射程，則需要研究炮彈運動軌跡與拋射角的關係。在忽略空氣阻力的情況下，拋體在拋射角 $heta_{0}=45^{circ}$ 時達到最大射程，這是我們熟知的結果；在考慮空氣阻力的情況下，拋體受到與其速度二次方成正比的阻力作用，其最大射程將與發射角呈現非線性關係，使得求極值的過程變得非常困難，難以用顯式形式給出最大射程及對應的射角。如前所述，炮彈在受到空氣阻力作用時的運動方程為

$egin{equation} left { egin{aligned} mfrac{mathrm{d}^{2}x }{mathrm{d} t^{2}}&=-B_{2}sqrt{v_{x}^{2}+v_{y}^{2}}cdot v_{x}\ mfrac{mathrm{d}^{2}y}{mathrm{d} t^{2}}&=-B_{2}sqrt{v_{x}^{2}+v_{y}^{2}}cdot v_{y}-mg end{aligned} ight. ag{13} end{equation}$

其中 $B_{2}$ 為阻尼係數。方程組 $(13)$ 為強耦合方程組，一般而言不能求出軌跡方程的顯式形式，需要採用數值方法求解。為了給出考慮空氣阻力情況下炮彈的最大射程及對應的拋射角，我們將在得到數值解的基礎上，近似求得炮彈的落地點，再通過搜索法給出達到最大水平射程時的拋射角。

由於數值計算的離散性，一般來說我們不能直接求得拋體高度剛好為零的點（即落地點），因此我們可以考慮採用近似的方法。正如前面提到的，在利用Euler法或Runge-Kutta法進行數值計算時，只要知道前一時刻 $t_{i}=iDelta t$ 拋體的位置 $x_{i}, y_{i}$ ，即可求得下一時刻 $t_{i+1}=(i+1)Delta t$ 拋體的位置 $x_{i+1}, y_{i+1}$ ，如此向後遞推，直至某一時刻 $t_{n+1}=(n+1)Delta t$ 的 $y$ 方向高度 $y_{n+1}<0$ 為止，這意味著拋體已經落地，無需再計算下去。根據上述分析，一定有

$egin{equation} y_{n}geq 0 ,qquad y_{n+1}<0 end{equation} ag{14}$

式 $(14)$ 表明落地點必定出現在 $(x_{n},y_{n})$ 和 $(x_{n+1},y_{n+1})$ 之間的軌跡上。我們可以利用落地點附近的數值解進行插值，求得插值多項式函數，進而通過求該函數的零點給出炮彈落地點的近似坐標。

考慮到線性插值的精度不高，我們採用Lagrange二次插值。選取落地點附近的三個數值解 $(x_{n-1},y_{n-1})$ ， $(x_{n},y_{n})$ 和 $(x_{n+1},y_{n+1})$ ，構造Lagrange基函數

$egin{equation} egin{aligned} cal L_{0}(x)&=frac{(x-x_{n})(x-x_{n+1})}{(x_{n-1}-x_{n})(x_{n-1}-x_{n+1})}\ cal L_{1}(x)&=frac{(x-x_{n-1})(x-x_{n+1})}{(x_{n}-x_{n-1})(x_{n}-x_{n+1})}\ cal L_{2}(x)&=frac{(x-x_{n-1})(x-x_{n})}{(x_{n+1}-x_{n-1})(x_{n+1}-x_{n})} end{aligned} end{equation} ag{15}$

進而得二次插值多項式函數

$egin{equation} y=y_{n-1}cal L_{0}(x)+y_{n}cal L_{1}(x)+y_{n+1}cal L_{2}(x) end{equation} ag{16}$

通過求該函數的零點即可得到落地點的近似坐標。我們利用上述方法從 $heta_{0}=1^{circ}$ 開始逐漸改變拋射角（角增量取為 $Delta heta_{0}=0.1^{circ}$ ）並計算各拋射角下落地點的橫坐標（即射程），順序搜索，直至某次計算得到的射程小於上一次的射程，結束搜索，得到最大射程及對應的拋射角。根據這一思想，計算得到達到最大射程時炮彈的軌跡如圖5所示，最大射程和對應的拋射角已標註於圖中。

圖5: 考慮阻力情況下炮彈達到最大射程時的運動軌跡

從圖中可以看出，在經Euler法和Runge-Kutta法算得的數值解的基礎上得到的最大射程和對應的拋射角存在一些微小的差別，在誤差的允許範圍內，我們姑且認為Runge-Kutta法得到的結果更接近真實情況，因此在考慮空氣阻力的情況下，當炮彈的拋射角 $heta_{0}=40.9^{circ}$ 時達到最大射程 $15031mathrm{m}$ 。利用同樣的方法，我們還可以算得忽略空氣阻力時的最大射程為 $25510mathrm{m}$ ，可見空氣阻力的存在使得最大射程減小約一半左右。這一結果對實際情況具有指導性的意義，由於炮彈在飛行過程中受到阻力，而阻力會顯著影響射程對拋射角的依賴關係，因此在實際中需要根據空氣阻力的具體條件來確定發射角以使炮彈能飛行更遠的距離。

炮彈對目標的精確打擊

在實際應用中，人們利用炮彈打擊的目標往往不是固定於地面的死靶，目標的豎直高度可能是任意的，因此不能只考慮打擊同高度的地面目標的情況，而是希望對不同高度的目標予以精確打擊。我們在這裡考察一種簡單的情況，即炮彈的發射速度固定為 $v_{0}=500 mathrm{m/s}$ ，加農炮處於坐標系原點，假設目標位於 $(15mathrm{km},3mathrm{km})$ 。下面我們將通過數值計算給出炮彈以何種發射角出射才能命中這一目標。

考慮存在空氣阻力的情況，利用Euler法和Runge-Kutta法求解運動方程 $(5)$ 。當發射角較大時，炮彈將能夠打到 $x=15mathrm{km}$ 附近的位置，由於數值計算的離散性，一般來說我們不能直接得到 $x$ 剛好為 $15mathrm{km}$ 的點（即目標點的 $x$ 方向坐標），仿照前面求射程時的辦法，我們採用Lagrange線性插值得到炮彈打到 $x=15mathrm{km}$ 時對應的豎直高度的近似值。利用上述方法從 $heta_{0}=1^{circ}$ 開始逐漸改變拋射角（角增量取為 $Delta heta_{0}=0.01^{circ}$ ）並計算各拋射角下炮彈打到 $x=15mathrm{km}$ 時對應的豎直高度和目標點的 $y$ 方向坐標的偏差 $Delta y$ ，如此順序搜索，直至某次計算得到的偏差 $Delta y$ 的值大於上一次的偏差值，結束搜索，得到擊中位於 $(15mathrm{km},3mathrm{km})$ 的目標所需的發射角。根據上述思想，得到炮彈對目標實現精確打擊時的運動軌跡如圖6所示，命中目標所需的發射角以及炮彈打擊點與目標點的豎直方向偏差值 $Delta y$ 已標註於圖中。

圖6: 考慮阻力情況下炮彈在恆定發射速度時實現對目標的精確打擊

我們在計算時採取時間步長 $Delta t=0.3mathrm{s}$ 以提高計算精度。可以發現，兩種數值方法給出基本一致的結果，如果要打擊位於 $(15mathrm{km},3mathrm{km})$ 的目標，應選取發射角 $heta_{0}=34^{circ}$ 左右。值得說明的是，我們在這裡把約化空氣阻力係數取為 $1.0 imes 10^{-5}mathrm{m^{-1}}$ 而不是前面一直採用的 $4.0 imes 10^{-5}mathrm{m^{-1}}$ ，我們發現在空氣阻力過大時，無論以何種發射角出射，炮彈都不能擊中目標。這是因為過大的空氣阻力過多地耗散了炮彈的能量，只有增大炮彈的發射速度才有可能擊中目標。由於我們固定了發射速度為 $500 mathrm{m/s}$ ，那麼在這一恆定初始速度下，如果空氣阻力增大到一定值時，無論發射角為何，炮彈都將無法命中目標。

結論

我們利用微分方程的數值解法Euler法和Runge-Kutta法研究了加農炮彈出膛後的運動軌跡，考察了忽略空氣阻力和考慮空氣阻力兩種情況下炮彈的運動狀態。我們的研究表明，空氣阻力的存在會顯著影響炮彈運動軌跡的形態，使其水平射程和最高點大幅度降低，空氣阻力還會影響炮彈的射程對拋射角的依賴關係，使其達到最大射程時對應的拋射角不同於無阻力的情形。此外，我們還考察了數值解的收斂性問題，研究了數值計算結果與時間步長的關係。我們指出，Euler法和Runge-Kutta法在處理炮彈運動問題時得到的解是收斂的和穩定的，這保證了我們前面作出的一系列結論的可靠性和正確性。最後我們研究了炮彈對目標的打擊問題，給出了一個計算例子。我們期望本文得到的研究結果能對炮彈發射相關的實際問題提供一定的幫助。

參考文獻

[1] 劉金遠, 段萍, 鄂鵬. 計算物理學[M]. 北京: 科學出版社, 2012.

[2] 郝亞非. 拋體運動軌跡的數值分析[J]. 物理通報. 2015(12): 10-12.