胡志棟++董喜斌

摘要 針對(duì)學(xué)生缺少足夠的高等數(shù)學(xué)知識(shí)解決工程實(shí)踐中實(shí)際問(wèn)題的現(xiàn)狀，為了在森林工程課堂教學(xué)中將高等數(shù)學(xué)知識(shí)與工程實(shí)踐更好地結(jié)合，基于實(shí)際的研究工作，對(duì)用于解決常微分方程的顯式歐拉法、隱式歐拉法、改進(jìn)歐拉法、龍格庫(kù)塔法和ODE45方法進(jìn)行了介紹和對(duì)比研究。結(jié)果表明，這些方法在解決常微分方程的數(shù)值解上是一致的，而且相互之間的誤差很小。介紹了解決偏微分方程的PDE 方法。在建模的過(guò)程中用到了有限元方法，在方程的列寫過(guò)程中用到了矩陣方法，通過(guò)理論分析可以看出，矩陣表示的是實(shí)際問(wèn)題與空間、時(shí)間的關(guān)系，因而可以從直觀的角度去理解抽象的矩陣知識(shí)。這些方法可以很好地幫助森林工程領(lǐng)域的學(xué)生在今后的工作中建立數(shù)學(xué)模型并加以分析，從而加強(qiáng)解決實(shí)際工程問(wèn)題的能力，也為森林工程課題今后的教學(xué)提供一定的參考。

關(guān)鍵詞 高等數(shù)學(xué)；森林工程；教學(xué)課堂；常微分方程（ODE）；偏微分方程（PDE）；矩陣

中圖分類號(hào) G542 文獻(xiàn)標(biāo)識(shí)碼 A 文章編號(hào) 1007-5739（2017）11-0272-04

1 研究背景概述

縱觀世界數(shù)學(xué)發(fā)展史，17—19世紀(jì)的英國(guó)、德國(guó)、法國(guó)等歐洲大國(guó)都是數(shù)學(xué)強(qiáng)國(guó)。英國(guó)的牛頓提出了微積分理論，不僅在數(shù)學(xué)領(lǐng)域引發(fā)一場(chǎng)革命[1]，同時(shí)也大量應(yīng)用于許多物理問(wèn)題的研究。法國(guó)一直擁有深厚的數(shù)學(xué)文化發(fā)展歷史，而德國(guó)的哥廷根于20世紀(jì)初成為世界數(shù)學(xué)的中心[2]。俄羅斯在數(shù)學(xué)領(lǐng)域的發(fā)展從19世紀(jì)開始，至蘇聯(lián)時(shí)期成為世界數(shù)學(xué)強(qiáng)國(guó)之一。特別是蘇聯(lián)于1958年成功發(fā)射了第一顆人造地球衛(wèi)星，也象征著蘇聯(lián)在與此相關(guān)的數(shù)學(xué)領(lǐng)域處于世界領(lǐng)先地位。此外，蘇聯(lián)重視基礎(chǔ)科學(xué)教育，也是其在基礎(chǔ)科學(xué)研究中具有雄厚實(shí)力的一個(gè)重要原因。第二次世界大戰(zhàn)前美國(guó)在數(shù)學(xué)上遠(yuǎn)落后于歐洲，但如今也已成為數(shù)學(xué)超級(jí)大國(guó)。一方面由于大批猶太裔數(shù)學(xué)家被迫移居美國(guó)，增強(qiáng)了美國(guó)的數(shù)學(xué)實(shí)力[3]；另一方面在蘇聯(lián)發(fā)射第一顆人造地球衛(wèi)星后，美國(guó)加強(qiáng)了對(duì)數(shù)學(xué)研究和數(shù)學(xué)教育的投入，使其迅速成為一個(gè)數(shù)學(xué)強(qiáng)國(guó)。

學(xué)界普遍認(rèn)為下一次科技革命將以人類3種新的“生存形式”為重要標(biāo)志，即網(wǎng)絡(luò)人（生活在網(wǎng)絡(luò)空間的虛擬人）、仿生人（高仿真智能人）和再生人（具有自然人特征的“復(fù)制人”）。并預(yù)計(jì)這次科技革命將在2020—2050年到來(lái)。數(shù)學(xué)對(duì)整個(gè)社會(huì)經(jīng)濟(jì)的發(fā)展起著至關(guān)重要的作用[4-6]，在前幾次科技革命中都起到了先導(dǎo)與支柱的作用，因而有理由相信數(shù)學(xué)必將成為下一次科技革命重要的推動(dòng)力之一[7]。由此可見數(shù)學(xué)在科技發(fā)展中的重要地位[8]，本文對(duì)于將高等數(shù)學(xué)知識(shí)應(yīng)用到森林工程教學(xué)課堂進(jìn)行研究，以期推動(dòng)數(shù)學(xué)在教育領(lǐng)域的綜合應(yīng)用。

為了更好地說(shuō)明如何將高等數(shù)學(xué)知識(shí)應(yīng)用到森林工程教學(xué)課堂，文中采用了基于氣穴氣泡研究的一個(gè)實(shí)例。在過(guò)去的幾個(gè)世紀(jì)中，氣穴問(wèn)題在森林工程領(lǐng)域一直被認(rèn)為是流體機(jī)械氣蝕問(wèn)題的主要起因之一。大量的試驗(yàn)和數(shù)值分析工作都是圍繞企圖找出氣蝕問(wèn)題的產(chǎn)生機(jī)理展開的。Rayleigh[9] 最先開始這方面的工作，其研究認(rèn)為圓形氣泡潰破時(shí)產(chǎn)生的局部高壓，即氣穴是產(chǎn)生氣蝕現(xiàn)象的原因。式（1）為R-ayleigh推導(dǎo)得到的經(jīng)典氣穴氣泡動(dòng)力學(xué)方程。

r■+■■2=■（1）

式（1）中：r為氣泡半徑（m）；PB為氣泡內(nèi)的壓力（Pa）；P∞為外界無(wú)窮遠(yuǎn)處的壓力（Pa）； ρ為流體的密度（kg/m3）；t為時(shí)間（s）。

2 研究方法

2.1 量綱分析方法

國(guó)際單位組織定義了7個(gè)最基本的單位，其他的單位都可以由這7個(gè)基本單位導(dǎo)出。這7個(gè)基本單位分別為長(zhǎng)度，米（m）；質(zhì)量，千克（kg）；時(shí)間，秒（s）；電流，安培（A）；熱力學(xué)溫度，開爾文（K）；發(fā)光強(qiáng)度，坎德拉（cd）；物質(zhì)的量，摩爾 （mol）。長(zhǎng)度、質(zhì)量和時(shí)間這3個(gè)基本單位是森林工程課堂中常用的單位。

量綱分析法又稱為因次分析法，是一種數(shù)學(xué)分析方法，其可以正確地分析各變量之間的關(guān)系，簡(jiǎn)化試驗(yàn)和成果整理，所以量綱分析是分析流體運(yùn)動(dòng)的有力工具。通過(guò)量綱分析可以檢查反映物理現(xiàn)象規(guī)律的方程在計(jì)量方面是否正確，甚至可提供尋找物理現(xiàn)象某些規(guī)律的線索。

應(yīng)用量綱分析法，式（1）可以寫成：

m■+■■■=■∝■■（2）

通過(guò)這種分析方法，推導(dǎo)過(guò)程的準(zhǔn)確性可以很容易得到保證，而且方程的特性和實(shí)質(zhì)也很容易發(fā)現(xiàn)。

2.2 ODE求解方法

在數(shù)學(xué)領(lǐng)域中，只包含一個(gè)獨(dú)立變量及其導(dǎo)數(shù)的單個(gè)或多個(gè)微分方程被稱為常微分方程（ODEs）。常微分方程廣泛應(yīng)用于數(shù)學(xué)、工程和科學(xué)等領(lǐng)域。在數(shù)學(xué)中“變化”是用導(dǎo)數(shù)和微分來(lái)描述的。種類繁多的導(dǎo)數(shù)、微分及函數(shù)通過(guò)方程聯(lián)系在一起，用來(lái)刻畫各種動(dòng)態(tài)變化的現(xiàn)象、進(jìn)化或演變過(guò)程。

如果式（1）中的壓差為定值，那么它屬于一個(gè)ODE問(wèn)題，然后對(duì)比研究顯式歐拉法（Forward Euler′s method）、隱式歐拉法（Backward Euler′s method）、改進(jìn)歐拉法（Modified Euler′s method）、ODE45方法和龍格庫(kù)塔法（Runge-Kutta method）在分析ODE問(wèn)題上的應(yīng)用[10]。

顯式歐拉法表達(dá)式：

ri+1=ri+h×f（ti，ri）（3）

隱式歐拉法表達(dá)式：

ri+1=ri+h×f（ti+1，ri+1）（4）

改進(jìn)歐拉法表達(dá)式：

ri+1=ri+■×[f（ti+1，ri+1）+f（ti，ri）]（5）

4階龍格庫(kù)塔法表達(dá)式：

ri+1=ri+■×[K1+2K2+2K3+K4]（6）

式（6）中參數(shù)的計(jì)算如下：

K1=fti，riK2=fti+■，ri+■K1K3=fti+■，ri+■K2K4=fti+h，ri+hK3（7）

式（3）～（7）中：h為步長(zhǎng)。

ODE45方法是MATLAB程序軟件包內(nèi)植的一條命令程序。

首先應(yīng)用上述方法分析單個(gè)氣泡的動(dòng)態(tài)生長(zhǎng)過(guò)程，包括氣泡半徑變化情況和生長(zhǎng)速率情況（圖1）。

再應(yīng)用上述方法分析單個(gè)氣泡的動(dòng)態(tài)潰滅過(guò)程，包括氣泡半徑變化情況和潰滅速率情況（圖2）。

通過(guò)對(duì)比可以看出，這5種方法在解決單個(gè)氣泡的動(dòng)態(tài)問(wèn)題上是統(tǒng)一的，它們之間的誤差也很小。

2.3 矩陣方法

矩陣是數(shù)學(xué)中最重要的基本概念之一，是代數(shù)學(xué)的一個(gè)主要研究對(duì)象，也是數(shù)學(xué)研究及應(yīng)用的一個(gè)重要工具。在數(shù)學(xué)中，矩陣是一個(gè)按照長(zhǎng)方陣列排列的數(shù)、符號(hào)或表達(dá)式組成的集合。在不知道矩陣如何得來(lái)的情況下分析矩陣是很抽象的。為了更好地使學(xué)生理解矩陣的應(yīng)用，通過(guò)兩充液平行平面的動(dòng)態(tài)開啟過(guò)程實(shí)例來(lái)加以解釋。

在一定的初始間隙值下，建立兩充液平行平面的動(dòng)態(tài)開啟模型，然后通過(guò)理論分析可以得到開啟過(guò)程中壓力分布的變化情況（圖3）。

首先針對(duì)沒(méi)有氣穴氣泡的情況，為了分析壓力分布的變化情況，針對(duì)第i個(gè)單元，運(yùn)用間隙流動(dòng)的理論公式和有限元方法得到理論公式，最后用矩陣的方法得到整個(gè)問(wèn)題域的數(shù)學(xué)模型：

■×■=■（8）

式（8）中：h（t）為開啟高度函數(shù)；η為液體的動(dòng)力黏度（Pa·s）。

1 0 0 L 02 -3 1 O 0 1 -2 1 O M 0 1 -2 1 M O O O O 0 L 0 1 -2 10 L 0 1 -1 p■p■p■p■Mp■p■=a1p0/a1h′（t）h′（t）h′（t）Mh′（t）h（t）（9）

式（9）用向量表示可以簡(jiǎn)寫為：

■■=a■■（10）

式（10）中：a■=■，■=1 0 0 L 02 -3 1 O0 1 -2 1 O M 0 1 -2 1M O O O O 0 L 0 1 -2 10 L 0 1 -1，

■=P■P■P■P■MP■P■，■=p■/a■h′（t）h′（t）h′（t）Mh′（t）h（t）。

通過(guò)MATLAB編程求解，就可以得到壓力分布與空間、時(shí)間的關(guān)系（圖4）。

2.4 PDE求解方法

在數(shù)學(xué)領(lǐng)域，如果一個(gè)微分方程中出現(xiàn)多元函數(shù)的偏導(dǎo)數(shù)，或者說(shuō)如果未知函數(shù)和幾個(gè)變量有關(guān)，而且方程中出現(xiàn)未知函數(shù)對(duì)應(yīng)幾個(gè)變量的導(dǎo)數(shù)，那么這種微分方程就是偏微分方程。

許多問(wèn)題用一個(gè)自變量的函數(shù)來(lái)描述已不夠，不少問(wèn)題由多個(gè)變量的函數(shù)來(lái)描述。比如，從物理角度來(lái)說(shuō)，物理量有不同的性質(zhì)，溫度、密度等是用數(shù)值來(lái)描述的，叫做純量；速度、電場(chǎng)的引力等不僅在數(shù)值上有不同，而且還具有方向，這些量叫做向量；物體在一點(diǎn)上的張力狀態(tài)描述出的量叫做張量。這些量不僅和時(shí)間有關(guān)系，而且和空間坐標(biāo)也有聯(lián)系，這就要用多個(gè)變量的函數(shù)來(lái)表示。

同樣地，針對(duì)兩充液平行平面的動(dòng)態(tài)開啟過(guò)程，Matlab軟件提供了一個(gè)求解PDE問(wèn)題的工具箱“pdepe”?；趐depd方法的壓力分布見圖5。在求解的時(shí)候，初始條件IC（Initial Condition）和邊界條件BC（Boundary Condition）是需要提前確定的[11-13]。

Cx，t，u，■■=x■-m■xmfx，t，u，■

+Sx，t，u，■（11）

令 m=0， x？綴[0，0.012 5]，t？綴[0，0.002] ，fx，t，u，■=■，

Cx，t，u，■=0，Sx，t，u，■=-■■，

IC：p（x，0）=1bar，BC：p（0，t）=1barp（0.012 5，t）=1bar。

通過(guò)圖4和圖5的對(duì)比研究可以看出，PDE求解和理論分析的結(jié)果是一致的，這為以后的理論分析提供了另外一條途徑。

如果將兩充液平行平面的開啟過(guò)程和氣穴氣泡的動(dòng)態(tài)生長(zhǎng)過(guò)程結(jié)合起來(lái)，則可得出理論分析模型，具體如圖6所示。

同時(shí)考慮流體的黏度，則含有氣穴氣泡的兩充液平行平面動(dòng)態(tài)開啟過(guò)程可以用式（12）和式（13）來(lái)描述。

r■+■■■2+4■=■（12）

1 0 0 … 02 -3 1 …0 1 -2 1 … … 0 1 -2 1… … … … … 0 … 0 1 -2 10 … 0 1 -1p■p■p■p■？噎p■p■=a1p■/a■h′（t）h′（t）h′（t） ？噎h(huán)′（t）h（t）-a2■0

r12■1

r22■2

r32■3

？噎

r2N-1■N-1

r2N■N（13）

式（13）用向量表示可以簡(jiǎn)寫為：A

■■=a1■-a2■（14）

式（14）中：a■=■，■=0r12■1r12■1r22■2r32■3？噎r2N-1■N-1r2■■N■N。

在每個(gè)單元內(nèi)的每個(gè)氣泡的整個(gè)生命周期變化情況可以獲得（圖7）。

3 結(jié)論

在2個(gè)充液平行平面的動(dòng)態(tài)開啟過(guò)程中，氣穴氣泡先是生長(zhǎng)而后潰滅。在這個(gè)問(wèn)題研究的基礎(chǔ)上，介紹和對(duì)比研究了顯式歐拉法、隱式歐拉法、改進(jìn)歐拉法、龍格庫(kù)塔法和ODE45方法在數(shù)值分析ODE問(wèn)題上的應(yīng)用。結(jié)果表明，這些方法在解決ODE問(wèn)題上是互相吻合的，在步長(zhǎng)足夠小的情況下它們之間的誤差也非常小[14-16]。

這些方法在MATLAB環(huán)境下編程也不復(fù)雜，這樣就可以結(jié)合其他的程序獲得復(fù)雜問(wèn)題的數(shù)值解。看起來(lái)很抽象的矩陣其實(shí)它表示的是向量之間在時(shí)間、空間上的關(guān)系，一旦這些關(guān)系確立以后，通過(guò)分析就可以得到很直觀的解釋。另外，還介紹了解決PDE方程的方法。在建立兩平行平面的動(dòng)態(tài)開啟模型時(shí)，使用的是有限元法。

在實(shí)際的森林工程實(shí)踐中，有很多問(wèn)題需要用微分方程來(lái)刻畫[17-20]，本研究旨在說(shuō)明如何將高等數(shù)學(xué)知識(shí)與工程實(shí)踐更好地結(jié)合，在今后的森林工程課堂教學(xué)中加強(qiáng)這方面的運(yùn)用，讓學(xué)生在今后的工作中有能力做出精確的分析，從而取得更大的發(fā)展。

4 參考文獻(xiàn)

[1] 劉鴻基，劉剛.數(shù)學(xué)發(fā)展簡(jiǎn)史及未來(lái)趨勢(shì)簡(jiǎn)論[J].科技信息，2012（4）：129-130.

[2] 鄧明立，張生春.19世紀(jì)德國(guó)數(shù)學(xué)發(fā)展原因初探[J].自然辯證法研究，2000（8）：45-49.

[3] 賴發(fā)孝.淺析高中數(shù)學(xué)在社會(huì)經(jīng)濟(jì)發(fā)展中的作用[J].經(jīng)貿(mào)實(shí)踐，2016（18）：36.

[4] 張奠宙.二十世紀(jì)數(shù)學(xué)發(fā)展一瞥[J].自然雜志，1982（3）：179-182.

[5] 冷螢征.探討數(shù)學(xué)在經(jīng)濟(jì)發(fā)展中的地位及作用[C]//“決策論壇：企業(yè)黨建與政工創(chuàng)新工作發(fā)展學(xué)術(shù)研討會(huì)”論文集（下）.中國(guó)武漢決策信息研究開發(fā)中心、決策與信息雜志社、北京大學(xué)經(jīng)濟(jì)管理學(xué)院，2016：1.

[6] 曹建美.數(shù)學(xué)在現(xiàn)代經(jīng)濟(jì)學(xué)及經(jīng)濟(jì)發(fā)展中的作用[N].山西經(jīng)濟(jì)日?qǐng)?bào)，2009-12-14（004）.

[7] 張恭慶.數(shù)學(xué)與國(guó)家實(shí)力（上）[J].紫光閣，2014（8）：76-78.

[8] 王汝發(fā).從20世紀(jì)數(shù)學(xué)發(fā)展再看數(shù)學(xué)與科學(xué)技術(shù)之關(guān)系[J].廣東工業(yè)大學(xué)學(xué)報(bào)（社會(huì)科學(xué)版），2003（1）：57-60.

[9] RAYLEIGH L.On the Pressure Developed in a Liquid during the Colla-pse of a Spherical Cavity[J].Philosophical Magazine，1917，34：94-98.

[10] E HAIRER，S P NORSETT，G Wanner.Solving Ordinary Differential Equ-ations I-Nonstiff Problems[M].Second Revised Edition.Springer press，1992.

[11] SUN Zhizhong.Numerical solutions for partial differential equations[M].Science Publishing House，2005.

[12] DAVID BLEECKER，GEORGE CSORDAS.Basic Partial Differential Eq-uations[M].International Press of Boston，1997.

[13] ANDREI D，POLYANIN.HANDBOOK OF LINEAR PARTIAL DIFFER-ENTIAL EQUATIONS for ENGINEERS and SCIENTISTS （First Edition）[M].CHAPMAN & HALL/CRC，2002.

[14] 杜世拔.森林經(jīng)營(yíng)系統(tǒng)工程數(shù)學(xué)模型的研究[J].林業(yè)資源管理，1982（4）：27-29.

[15] 董斌，陳立平，錢國(guó)英.基于遙感的層次分析法和模糊數(shù)學(xué)模型綜合評(píng)價(jià)森林資源生態(tài)適宜性[J].自然資源學(xué)報(bào)，2011（3）：468-476.

[16] 吳發(fā)云，孫濤，王小昆.基于數(shù)學(xué)形態(tài)學(xué)的森林資源數(shù)據(jù)更新技術(shù)研究[J].林業(yè)資源管理，2009（5）：105-108.

[17] 成子純，呂勇.森林生物數(shù)學(xué)的發(fā)展與展望[J].中南林業(yè)調(diào)查規(guī)劃，1997（2）：47-49.

[18] 張映堂，郭安.云南楚雄州林火分布的數(shù)學(xué)模型及森林防火期的劃分[J].西南林學(xué)院學(xué)報(bào)，1996（2）：100-104.

[19] 劉春玲，路紫.數(shù)學(xué)方法在森林生態(tài)旅游區(qū)開發(fā)中的具體應(yīng)用[J].經(jīng)濟(jì)地理，2001（1）：118-120.

[20] 胥曉，蘇智先，黎云祥，等.嘉陵江流域南充金城山森林群落的模糊數(shù)學(xué)分析[J].四川師范學(xué)院學(xué)報(bào)（自然科學(xué)版），1999（2）：182-189.