一維輸運(yùn)方程的編程求解及可視化

黃海深，萬秋紅，李 平，周庭艷，吳 波

（遵義師范學(xué)院 物理與電子科學(xué)學(xué)院，貴州 遵義)

引言

數(shù)學(xué)物理方法是物理學(xué)本科班的專業(yè)必修課程，按照人才培養(yǎng)方案設(shè)計(jì)的理念，重點(diǎn)加強(qiáng)對學(xué)生自學(xué)能力的培養(yǎng)，學(xué)生應(yīng)當(dāng)掌握課程中的復(fù)變函數(shù)論和數(shù)學(xué)物理方程的基礎(chǔ)知識、基本研究方法和思想，并能夠在今后的學(xué)習(xí)和實(shí)踐中加以應(yīng)用。

數(shù)學(xué)物理方法是以研究物理問題為目標(biāo)的數(shù)學(xué)理論和數(shù)學(xué)方法，對物理現(xiàn)象建立數(shù)學(xué)模型，尋求物理現(xiàn)象的數(shù)學(xué)描述和詮釋。數(shù)學(xué)物理方法主要包括復(fù)變函數(shù)論、特殊函數(shù)論和數(shù)學(xué)物理方程等內(nèi)容，涉及到許多偏微分方程和特殊方程。數(shù)學(xué)物理方法以普通物理和高等數(shù)學(xué)為基礎(chǔ)，內(nèi)容具有“繁、雜、難”的特點(diǎn)，加上學(xué)生基礎(chǔ)薄弱，學(xué)生學(xué)習(xí)起來比較困難，成為學(xué)生眼中的“天書”，造成學(xué)習(xí)主動性不高，學(xué)習(xí)效果較差[1]。

計(jì)算科學(xué)的發(fā)展和計(jì)算機(jī)處理能力的大幅提升，為數(shù)學(xué)物理中的許多問題提供了另一個(gè)解決思路，并發(fā)展出了一個(gè)新的學(xué)科，即計(jì)算物理學(xué)。計(jì)算物理學(xué)是以計(jì)算機(jī)技術(shù)為工具和手段，運(yùn)用數(shù)學(xué)方法解決復(fù)雜物理問題的一門應(yīng)用科學(xué)，它是計(jì)算機(jī)科學(xué)、數(shù)學(xué)和物理學(xué)相結(jié)合的學(xué)科。目前計(jì)算物理已經(jīng)成為對復(fù)雜體系的物理規(guī)律、物理性質(zhì)進(jìn)行研究的重要手段，是物理學(xué)的重要分支之一，與理論物理和實(shí)驗(yàn)物理有著密切的關(guān)系，對物理學(xué)的發(fā)展起著越來越大的作用。

計(jì)算物理首先對復(fù)雜物理問題建立物理模型，然后用數(shù)學(xué)理論導(dǎo)出方程，最后利用數(shù)值計(jì)算方法來求解和仿真，可以很直觀地對問題進(jìn)行分析，因此計(jì)算物理在數(shù)學(xué)物理方法課程的許多章節(jié)中發(fā)揮作用。本文以一維輸運(yùn)方程為例，研究有限差分方法在偏微分方程的數(shù)值求解和可視化方面的應(yīng)用，為建立數(shù)學(xué)物理方法課程的數(shù)值求解及可視化資源庫打下基礎(chǔ)。

1 輸運(yùn)方程的建立

當(dāng)物體的溫度不均勻時(shí)，熱量就會從溫度高的地方向溫度低的地方轉(zhuǎn)移，這種熱傳導(dǎo)現(xiàn)象遵從熱傳導(dǎo)定律，即傅里葉定律：

式中：u 為溫度，它是位置r 和時(shí)間t 的函數(shù)；?u為溫度的梯度，表示u 的不均勻程度；q 為單位時(shí)間內(nèi)流過單位橫截面積的熱量，即熱流密度，比例系數(shù)k 為熱傳導(dǎo)系數(shù)。

根據(jù)熱傳導(dǎo)定律及能量守恒定律，可推導(dǎo)出物體的溫度滿足以下偏微分方程[2]：

當(dāng)物質(zhì)的濃度不均勻時(shí)，物質(zhì)的擴(kuò)散現(xiàn)象也遵循與傅里葉定律相似的擴(kuò)散定律，即斐克定律，因此熱傳導(dǎo)方程和物質(zhì)的擴(kuò)散方程具有相同的形式：

此類方程即輸運(yùn)方程。為簡單起見，本文只考慮比較簡單的一維的情況，即：

2 輸運(yùn)方程的求解

對于齊次輸運(yùn)方程，如果邊界條件也是齊次的，一般可用分離變量方法進(jìn)行求解。例如對于定解問題：

式(5)化為

關(guān)于x 部分的解已經(jīng)求出，將本征值代入式(11)可得：

此為關(guān)于t 部分的解。考慮初始條件式(8)后，可求出定解問題的解為[3]：

此題的初始條件是三角函數(shù)，而且正好為式(11)的本征函數(shù)，所以最后求解特解時(shí)相對容易些，如果不是本征函數(shù)，則需要采用傅里葉變換求解展開系數(shù)，求解過程則更為復(fù)雜。對于一般的輸運(yùn)方程，求解過程是非常復(fù)雜的，甚至求不出解析解。因此，接下來我們采用有限差分方法求輸運(yùn)方程的數(shù)值解并對數(shù)值進(jìn)行可視化。

3 數(shù)值求解及可視化

有限差分方法是求解偏微分方程的一種常用方法，它使用差商代替偏微分方程中的偏導(dǎo)數(shù)，這樣偏微分方程就變成了代數(shù)方程，連續(xù)方程實(shí)現(xiàn)了離散化，再聯(lián)合邊界條件和初值條件即可求解出方程的數(shù)值解。

下面是采用有限差分方法對一維輸運(yùn)方程進(jìn)行數(shù)值求解并進(jìn)行可視化。解為u(x,t)=ex+t輸運(yùn)方程為：

選用此方程的目的是其解析解已知，方便評估我們的求解方法和過程的可靠性。

則式(20-23)可離散為[4]

至此通過C++或其他語言編程可很方便地求解出數(shù)值解，將數(shù)值導(dǎo)入到作圖軟件中即可畫出函數(shù)圖像，對輸運(yùn)方程的解進(jìn)行可視化。例如取M=50、N=5000 時(shí)，圖像如圖1 所示。由圖1 可以直觀地看出，函數(shù)值隨時(shí)間是增加的，其原因是邊界處函數(shù)值是隨時(shí)間增加的，相當(dāng)于兩端有兩個(gè)熱源不斷地對體系進(jìn)行加熱，所以溫度一直在上升。

圖1 M=50、N=5000 時(shí)式(20-23)的圖像

有限差分方法是使用有限大小的間隔對求解區(qū)域時(shí)行劃分，并用差商代替偏導(dǎo)數(shù)，因此迭代的次數(shù)越多，誤差越大，即時(shí)間越長，誤差越大，如圖2 所示。在本例中，t=1 時(shí)最大的誤差在x=0.5 附近，約為3.48×10-5，邊界處的值已知，誤差為0。除了迭代次數(shù)，誤差還與網(wǎng)格大小相關(guān)，網(wǎng)格劃分得越細(xì)，誤差越小。但也不是越細(xì)越好。網(wǎng)格越細(xì)，計(jì)算的數(shù)據(jù)越多，處理起來越困難。而且許多作圖軟件能處理的精度有限，過高的精度并不能體現(xiàn)出來。

圖2 M=50、N=5000、t=1 時(shí)，采用有限差分方法求解式(20-23)的誤差圖像

對于式(5-8)的定解問題，取D=1，M=90，N=5000，圖像如圖3 所示。由圖3 看出，函數(shù)值隨時(shí)間衰減得很快，t=1 時(shí)，最大值已經(jīng)衰減到10-11的數(shù)量級。由邊界條件可知邊界的值固定為0，對于熱傳導(dǎo)問題相當(dāng)于兩端有溫度為0 的熱浴，經(jīng)過一段很短的時(shí)間，體系的溫度將接近熱浴的溫度，即趨于0。

圖3 M=90、N=5000、D=1 時(shí)，式(5-8)的圖像

再看一個(gè)齊次方程的例子：

初始條件為0，邊界條件為關(guān)于t 的正弦函數(shù)，兩端的符號相反。此方程的邊界條件是非齊次的，求解非常麻煩，采用有限差分方法求出數(shù)值解，由origin 軟件畫出的函數(shù)圖像如圖4 所示。由圖4 可以看出，兩端都以正弦函數(shù)變化，向中間輸運(yùn)。因兩端符號相反，造成x＜0.5 時(shí)函數(shù)始終大于0，x＞0.5 時(shí)函數(shù)始終小于0，而x=0.5 時(shí)函數(shù)始終為0。從圖中可以看出，輸運(yùn)方程與波動方程的圖像有著本質(zhì)區(qū)別，后者的圖像遵循疊加原理，而前者不遵循疊加原理，而會相互抵消。

圖4 M=50、N=5000 時(shí)，式(31-34)的圖像

4 結(jié)論

本文首先簡單介紹了一維輸運(yùn)方程的建立和基于分離變量法的解析解求解方法，然后采用有限差分方法對一維輸運(yùn)方程進(jìn)行了離散化處理，最后經(jīng)過編程求出了兩個(gè)定解問題的數(shù)值解，并進(jìn)行了可視化。該方法一旦編程成功，改變邊界條件和初始條件也可很方便地求解出方程的數(shù)值解并進(jìn)行畫圖。該方法可移植性較高，直觀性較強(qiáng)，可降低學(xué)生學(xué)習(xí)輸運(yùn)方程的難度，提高其學(xué)習(xí)興趣。