基于Python求解偏微分方程的有限差分法

王登岳　張宏偉

（西安石油大學(xué)電子工程學(xué)院光電油氣測井與檢測教育部重點(diǎn)實(shí)驗(yàn)室， 陜西 西安 710065）

摘 要： 偏微分方程的求解是很多科學(xué)技術(shù)問題的關(guān)鍵難點(diǎn)。隨著計(jì)算機(jī)性能的不斷提高，數(shù)值解法能夠解復(fù)雜的偏微分方程并將計(jì)算結(jié)果圖形化。相對(duì)于昂貴的科學(xué)計(jì)算軟件，Python是一種免費(fèi)的面向?qū)ο?、?dòng)態(tài)的程序設(shè)計(jì)語言。有限差分法以其概念清晰，方法簡單、直觀等特點(diǎn)在偏微分方程的求解中得到了廣泛的應(yīng)用。文章對(duì)矩形區(qū)域的拉普拉斯方程進(jìn)行數(shù)值求解，采用Numpy對(duì)有限差分法進(jìn)行計(jì)算，運(yùn)用Matplotlib繪制等值線，輸出迭代次數(shù)以及誤差。

關(guān)鍵詞： 偏微分方程； Python； 數(shù)值解法； 有限差分法

中圖分類號(hào)：TP3 文獻(xiàn)標(biāo)志碼：A 文章編號(hào)：1006-8228（2016）11-14-03

Python programmed finite difference method for solving partial differential equations

Wang Dengyue， Zhang Hongwei

（Xi'an Shiyou University， Electric Engineering College， Key Laboratory of Photo Electricity Gas and Oil Detecting of Ministry of Education， Xi'an， Shaanxi 710065， China）

Abstract： To solve the partial differential equations （PDE） is a key difficult point in many scientific and technical problems. With the development of computer performance， numeric solution can solve many sophisticated PDE and visualize the numeric results. Rather than the expensive science computing software， Python is a free object-oriented language， dynamic programming language. Finite difference method （FDM） is widely used for its clear， simple and intuitive. Laplace problem in a rectangular area is solved numerically in the article， computed through FDM with the Numpy library， visualized through plotting the contour by the Matplotlib library， and the number of iteration and the error are given.

Key words： partial difference equation； Python； numeric solution； finite difference method

0 引言

在數(shù)學(xué)中，偏微分方程是包含多變量和它們的偏導(dǎo)數(shù)在內(nèi)的微分方程。偏微分方程通常被用來求解聲、熱、靜態(tài)電場、動(dòng)態(tài)電場、流體、彈性力學(xué)或者量子力學(xué)方面的問題[1]。這些現(xiàn)象能夠被模式化的偏微分方程描述，正如一維動(dòng)態(tài)系統(tǒng)通常會(huì)用常微分方程描述。為了更深入地理解上述各種現(xiàn)象，求解偏微分方程成為理解以及解釋上述現(xiàn)象的關(guān)鍵。

1 Python及相關(guān)模塊簡介

Python是一種面向?qū)ο?、?dòng)態(tài)的程序設(shè)計(jì)語言。具有非常簡潔而清晰的語法，適合完成各種高層任務(wù)。它既可用來快速開發(fā)程序腳本，也可用來開發(fā)大規(guī)模的軟件。隨著Numpy，SciPy，Matplotlib等眾多程序庫的開發(fā)，Python同樣適合于做科學(xué)計(jì)算以及繪制高質(zhì)量的2D和3D圖像。與科學(xué)計(jì)算領(lǐng)域的商業(yè)軟件Matlab相比，Python是一門通用的程序設(shè)計(jì)語言，比Matlab所采用的腳本語言的應(yīng)用更廣泛，有更多的程序庫的支持。

Numpy是使用Python進(jìn)行科學(xué)計(jì)算的基礎(chǔ)包：①它的基本類型是N-維陣列對(duì)象；②提供了功能強(qiáng)大的函數(shù)；③可以提供C/C++與Fortran代碼的接口；④強(qiáng)大的線性代數(shù)計(jì)算，傅里葉變換以及隨機(jī)數(shù)計(jì)算[2]。SciPy由一系列的數(shù)值計(jì)算和特定領(lǐng)域的工作箱構(gòu)成，常用的工具箱包括信號(hào)處理，優(yōu)化算法以及統(tǒng)計(jì)等。

Matplotlib能夠產(chǎn)生不同格式高質(zhì)量的圖片，該軟件可以在不同的平臺(tái)上使用，例如采用Python腳本，提供python或者Ipython的命令行接口（類似于Matlab或Mathematica）等等[3]。Matplotlib能夠繪制直方圖、功率譜圖、殘差圖以及散點(diǎn)圖等等。本次我們采用Matplotlib中的簡單的繪圖對(duì)象Pyplot，Pyplot提供了一個(gè)類Matlab的接口，可以通過Pyplot完全控制線型、字體特征和坐標(biāo)軸特征等。

2 PDE解法簡介

求解偏微分方程可以通過解析方法求解，解析法包括建立和求解偏微分方程，嚴(yán)格求解偏微分方程的經(jīng)典方法是分離變量法。解析法的優(yōu)點(diǎn)是可以求解出表示為已知函數(shù)的顯式，從而得到精確的結(jié)果。但解析法存在嚴(yán)重的缺點(diǎn)，只有在為數(shù)不多的坐標(biāo)系中進(jìn)行分離變量，對(duì)更復(fù)雜的模型很難通過解析計(jì)算得到。隨著計(jì)算機(jī)性能的日益提升，數(shù)值計(jì)算在近年來取得來巨大的進(jìn)步，不僅能夠?qū)ΤＲ?guī)的偏微分進(jìn)行數(shù)值驗(yàn)證，而且也能夠?qū)δＰ蛷?fù)雜的偏微分方程進(jìn)行計(jì)算。常用的數(shù)值方法主要包括有限差分法，有限元法以及有限體積法。數(shù)值方法適合求解模型比較復(fù)雜的偏微分方程[4]。

3 Python解PDE流程

在數(shù)值計(jì)算方法中，有限差分法是應(yīng)用最早的一種方法，直到今天，它仍然以簡單、直觀的特點(diǎn)而應(yīng)用廣泛[5]。不論是常微分方程還是偏微分方程，對(duì)初值問題或者邊值問題、橢圓型、雙曲型或拋物型二階線性方程，以及高階或非線性方程，通常都可以采用有限差分法將它們轉(zhuǎn)化為代數(shù)方程組，再借助計(jì)算機(jī)求其數(shù)值解。

本文結(jié)合一個(gè)簡單的電場的例子，給出求解偏微分方程的數(shù)值解法的過程。一個(gè)長直接地金屬槽中的電場，其側(cè)壁與底面電位均為0，頂蓋電位的相對(duì)值為100。對(duì)于此槽中間區(qū)段的電場分析，可以理想化為二維平行平面場問題。選定直角坐標(biāo)系如圖1所示，槽內(nèi)電位函數(shù)滿足拉普拉斯方程，構(gòu)成如下的第一類邊值問題：

按照有限差分法的計(jì)算步驟，上述給出的例子的解題過程如下。

⑴ 場域離散化：用正方形網(wǎng)格對(duì)圖示區(qū)域進(jìn)行剖分，步距為1，剖分為100個(gè)點(diǎn)。

⑵ 給出高斯-塞德爾迭代法的差分方程：

其中程序中采用的不是單點(diǎn)迭代計(jì)算，而是將所有的點(diǎn)分成如圖2所示的矩陣塊，圖中陰影部分矩陣代表N-2維矩陣，進(jìn)行迭代計(jì)算。

給定邊界條件以及內(nèi)點(diǎn)的初值，側(cè)壁與底面電位均為0，頂蓋電位的相對(duì)值為100。

給定迭代解收斂判定指標(biāo)：當(dāng)每一點(diǎn)相鄰兩次迭代值得絕對(duì)誤差小于可以接收的范圍時(shí)；迭代次數(shù)過多大于所能接受的范圍時(shí)。當(dāng)上述兩者任意一個(gè)條件滿足，則停止迭代，輸出計(jì)算結(jié)果，迭代次數(shù)n=5991，誤差m=0.001，表明經(jīng)過5991次后誤差已經(jīng)達(dá)到我們所需的精度。

繪制電場分布的等值線圖，所得的結(jié)果如圖3所示。

上述計(jì)算的整個(gè)過程如圖4所示。

4 結(jié)束語

偏微分方程的求解在科學(xué)和工程計(jì)算中占據(jù)重要地位。雖然Matlab能夠進(jìn)行大部分PDE方程求解，但是價(jià)格昂貴，對(duì)三維的PDE方程組的支持并不完善，而且因其閉源軟件屬性，對(duì)軟件的擴(kuò)展使用必須經(jīng)過廠商認(rèn)證之后，廠商才會(huì)添加其接口。Python進(jìn)行科學(xué)計(jì)算時(shí)只需要根據(jù)功能選擇好相關(guān)的庫即可，Numpy模塊以及Matplotlib模塊中的pyplot子模塊便能夠解決矩形區(qū)域的拉普拉斯方程求解以及將其圖形化。現(xiàn)實(shí)中的求解問題通常并不是針對(duì)規(guī)則區(qū)域，因此對(duì)于不同類型的模型，建模問題還有待進(jìn)一步研究，并將其他軟件中的模型導(dǎo)入作為更深層次研究的重點(diǎn)內(nèi)容。

參考文獻(xiàn)（References）：

[1] Partial differential equation. [EB/OL]. [2016-08-11].

https：//en.wikipedia.org/wiki/Partial_differential_equation

[2] Numpy.[EB/OL] [2016-08-11].http：//www.numpy.org/

[3] Matplotlib.[EB/OL] [2016-08-11].http：//matplotlib.org/

index.html

[4] 王秉中.計(jì)算電磁學(xué)[M].科學(xué)出版社，2002.

[5] 倪光正，楊仕友，邱捷.工程電磁場數(shù)值計(jì)算[M].機(jī)械工業(yè)出

版社，2010.