輸運(yùn)方程數(shù)值解法分析

王海林，魏丙濤，劉天祥，桂朝覲，高全歸

（1.玉溪師范學(xué)院 理學(xué)院，云南 玉溪 653100；2.文山學(xué)院 數(shù)理系，云南 文山 663000）

輸運(yùn)方程是數(shù)學(xué)物理方程中的一類最基本的方程，在許多情況下，由于邊界條件比較復(fù)雜，或者由于方程比較復(fù)雜，或由于其他各種原因，難以嚴(yán)格解出［1］。因此，近似解法得以發(fā)展，只要近似程度能夠滿足實(shí)際需要，近似解的價(jià)值并不低于嚴(yán)格解，隨著電子計(jì)算機(jī)的日漸普及，用數(shù)值方法求解方程的近似解變得非常方便。

數(shù)值解法求解偏微分方程的方法有多種，應(yīng)用比較普遍的有有限差分法，變分法和有限元方法，有限差分法一般應(yīng)用于求解含時(shí)間問題比較方便。解同樣的離散微分方程，采用好的算法與采用一般的算法的計(jì)算效果往往相差很大，采用好的算法不但能使求解過程數(shù)值穩(wěn)定、數(shù)值解的精度得到提高，而且能數(shù)十倍、數(shù)百倍地節(jié)省計(jì)算工作量［2］。本文將通過不同的數(shù)值算法求解相同的定解問題，給出數(shù)值方法求解輸運(yùn)方程的一些有趣的結(jié)果。

1 定解問題

在數(shù)學(xué)物理方程中，輸運(yùn)方程有很多，比如擴(kuò)散方程、熱傳導(dǎo)方程等。在此，用一個(gè)比較簡(jiǎn)單的方程，來考察差分格式對(duì)數(shù)值結(jié)果的影響，以下面的定解問題為例：

當(dāng)然上述定解問題可以用分離變量的方法求得其解析解為：

于是可以給出解析解的函數(shù)圖像，這里給出函數(shù)值u隨x和t變化的函數(shù)關(guān)系圖（見圖1）：

從圖1中可以看出，隨著時(shí)間的增加，函數(shù)值逐漸衰減，函數(shù)值隨空間坐標(biāo)按正弦規(guī)律變化，這與定解問題給出的解析解是一致的，需要說明的是在得到解析解的過程中已經(jīng)考慮了初始條件和邊界條件。

圖1 函數(shù)值u隨x和t變化的解析解的函數(shù)圖像

2 數(shù)值方法

2.1 方程及初邊界條件的離散化

對(duì)上述方程的離散化方法有很多，在此對(duì)方程采用兩種離散方法：向前差分和向后差分，離散后向前差分格式的方程可以寫為［3］：

這里使用U（x，tn）僅僅是為了和微分方程中的函數(shù)值u相區(qū)分，于是（4）式可以改寫為

對(duì)（6）式等式左右兩邊同時(shí)做傅里葉積分可得：

上式可以化簡(jiǎn)為：

離散后的向后差分格式的方程可以寫為［5］：

這個(gè)差分格式是一個(gè)絕對(duì)穩(wěn)定的格式，雖然這樣的差分格式的截?cái)嗾`差也為，但是，這樣的離散在實(shí)際的計(jì)算中會(huì)存在比較大的誤差。

由于定解問題中的邊界全為0，這樣邊界條件很好處理，在邊界處直接取0，但初始條件需要離散化，初始條件取為：

有了上面的離散化后就可以進(jìn)行數(shù)值求解了。數(shù)值求解的過程中使用Visual C++6.0作為編譯器，圖形用Origin8.5畫出。

2.2 不同差分格式對(duì)結(jié)果的影響

圖2和圖3分別運(yùn)用向前差分的方法計(jì)算了t=0.5時(shí)的定解問題的解，圖2的求解中取a2λ=0.5，剛好能滿足穩(wěn)定性條件，圖3不滿足穩(wěn)定性條件，在圖3中取a2λ=1。可見，計(jì)算過程中不滿足穩(wěn)定性條件a2λ≤0.5時(shí)，數(shù)值解法是不穩(wěn)定的，不穩(wěn)定的解沒有意義。表1中給出了向前差分格式對(duì)應(yīng)的不同網(wǎng)格比得到的數(shù)值解。由表1中可以看出，當(dāng)a2λ=1時(shí)，函數(shù)值得振蕩特別嚴(yán)重，這是不穩(wěn)定的表現(xiàn)。

圖2 t=0.5時(shí)網(wǎng)格比為1/2時(shí)的函數(shù)圖

圖3 t=0.5時(shí)網(wǎng)格比大于1/2時(shí)的函數(shù)圖

圖4 t=0.1時(shí)解析解和兩種差分結(jié)果的比較

圖5 t=0.5時(shí)解析解和兩種差分結(jié)果的比較

圖4給出了t=0.1時(shí)解析解和兩種數(shù)值解法的比較，由圖4可以看出，三種解法的誤差很小，而向前差分和解析解的結(jié)果更加接近。圖五給出了t=0.5時(shí)解析解與兩種差分格式的數(shù)值解的比較，其中點(diǎn)線為解析解的結(jié)果，實(shí)線為向前差分格式的數(shù)值結(jié)果，虛線為向后差分的數(shù)值結(jié)果。有趣的是向后差分格式雖然在任何情況下都是穩(wěn)定的，但是，誤差比較大，隨著數(shù)值計(jì)算中不斷的迭代，誤差越來越大，而向前差分只要滿足穩(wěn)定性條件，誤差比向后差分要小的多。這是容易理解的，因?yàn)樵跀?shù)值計(jì)算中，向前差分是由前一個(gè)時(shí)間層上的三個(gè)點(diǎn)來計(jì)算下一個(gè)時(shí)間層上一個(gè)點(diǎn)的值，而向后差分是在兩個(gè)時(shí)間層上建立線性方程組來求解，因而，其誤差更大。表2中分別給出了解析解、向前差分格式和向后差分格式在t=0.1時(shí)及t=0.5時(shí)的數(shù)值結(jié)果，從表2中的數(shù)值解的結(jié)果可以看出當(dāng)t=0.1時(shí)數(shù)值解的誤差較小，但當(dāng)t=0.5時(shí)向后差分的誤差就變大了。需要指出的是解析解用數(shù)值形式表出的時(shí)候依然是存在誤差的，主要是由于π的取值是存在誤差，但這個(gè)誤差很小，在很多實(shí)際應(yīng)用中時(shí)可以忽略的。

表1 向前差分格式對(duì)應(yīng)不同的網(wǎng)格比得到的數(shù)值解

（續(xù)表1）

表2 解析解、向前差分格式和向后差分格式在t=0.1時(shí)及t=0.5時(shí)的數(shù)值

（續(xù)表2）

（續(xù)表2）

3 結(jié)論與展望

通過上面的討論，可以看出，輸運(yùn)方程的數(shù)值解法需要考慮數(shù)值解的穩(wěn)定性，對(duì)于本文中給出的兩種差分方法，向前差分是條件穩(wěn)定的，其穩(wěn)定條件為a2λ≤0.5，而向后差分是絕對(duì)穩(wěn)定的，在實(shí)際的數(shù)值計(jì)算中，兩種差分格式如果對(duì)于空間變量x取相同的步長(zhǎng)，向后差分對(duì)時(shí)間步長(zhǎng)沒有嚴(yán)格的限制，但是，向前差分就必須對(duì)時(shí)間步長(zhǎng)加以限制，這樣導(dǎo)致時(shí)間步長(zhǎng)需要取的很小，就會(huì)大大提升計(jì)算量，而向后差分可以將時(shí)間步長(zhǎng)取的相對(duì)較長(zhǎng)，但這樣的差分格式的計(jì)算精度不是很高，于是在很多需要高精度高效率的工程技術(shù)問題中，上面的方法實(shí)用性就變得比較有限了，所以我們需要尋找一些高效率高精度的計(jì)算方法來滿足實(shí)際需求，比如加權(quán)隱格式，多層網(wǎng)格方法等，這將是我們接下來需要開展的工作。

［1］梁昆淼.數(shù)學(xué)物理方法［M］.北京：高等教育出版社，2002：459.

［2］徐長(zhǎng)發(fā)，李紅.微分方程數(shù)值解法［M］.武漢：華中科技大學(xué)出版社，2002：381.

［3］陸金甫，關(guān)治.偏微分方程數(shù)值解法［M］.北京：清華大學(xué)出版社，2007：82.

［4］Mitchell A R, Griffiths D F.The Finite Difference Methed in Partial Dofferential Equations［M］.New York： John Wiley & Sons, 1967:259.

［5］張文生.科學(xué)計(jì)算中的偏微分方程有限差分格式［M］.北京：高等教育出版社，2006：236.