黃海深,萬秋紅,李 平,周庭艷,吳 波
(遵義師范學(xué)院 物理與電子科學(xué)學(xué)院,貴州 遵義)
數(shù)學(xué)物理方法是物理學(xué)本科班的專業(yè)必修課程,按照人才培養(yǎng)方案設(shè)計(jì)的理念,重點(diǎn)加強(qiáng)對(duì)學(xué)生自學(xué)能力的培養(yǎng),學(xué)生應(yīng)當(dāng)掌握課程中的復(fù)變函數(shù)論和數(shù)學(xué)物理方程的基礎(chǔ)知識(shí)、基本研究方法和思想,并能夠在今后的學(xué)習(xí)和實(shí)踐中加以應(yīng)用。
數(shù)學(xué)物理方法是以研究物理問題為目標(biāo)的數(shù)學(xué)理論和數(shù)學(xué)方法,對(duì)物理現(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í)主動(dòng)性不高,學(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)成為對(duì)復(fù)雜體系的物理規(guī)律、物理性質(zhì)進(jìn)行研究的重要手段,是物理學(xué)的重要分支之一,與理論物理和實(shí)驗(yàn)物理有著密切的關(guān)系,對(duì)物理學(xué)的發(fā)展起著越來越大的作用。
計(jì)算物理首先對(duì)復(fù)雜物理問題建立物理模型,然后用數(shù)學(xué)理論導(dǎo)出方程,最后利用數(shù)值計(jì)算方法來求解和仿真,可以很直觀地對(duì)問題進(jìn)行分析,因此計(jì)算物理在數(shù)學(xué)物理方法課程的許多章節(jié)中發(fā)揮作用。本文以一維輸運(yùn)方程為例,研究有限差分方法在偏微分方程的數(shù)值求解和可視化方面的應(yīng)用,為建立數(shù)學(xué)物理方法課程的數(shù)值求解及可視化資源庫打下基礎(chǔ)。
當(dāng)物體的溫度不均勻時(shí),熱量就會(huì)從溫度高的地方向溫度低的地方轉(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)方程。為簡(jiǎn)單起見,本文只考慮比較簡(jiǎn)單的一維的情況,即:
對(duì)于齊次輸運(yùn)方程,如果邊界條件也是齊次的,一般可用分離變量方法進(jìn)行求解。例如對(duì)于定解問題:
式(5)化為
關(guān)于x 部分的解已經(jīng)求出,將本征值代入式(11)可得:
此為關(guān)于t 部分的解??紤]初始條件式(8)后,可求出定解問題的解為[3]:
此題的初始條件是三角函數(shù),而且正好為式(11)的本征函數(shù),所以最后求解特解時(shí)相對(duì)容易些,如果不是本征函數(shù),則需要采用傅里葉變換求解展開系數(shù),求解過程則更為復(fù)雜。對(duì)于一般的輸運(yùn)方程,求解過程是非常復(fù)雜的,甚至求不出解析解。因此,接下來我們采用有限差分方法求輸運(yùn)方程的數(shù)值解并對(duì)數(shù)值進(jìn)行可視化。
有限差分方法是求解偏微分方程的一種常用方法,它使用差商代替偏微分方程中的偏導(dǎo)數(shù),這樣偏微分方程就變成了代數(shù)方程,連續(xù)方程實(shí)現(xiàn)了離散化,再聯(lián)合邊界條件和初值條件即可求解出方程的數(shù)值解。
下面是采用有限差分方法對(duì)一維輸運(yùn)方程進(jìn)行數(shù)值求解并進(jìn)行可視化。解為u(x,t)=ex+t輸運(yùn)方程為:
選用此方程的目的是其解析解已知,方便評(píng)估我們的求解方法和過程的可靠性。
則式(20-23)可離散為[4]
至此通過C++或其他語言編程可很方便地求解出數(shù)值解,將數(shù)值導(dǎo)入到作圖軟件中即可畫出函數(shù)圖像,對(duì)輸運(yùn)方程的解進(jìn)行可視化。例如取M=50、N=5000 時(shí),圖像如圖1 所示。由圖1 可以直觀地看出,函數(shù)值隨時(shí)間是增加的,其原因是邊界處函數(shù)值是隨時(shí)間增加的,相當(dāng)于兩端有兩個(gè)熱源不斷地對(duì)體系進(jìn)行加熱,所以溫度一直在上升。
圖1 M=50、N=5000 時(shí)式(20-23)的圖像
有限差分方法是使用有限大小的間隔對(duì)求解區(qū)域時(shí)行劃分,并用差商代替偏導(dǎo)數(shù),因此迭代的次數(shù)越多,誤差越大,即時(shí)間越長(zhǎng),誤差越大,如圖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)的誤差圖像
對(duì)于式(5-8)的定解問題,取D=1,M=90,N=5000,圖像如圖3 所示。由圖3 看出,函數(shù)值隨時(shí)間衰減得很快,t=1 時(shí),最大值已經(jīng)衰減到10-11的數(shù)量級(jí)。由邊界條件可知邊界的值固定為0,對(duì)于熱傳導(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ù),兩端的符號(hào)相反。此方程的邊界條件是非齊次的,求解非常麻煩,采用有限差分方法求出數(shù)值解,由origin 軟件畫出的函數(shù)圖像如圖4 所示。由圖4 可以看出,兩端都以正弦函數(shù)變化,向中間輸運(yùn)。因兩端符號(hào)相反,造成x<0.5 時(shí)函數(shù)始終大于0,x>0.5 時(shí)函數(shù)始終小于0,而x=0.5 時(shí)函數(shù)始終為0。從圖中可以看出,輸運(yùn)方程與波動(dòng)方程的圖像有著本質(zhì)區(qū)別,后者的圖像遵循疊加原理,而前者不遵循疊加原理,而會(huì)相互抵消。
圖4 M=50、N=5000 時(shí),式(31-34)的圖像
本文首先簡(jiǎn)單介紹了一維輸運(yùn)方程的建立和基于分離變量法的解析解求解方法,然后采用有限差分方法對(duì)一維輸運(yùn)方程進(jìn)行了離散化處理,最后經(jīng)過編程求出了兩個(gè)定解問題的數(shù)值解,并進(jìn)行了可視化。該方法一旦編程成功,改變邊界條件和初始條件也可很方便地求解出方程的數(shù)值解并進(jìn)行畫圖。該方法可移植性較高,直觀性較強(qiáng),可降低學(xué)生學(xué)習(xí)輸運(yùn)方程的難度,提高其學(xué)習(xí)興趣。