劉相國,徐富強
( 巢湖學院 數學系,安徽 巢湖 238000 )
一類拋物型方程的有限差分法
劉相國,徐富強
( 巢湖學院 數學系,安徽 巢湖 238000 )
利用有限差分法求解了拋物型方程邊值問題,得到了相應的穩定性分析,并進行了數值模擬。模擬結果表明該方法是可行的、有效的。
拋物型方程; 有限差分法; 穩定性
現代科學、技術、工程中的大量數學模型都可用偏微分方程來描述,很多近代自然科學的基本方程本身就是偏微分方程。人們一直用微分方程來描述、解釋遇見的各種自然現象,但絕大部分偏微分方程的定解不能以實用的解析形式來表示,這就給一些自然現象的定量描述帶來了很大的困難。隨著計算機的出現和發展,一門新興的學科——微分方程的數值[1][2][3]方法得到了前所未有的發展和應用,并形成了一個專門的數學分支,成為理論和應用研究的熱點課題。本文利用有限差分法研究了熱傳導方程邊值問題,[4]并進行了穩定性分析[5][6][7]和相應的數值模擬,取得了較滿意的數值試驗結果。
設一類最簡單的拋物型方程——熱傳導方程,表示為:

熱傳導方程是初始溫度分布函數為f(x),端點有恒溫c1和c2的絕緣桿上溫度的數學模型。盡管通過傅里葉級數可得到方程的解析解,但這里將它作為求解拋物型方程數值解的一個原型。


設行j的近似值已知,通過式(7)可得到網格中的行j+1的近似值。注意此公式是根據ui?1,j、ui,j和


圖1 在區間R內求解ut(x,t )=c2uxx(x,t)網格

圖2:向前差分空格樣板
差分方程(7)的精度為O(k)+O(h2)。因為項O(k)隨著k趨近于零而線性減小,所以使它的值小可得到更好的近似值。然而,穩定性需求需要考慮更多。設網格的解精度不夠,而且Δx=h0和Δt=k0必須減小。為簡單起見,設新的x增量為Δx=h1=h0/2。如果采用同樣的r,則k1必須滿足:

這使得沿x軸和t軸的網格點的數量分別增加到2倍和4倍。這樣當減少網格大小時,會增加8倍的計算量。這樣的努力是不提倡的,需要開發更有效的無穩定性限制的方法。通過增加方法的復雜性可達到所需的無條件穩定性。
由John Crank和Phyllis Nicholson發明的隱式差分格式是基于求解網格中在行之間點(x,t+k/2)處的式(1)的數值近似解。而且求解ut(x,t+k2)的近似值公式是從中心差分方程中得到,表示為:



其中i=2,3,…,n?1。式(15)右邊的項都是已知的。因此,式(15)可形成對角線性方程組AX=B。式(15)中使用了 6個點,并結合了數值差分基于的中間網格點,如圖3所示,有時通過使r=1來實現式(15)。在這種情況下,沿t軸的增量為Δt=k=h2/c2,同時式(15)可簡化為:

當Crank-Nicholson法用計算機實現時,可通過直接解法或迭代法得到線性方程組AX=B。

圖3:Crank-Nicholson法的空格樣板
拋物型方程的差分格式在實際應用時,都取逐層計算的形式。當初始層上的解有誤差時,必須逐層傳播,影響以后各層的解,研究這個誤差傳播的規律是穩定性討論的任務。下面以齊次邊界條件熱傳導方程(1)、(2)、(3)的差分格式為模型,進行敘述和討論:



用向前差分法以及C-N格式求解熱傳導方程:

例1用向前差分法求解(22)如下:
(1)當r=0.1、t=0.1、h=0.1(m=10、N=100)時取部分節點上的運行結果(表1)及數值解剖分圖(圖4)如下:

表 1 當r=0.1、t=0.1、h=0.1(m=10、N=100)時部分節點上的運行結果
(2)當r=0.5、t=0.1、h=0.1(m=10、N=20)時取部分節點上的運行結果(表2)及數值解剖分圖(圖5)如下:

表2 當r=0.5、t=0.1、h=0.1(m=10、N=20)時部分節點上的運行結果

圖4 數值解剖分圖

圖5 數值解剖分圖
例2用Crank-Nicholson法求解(22)如下:
<1>當r=0.1、t=0.1、h=0.1(m=10、N=100)時,取部分節點上的運行結果,見表3。

表 3 當r=0.1、t=0.1、h=0.1(m=10、N=100)時部分節點上的運行結果

表4 當r=0.5、t=0.1、h=0.1(m=10、N=100)時部分節點上的運行結果
本文利用有限差分法研究了一類拋物型偏微分方程邊值問題,并建立了相應的穩定性分析。通過數值模擬可以看出,有限差分法在求解此類問題是有效性和可行的。并且該法求得的數值解有較高的精度。
[1] Carlos Castro, Sorin Micu. Boundary controllability of a linear semi-discrete 1-D wave equation derived from a mixed finite element method[J]. Numer. Math. (2006) 102:413–462.
[2] Partha Roy Chaudhuri. Sourabh Roy. Analysis of arbitrary index profile planar optical waveguides and multilayer nonlinear structures: a simple finite difference algorithm[J].Opt. Quant. Electron. (2007) 39: 221–237.
[3] LUO Yun-ju, LIU Dong-yan, LIU Xin-rong. Finite element numerical simulation for the hydrodynamic field evolution of geothermal water in the nanwenquan anticline in chongqing in china[J]. Journal of Hydrodynamics, 2006,18(4): 443-448.
[4] 周順興.解拋物型偏分方程的高精度差分格式[J].計算數學,1982,4(2):204-203.
[5] 金承日.解拋物型方程的高精度顯式差分格式[J].計算數學,1991,(1): 38-44.
[6] 李榮華.偏微分方程數值解法[M].北京:高等教育出版社,2005.
[7] 范德輝等.對流擴散方程差分格式穩定性分析[J].暨南大學學報(自然科學版),2006,27(1):24-29.
A Finite Difference Method for Solving Parabolic Equation
LIU Xiang-guo, XU Fu-qiang
( Department of Mathematics, Chaohu University, Chaohu, Anhui 238000, China )
This paper uses finite difference method to solve the boundary value problem of parabolic equation,obtaining the corresponding stability analysis and numerical simulation. Simulation results show that the method is feasible and effective.
parabolic equation;finite difference;stability
(責任編輯 李宗寶)
H319 < class="emphasis_bold">文獻標識碼:A
A
1673-9639 (2011) 04-0139-05
2011-01-03
巢湖學院科研基金項目(編號:XLY-201006)。
劉相國(1979-),男,寧夏人,巢湖學院數學系教師。
徐富強(1984-),安徽宣城人,研究方向:智能計算。