有限差分法求解偏微分方程

1、有限差分法求解偏微分方程摘要：本文主要使用有限差分法求解計(jì)算力學(xué)中的系統(tǒng)數(shù)學(xué)模型，推導(dǎo)了有限差分法的理論基礎(chǔ)，并在此基礎(chǔ)上給出了部分有限差分法求解偏微分方程的算例驗(yàn)證了推導(dǎo)的正確性及操作可行性。關(guān)鍵詞：計(jì)算力學(xué)，偏微分方程，有限差分法Abstract：This dissertation mainly focuses on solving the mathematic model of computation mechanics with finite-difference method. The theoretical basis of finite-difference is derived

2、 in the second part of the dissertation, and then I use MATLAB to program the algorithms to solve some partial differential equations to confirm the correctness of the derivation and the feasibility of the method.Key words：Computation Mechanics, Partial Differential Equations, Finite-Difference Meth

3、od1 引言機(jī)械系統(tǒng)設(shè)計(jì)常常需要從力學(xué)觀點(diǎn)進(jìn)行結(jié)構(gòu)設(shè)計(jì)以及結(jié)構(gòu)分析，而這些分析的前提就是建立工程問題的數(shù)學(xué)模型。通過對(duì)機(jī)械系統(tǒng)應(yīng)用自然的基本定律和原理得到帶有相關(guān)邊界條件和初始條件的微分積分方程，這些微分積分方程構(gòu)成了系統(tǒng)的數(shù)學(xué)模型。求解這些數(shù)學(xué)模型的方法大致分為解析法和數(shù)值法兩種，而解析法的局限性眾所周知，當(dāng)系統(tǒng)的邊界條件和受載情況復(fù)雜一點(diǎn)，往往求不出問題的解析解或近似解。另一方面，計(jì)算機(jī)技術(shù)的發(fā)展使得計(jì)算更精確、更迅速。因此，對(duì)于絕大多數(shù)工程問題，研究其數(shù)值解法更具有實(shí)用價(jià)值。對(duì)于微分方程而言，主要分為差分法和積分法兩種，本論文主要討論差分法。2 有限差分法理論基礎(chǔ)2.1 有限差分法的基本

4、思想當(dāng)系統(tǒng)的數(shù)學(xué)模型建立后，我們面對(duì)的主要問題就是微分積分方程的求解?；舅枷胧怯秒x散的只含有限個(gè)未知量的差分方程組去近似地代替連續(xù)變量的微分方程和定解條件，并把差分方程組的解作為微分方程定解問題的近似解。將原方程及邊界條件中的微分用差分來近似，對(duì)于方程中的積分用求和或及機(jī)械求積公式來近似代替，從而把原微分積分方程和邊界條件轉(zhuǎn)化成差分方程組。有限差分法求解偏微分方程的步驟主要有以下幾步：n 區(qū)域離散，即把所給偏微分方程的求解區(qū)域細(xì)分成由有限個(gè)格點(diǎn)組成的網(wǎng)格，這些離散點(diǎn)稱作網(wǎng)格的節(jié)點(diǎn)；n 近似替代，即采用有限差分公式替代每一個(gè)格點(diǎn)的導(dǎo)數(shù)；n 逼近求解，換而言之，這一過程可以看作是用一個(gè)插值多項(xiàng)式

5、及其微分來代替偏微分方程的解的過程。從原則上說，這種方法仍然可以達(dá)到任意滿意的計(jì)算精度。因?yàn)榉匠痰倪B續(xù)數(shù)值解可以通過減小獨(dú)立變量離散取值的間格，或者通過離散點(diǎn)上的函數(shù)值進(jìn)行插值計(jì)算來近似得到。理論上，當(dāng)網(wǎng)格步長(zhǎng)趨近于零時(shí)，差分方程組的解應(yīng)該收斂于精確解，但由于機(jī)器字節(jié)的限制，網(wǎng)格步長(zhǎng)不可能也沒有必要取得無限小，那么差分法的收斂性或者說算法的穩(wěn)定性就顯得至關(guān)重要。因此，在運(yùn)用有限差分法時(shí)，除了要保證精度外，還必須要保證其收斂性。2.2 系統(tǒng)微分方程的一般形式（1）由于大多數(shù)工程問題都是二維問題，所以得到的微分方程一般都是偏微分方程，對(duì)于一維問題得到的是常微分方程，解法與偏微分方程類似，故為了不是

6、一般性，這里只討論偏微分方程。由于工程中高階偏微分較少出現(xiàn)，所以本文僅僅給出二階偏微分方程的一般形式，對(duì)于高階的偏微分，可進(jìn)行類似地推廣。二階偏微分方程的一般形式如下：Axx+Bxy+Cyy=f(x,y,x,y)其中，為彈性體上的某一特征物理量（連續(xù)函數(shù)）。當(dāng)A、B、C都是常數(shù)時(shí)，（1）式稱為準(zhǔn)線性，有三種準(zhǔn)線性方程形式：n 如果=B2-4AC<0，則稱為橢圓型方程；n 如果=B2-4AC=0，則稱為拋物型方程；n 如果=B2-4AC>0，則稱為雙曲型方程。橢圓型方程主要用來處理穩(wěn)態(tài)或靜態(tài)問題，如熱傳導(dǎo)等問題；拋物線方程主要用來處理瞬態(tài)問題，如滲透、擴(kuò)散等問題；雙曲型方程主要用來處

7、理振動(dòng)問題，如玄震動(dòng)、薄膜震動(dòng)等問題。除了上述微分方程外，必須給出定解條件，通常有如下三類：n 第一類邊界條件（Dirichlet條件）：|=(x,y)；n 第二類邊界條件（Neumann條件）：n|=1(x,y)；n 第三類邊界條件（Robin條件）：n+(x,y)|=(x,y)；其中，為求解域的邊界，n為的單位外法矢，(x,y)|0 。第二類和第三類邊界條件統(tǒng)稱為導(dǎo)數(shù)邊界條件。2.3 有限差分方程的數(shù)學(xué)基礎(chǔ)2.3.1 一元函數(shù)導(dǎo)數(shù)的差分公式一個(gè)函數(shù)在x點(diǎn)上的導(dǎo)數(shù)，可以近似地用它所臨近的兩點(diǎn)上的函數(shù)值的差分來表示。函數(shù)fx在x=x0處的泰勒展式如下：fx=n=01n!fn(x0)x-x0n（

8、2） =fx0+f'xx-x0+12!f''xx-x02+13!f(3)xx-x03+（3）對(duì)一個(gè)單變量函數(shù)fx，以步長(zhǎng)x=h將a,b區(qū)間等距劃分，我們得到一系列節(jié)點(diǎn)：x0=a，x1=x0+x=a+h，x2=x1+x=a+2h，xi=xi-1+x=a+ih，xn=b （n=b-ah），然后求出 fx在這些節(jié)點(diǎn)上的近似值。與節(jié)點(diǎn)xi相鄰的節(jié)點(diǎn)有xi-h和xi+h，因此在點(diǎn)xi處可以構(gòu)造如下形式的展開式：fxi-h=fxi-f'xih+12!f''xih2+R2(x)（4）fxi+h=fxi+f'xih+12!f''xih2+

9、R2(x)由式（3）和式（4）可得到：（5）n 一階向前差分：f'xi=fxi+h-fxih（6）n 一階向后差分：f'xi=fxi-fxi-hh（7）n 一階中心差分：f'xi=fxi+h-fxi-h2h不妨，記fi=f(xi)，則式（5）、（6）、（7）分別簡(jiǎn)寫為：（8）n 一階向前差分：fi'=fi+1-fih（9）n 一階向后差分：fi'=fi-fi-1h（10）n 一階中心差分：fi'=fi+1-fi-12h根據(jù)式（8）、式（9）和式（10），可得二階差分：（11）n 二階向前差分：fi''=fi+1'-fi&#

10、39;h=fi+2-2fi+1+fih2（12）n 二階向后差分：fi''=fi'-fi-1'h=fi-2fi-1+fi-2h2（13）n 二階中心差分：fi''=fi+1'-fi-1'2h=fi+2-2fi+fi-24h2差分公式（13）是以相隔2h的兩結(jié)點(diǎn)處的函數(shù)值來表示中間結(jié)點(diǎn)處的一階導(dǎo)數(shù)值，可稱為中點(diǎn)導(dǎo)數(shù)公式。式（11）和式（12）是以相鄰三結(jié)點(diǎn)處的函數(shù)值來表示一個(gè)端點(diǎn)處的一階導(dǎo)數(shù)值，可稱為端點(diǎn)導(dǎo)數(shù)公式。應(yīng)當(dāng)指出：中點(diǎn)導(dǎo)數(shù)公式與端點(diǎn)導(dǎo)數(shù)公式相比，精度較高。因?yàn)榍罢叻从沉私Y(jié)點(diǎn)兩邊的函數(shù)變化，而后者卻只反映了結(jié)點(diǎn)一邊的函數(shù)變化

11、。因此，我們總是盡可能應(yīng)用前者，而只有在無法應(yīng)用前者時(shí)才不得不應(yīng)用后者。（14）但是，由于式（11）中的各階導(dǎo)數(shù)均使用的是向前差分，導(dǎo)致用到的節(jié)點(diǎn)不相鄰，同時(shí)為了均衡誤差，將節(jié)點(diǎn)i處用到的一階差分換成向后差分，則式（11）修正為：fi''=fi+1'-fi'h=fi+1-fih-fi-fi-1hh=fi+1-2fi+fi-1h2同理，根據(jù)上述推導(dǎo)過程，可得到任意階的差分公式：（15）n n階向前差分：fi(n)=fi+1(n-1)-fi(n-1)h（16）n n階向后差分：fi(n)=fi(n-1)-fi-1(n-1)h（17）n n階中心差分：fi(n)=fi

12、+1(n-1)-fi-1(n-1)2h說明，上述公式中各節(jié)點(diǎn)處前一階導(dǎo)數(shù)的代入可能存在不一致，可能是向前差分、向后差分或者中心差分，從而使最終的公式在系數(shù)上存在差別。當(dāng)然，也可以對(duì)各相鄰節(jié)點(diǎn)進(jìn)行需要階數(shù)的泰勒展開，從而建立方程組直接求各階導(dǎo)數(shù)。2.3.2 微分方程轉(zhuǎn)化為線性方程ym（18）由于三種類型的微分方程解法類似，故這里僅以橢圓型微分方程為例，將微分方程轉(zhuǎn)化為代數(shù)方程，對(duì)于雙曲型和拋物型方程依次類推即可。不妨記：2u=uxx+uyy（2稱為拉普拉斯算子），fx,y和g(x,y)是求解域上的連續(xù)函數(shù)。假設(shè)求解區(qū)域?yàn)椋篟=x,y:0xa,0yb,ba=m/n，將求解區(qū)域劃分成(n-1)

13、15;(m-1)個(gè)網(wǎng)格，其中：a=nh,b=mh，如圖1所示。記fi,j=f(xi,yj)，則根據(jù)式（14）可得到：2u=uxx+uyy=ui+1,j-2ui,j+ui-1,jh2+ui,j+1-2ui,j+ui,j-1h2 =ui+1,j+ui-1,j+ui,j+1+ui,j-1-4ui,jh2+O(h2)yj+1yjyj-1y1 x1 xi-1 xi xi+1 xn圖1 五點(diǎn)差分公式式（18）也稱為五點(diǎn)差分公式，同理根據(jù)式（12）和式（13）可分別得到向前差分公式（19）和向后差分公式（20），如圖（2所示）。（19）n 向前差分2u=uxx+uyy=ui+2,j-2ui+1,j+ui,j

14、h2+ui,j+2-2ui,j+1+ui,jh2 =ui+2,j+ui,j+2+2ui,j-2ui+1,j-2ui,j+1h2+O(h2)ymym（20）n 向后差分2u=uxx+uyy=ui,j-2ui-1,j+ui-2,jh2+ui,j-2ui,j-1+ui,j-2h2 =ui-2,j+ui,j-2+2ui,j-2ui-1,j-2ui,j-1h2+Oh2x1xi-2xi-1xixi+1xi+2xixnxnx1y1yj-2yj-1yjyjy1yj+1yj+2 ui-2,j2ui,jui,j+2ui,j+1圖2 向前差分（左）和向后差分（右）-2ui-1,j-2ui,j-1ui,j-2-2ui

15、,j+12ui,j-2ui+1,jui+2,j-4ui,jui,j-1ui+1,jui-1,j 圖3 中心差分、向前差分和向后差分的拉普拉斯算子表示（21）利用中心差分公式（18），由于式（18）在點(diǎn)x,y=(xi,yi)處具有二階精度（Oh2），所以式（18）可近似改寫成下式：2ui,jui+1,j+ui-1,j+ui,j+1+ui,j-1-4ui,jh2根據(jù)橢圓方程的具體形式可以將其分為以下三種形式：n 拉普拉斯（Laplace）方程：2u=0n 泊松（Poison）方程：2u=g(x,y)n 赫耳墨次（Helmholtz）方程：2u+fx,yu=g(x,y)根據(jù)式（21），可建立三種不同

16、形式橢圓方程的代數(shù)方程如下：n 拉普拉斯方程：2u=00=2ui,jui+1,j+ui-1,j+ui,j+1+ui,j-1-4ui,jh2化簡(jiǎn)后得到拉普拉斯方程的計(jì)算公式：（22）ui+1,j+ui-1,j+ui,j+1+ui,j-1-4ui,j=0（23）n 泊松方程：2u=gx,yui+1,j+ui-1,j+ui,j+1+ui,j-1-4ui,j-h2gi,j=0n 赫耳墨次方程：2u+fx,yu=g(x,y)（24）ui+1,j+ui-1,j+ui,j+1+ui,j-1+(h2fi,j-4)ui,j-h2gi,j=02.3.3 建立有限差分方程組根據(jù)式（22）（24）建立方程組，但是需要

17、知道對(duì)應(yīng)的邊界條件才能使方程組存在定解，根據(jù)2.2中可知，邊界條件一般分為狄利克雷邊界條件和導(dǎo)數(shù)邊界條件兩種，下面分別給出這兩種邊界條件的有限差分方程組的建立過程：n 狄利克雷邊界條件：|=(x,y)對(duì)于狄氏條件而言，給出了邊界上各節(jié)點(diǎn)出的函數(shù)計(jì)算公式，直接代入節(jié)點(diǎn)值(xi,yi)計(jì)算即可，如下所示為矩形區(qū)域的邊界點(diǎn)計(jì)算：u(x1,yj)=u1,j=(x1,yj)(1jm)（左邊界）（25）u(xn,yj)=un,j=(xn,yj)(1jm)（右邊界）u(xj,y1)=uj,1=(xj,y1)(1jn)（下邊界）u(xj,yn)=uj,n=(xj,yn)(1jn)（上邊界）n 導(dǎo)數(shù)邊界條件：u

18、(x,y)N=0（26）以右邊界點(diǎn)為例，對(duì)于右邊界點(diǎn)x=xn，根據(jù)Neumann條件可得下式：u(x,y)N=u(xn,yj)x=uxxn,yj=0對(duì)于拉普拉斯方程，根據(jù)計(jì)算公式（22），對(duì)于邊界上的點(diǎn)(xn,yj)可得：（27）un+1,j+un-1,j+un,j+1+un,j-1-4un,j=0（28）顯然，上式中的un+1,j在求解域外，是未知量。根據(jù)中心差分公式（10）可得到：uxxn,yjun+1,j-un-1,j2h根據(jù)式（28）可得到逼近表示：un+1,jun-1,j，并且具有2階逼近精度，代入式（27）可得下式：（29）2un-1,j+un,j+1+un,j-1-4un,j=0

19、同理，對(duì)于其它邊界可獲得如下邊界方程：（30）2ui,2+ui-1,1+ui+1,1-4ui,1=0（下邊界）2ui,m-1+ui-1,m+ui+1,m-4ui,m=0（上邊界）2u2,j+u1,j+1+u1,j-1-4u1,j=0（左邊界）圖4 Neumann條件算子對(duì)于泊松方程和赫耳墨次方程同樣根據(jù)上述方法，獲得邊界條件的線性方程，然后將這些方程添加到式（22）（24）所建立的方程組中，從而建立起(n-1)個(gè)(m-1)元的線性方程組，解該方程組即可獲得各節(jié)點(diǎn)的函數(shù)值。對(duì)于上述過程建立的線性方程組的求解，可采用多種方法，比如Jacob迭代法、Gauss-Seidel迭代法、超松弛迭代法（SO

20、R法）、高斯消元法等方法求解。2.4 有限差分法的收斂性和穩(wěn)定性由于迭代法必須保證收斂性，所以在解有限差分方程組時(shí)還應(yīng)保證其收斂性，也就是通常所說的算法穩(wěn)定性。有限差分法的算法穩(wěn)定性可以通過特征值方法、傅里葉變換（馮諾依曼條件）以及能量估計(jì)等方法來判斷，下面給出常用的馮諾依曼條件：n 向前差分：r1，絕對(duì)收斂；n 向后差分：r>0，絕對(duì)收斂；n 中心差分：對(duì)任何的r對(duì)不收斂；假設(shè)求解域內(nèi)x方向網(wǎng)格劃分的步長(zhǎng)為h，y方向網(wǎng)格劃分的步長(zhǎng)為k，將偏微分方程化為標(biāo)準(zhǔn)形式，具體來說標(biāo)準(zhǔn)形式如下：（31）n 雙曲方程：2uy2=c22ux2(c>0)對(duì)于式（31）所示的雙曲方程，馮諾依曼條件為

21、：r=ckh.（32）n 拋物方程：uy=a2ux2對(duì)于式（32）所示的拋物方程，馮諾依曼條件為：r=ckh2.（33）n 橢圓方程：c22ux2+2uy2=0對(duì)于式（33）所示的橢圓方程，馮諾依曼條件為：r=ckh.（34）為了使算法在任何情況下都能保持穩(wěn)定性，去掉對(duì)網(wǎng)格劃分的馮諾依曼條件，通常采用隱式方案，對(duì)五點(diǎn)差分公式中的節(jié)點(diǎn)所在的行做差分，然后把這些差分的加權(quán)作為中心點(diǎn)的差分值，則拉普拉斯算子可修正為：uxx=ui+1,j+1-2ui,j+1+ui-1,j+1h2+1-2ui+1,j-2ui,j+ui-1,jh2+ui+1,j-1-2ui,j-1+ui-1,j-1h2（01）利用式（3

22、4）進(jìn)行計(jì)算時(shí)，穩(wěn)定性沒有任何限制。取不同的值得到不同的差分公式，通常取=1/4.（35）為了提高計(jì)算精度，很明顯的一個(gè)措施就是網(wǎng)格細(xì)分，但是由于隨著網(wǎng)格步長(zhǎng)的減小，未知量的數(shù)目將會(huì)呈指數(shù)增長(zhǎng)，網(wǎng)格劃分太細(xì)會(huì)導(dǎo)致計(jì)算量過于龐大而無法計(jì)算。通常，我們可以通過提高逼近的精度，采用更高精度的差分公式，例如對(duì)公式（21）進(jìn)行修改，可得到九點(diǎn)差分公式：2ui,j16h2ui+1,j+1+ui-1,j+1+ui+1,j-1+ui-1,j-1+4ui+1,j+4ui-1,j+4ui,j-1+4ui,j+1-20ui,j+O(h4)ui-1,j+14ui-1,jui-1,j-14ui+1,j-20ui,j4u

23、i-1,j4ui,j+1ui-1,j+1ui+1,j+1xi+1xixi-1x1yj+1yj-1yjy1ym圖5 九點(diǎn)差分公式3 有限差分法求解實(shí)例根據(jù)上述推導(dǎo)有限差分法理論，對(duì)于不同類型的偏微分方程建立有限差分方程組，采用mat lab編程給出一些計(jì)算實(shí)例如下：1. 橢圓型方程n 拉普拉斯方程：2u=0；求解域：R=x,y:0x1.5,0y1.5下面分別給出拉普拉斯方程在不同的邊界條件下的解。a) 狄利克雷邊界條件：下邊界：ux,0=x4上邊界：ux,1.5=x4-13.5x2+5.0625左邊界：u0,y=y4右邊界：u1.5,y=y4-13.5y2+5.0625圖6 狄利克雷邊界條件下拉普拉斯方程的解b) Neumann邊界條件：u(x,y)N=0下邊界：ux,0=0上邊界：ux,1.5=400左邊界：u0