1、擴散方程的差分解法在研究熱傳導(dǎo)過程、擴散過程、邊界層現(xiàn)象時，我們常常遇到拋物型方程，這類方程中最典型、最簡單的就是熱傳導(dǎo)方程。熱傳導(dǎo)方程中的自變量中包括時間t，它是描述一種隨時間變化的物理過程，即所謂不定?，F(xiàn)象。這類問題的基本定解問題應(yīng)是初值問題，即在初始時刻（t=0）時給定定解條件，求解t>0時的解。本文主要運用有限差分法對一維擴散方程進行求解，并對差分解的適定性、相容性、收斂性及穩(wěn)定性進行分析，同時與解析解進行對比。1.擴散方程一維擴散方程為：（1）式中，u為因知量，為擴散系數(shù)，x為坐標(biāo)，t為時間。其定解條件如下： 初始條件：（2） 邊界條件：（3）一般假定函數(shù)，滿足連接條件，即，。

2、2.有限差分法有限差分法是數(shù)值計算解微分方程古老的方法之一，也是系統(tǒng)化地、數(shù)值地求解數(shù)學(xué)物理方法的方程。其控制方程中的導(dǎo)數(shù)用離散點上函數(shù)值的差商代替。差分格式可以分為顯格式和隱格式。所謂顯格式是指在任一結(jié)點上因變量在新是時間層上的值可以通過之前的時間層上相鄰結(jié)點變量的值顯式解出來。由于這些層的變量值是已知的，當(dāng)時間向前推進時，空間點上的新的變量值就只需逐點計算就行了，因此顯格式計算起來比較省事。隱格式則是指任一結(jié)點上變量在新的時間層的值，不能通過之前的時間層上相鄰結(jié)點的值顯式解出來，它不僅與之前的時間層上的已知值有關(guān)，而且也與新時間層的相鄰結(jié)點的變量值有關(guān)。因而一個差分方程常常包括幾個相鄰結(jié)點

3、上的未知數(shù)，未知數(shù)的個數(shù)取決于格式的構(gòu)成形式。為了解出這些未知數(shù)需要聯(lián)立新的方程，而每引進一個新的方程往往又同時引進了新的未知數(shù)。因此，隱格式總是伴隨著求解巨大的代數(shù)方程組。隱格式的主要缺點是計算工作量大，因而不如顯格式計算得快，但這只是就時間步長一樣的情況而言的。隱格式的主要優(yōu)點是時間步長可以比顯格式能夠采用的最大步長大很多。顯格式的時間步長受到穩(wěn)定性條件的限制，而隱格式則幾乎不受限制。3.方程的離散3.1 顯格式采用時間前差及第n時間層的空間中心差，得一維擴散方程的顯格式解：（4）即（5）式中，3.2 全隱格式采用時間前差及第（n+1）時間層的空間中心差，得一維擴散方程的全隱格式解：（6）

4、4.差分解的基本問題 差分解的基本問題包括：適定性、相容性、收斂性和穩(wěn)定性四個方面。4.1 適定性 在用差分方程作微分方程數(shù)值解時，首先，要求微分方程的問題是適定的。所謂適定性問題是指這一微分方程在一定的初始和邊界條件下要有唯一解，并且在初始條件和邊界條件稍有改變時，微分方程的解也只是稍有偏離。從數(shù)學(xué)的角度講，若微分方程的解存在并且是唯一的，同時連續(xù)依賴于數(shù)據(jù)（初始條件、邊界條件），則問題是適定的。對于拋物線方程，這是初值問題，要求給出初始條件，并且在區(qū)域邊界上給出邊界條件。在本文的一維擴散問題中，除了給出t=0時的初值外，應(yīng)在左、右兩端邊界，即及，各給定一個邊界條件，第一類邊界條件u的值，或

5、給定第二類邊界條件的值或給定第三類邊界條件u與的組合。 由于拋物型方程具有擴散性質(zhì)，它往往是隨時間擴散的，u值將不斷因擴散而衰減。 在本文的一維擴散方程中，給定了初始條件，同時在區(qū)域的左、右端邊界給定了邊界條件，滿足適定性要求。4.2 相容性 將一個偏微分方程用差分格式化為相應(yīng)差分方程，當(dāng)步長和趨近于零時，這個差分方程應(yīng)當(dāng)收斂于原微分方程，也就是說，相應(yīng)的差分方程和微分方程之間的截斷誤差在任一時刻任一網(wǎng)格點上均應(yīng)趨近于零，這樣的差分方程和微分方程才是相容的。 對一維擴散方程：（7）采用顯格式差分格式，令，則（8）用Taylor級數(shù)展開代入上述差分方程中，則有（9）當(dāng)，時，上式化為（10）可見，

6、此差分格式所構(gòu)成的差分方程與原來的微分方程是相容的，故該顯格式為相容格式。采用全隱格式差分格式，令，則（11）用Taylor級數(shù)展開代入上述差分方程中，則有（12）當(dāng)，時，上式化為（13）可見，此差分格式所構(gòu)成的差分方程與原來的微分方程是相容的，故該全隱格式為相容格式。4.3 穩(wěn)定性 差分法計算中所產(chǎn)生的誤差（舍入誤差，參數(shù)誤差等）隨時間衰減或不增大，則稱離散格式是穩(wěn)定的，反之，則是不穩(wěn)定的。分析差分方程穩(wěn)定性有不同的方法，如矩陣方法，諧波分析法等。下面用諧波分析法對以上兩種離散格式的穩(wěn)定性進行分析：（1）顯格式顯格式的差分方程為（14）即（15）其誤差方程為（16）任取一k次諧波分量（17）

7、則（18）誤差放大因子為（19）要滿足穩(wěn)定性條件，則要求對所有的k值均有，須，即。因此，一維擴散方程顯格式的穩(wěn)定性條件為（20）（2）全隱格式全隱格式的差分方程為（21）即（22）其誤差方程為（23）任取一k次諧波分量（24）則，（25）則誤差方程為（26）誤差放大因子為（27）要滿足穩(wěn)定性條件，則要求對所有的k值均有。從（28）式中可以看出，當(dāng)（即）時，恒成立。因此，全隱格式是無條件穩(wěn)定的。4.4 收斂性 如果差分方程的解為，微分方程的解為，若當(dāng)，時，差分方程的解與微分方程的解之差（28）則稱差分格式是收斂的。 拉克斯（Lax）等價定理指出，如果問題是適定的，并且差分格式滿足相容性條件，那么

8、差分格式的穩(wěn)定性就是該格式收斂性的充分而必要的條件。 由4.1、4.2和4.3的分析可知，顯格式和全隱格式都滿足適定性、相容性及穩(wěn)定性的條件，因而這兩種格式滿足收斂性要求。5.差分方程的求解5.1 初始條件及邊界條件由以上對一維擴散問題的分析，可知，求解一維擴散方程需給定初始條件及邊界條件。在本文計算中，取，。初始條件（時）（29）邊界條件為（30）其初始時刻（）時的u分布如圖1所示，x=0m處u隨時間變化情況如圖2所示，x=10m處u隨時間變化情況如圖3所示。圖1 初始時刻u分布圖圖2 x=0處u隨時間變化圖圖3 x=10m處u隨時間變化圖5.2 空間步長及時間步長 通過對差分格式的穩(wěn)定性分

9、析知，顯格式空間步長與時間步長間應(yīng)滿足一定的關(guān)系，即式（21）。本文選用的步長如表1所示，分別用顯格式和全隱格式進行數(shù)值計算，差分網(wǎng)格如圖4所示。表1 時間步長及空間步長取值表時間步長/s空間步長/m時間步數(shù)計算時長/s空間步數(shù)0.0030.150000015001000.3圖4 差分網(wǎng)格示意圖5.2 求解方法 對于顯格式，有，若n時間層上的u已知，則可直接求解出（n+1）時間層上的u值。因此，由于給出了初始條件，顯格式無需迭代，直接從t=0時刻開始，逐一計算下一個時間層上的u值，便可求解出各時間層上各空間點的u值。對于全隱格式，由于與，有聯(lián)系，不能直接求解，必須聯(lián)解代數(shù)方程組。此時方程組為三

10、對角方程組，可采用追趕法進行求解。5.3 計算程序 5.3.1顯格式!-顯格式求解擴散方程-!-參數(shù)設(shè)置-！parameter (nt=20001,nx=101) !nt:時間節(jié)點數(shù);nx:空間節(jié)點數(shù)dimension u(nx,nt)Sx=10.0!計算長度（m）dt=0.003!時間步長（s）dx=0.1!空間步長（m）alfa=1.0!擴散系數(shù)!-!open(1,file='顯格式沿程變化.txt')open(2,file='顯格式隨時間變化.txt')r=alfa*dt/dx*2!-邊界條件-!do n=1,nt u(1,n)=0 u(nx,n)=0en

11、ddo!-!-初始條件-!do j=1,nx x=(j-1)*dx if (x.le.sx/2) then u(j,1)=2*x /sx else u(j,1)=2*(1-x/sx) endifenddo!-!-求解差分方程-!do n=2,nt do j=2,nx-1 u(j,n)=u(j,n-1)+r*(u(j-1,n-1)-2*u(j,n-1)+u(j+1,n-1) enddoenddo!-!-每隔15s輸出沿程變化-!do n=1,nt t=(n-1)*dt if (mod(int(t*1000),15000).eq.0) then write(1,*) 't=',t,

12、's' do j=1,nx write(1,*) (j-1)*dx,u(j,n) enddo endifenddo!-!-每隔2.5m輸出隨時間變化-!do j=1,nx x=dx*(j-1) if (mod(int(x*10),25).eq.0) then write(2,*) 'x=',x,'m' do n=1,nt,200 write(2,*) (n-1)*dt,u(j,n) enddo endifenddo!-!end 5.3.2全隱格式!-全隱格式求解擴散方程-!-參數(shù)設(shè)置-！!implicit double precision (a-

13、h,o-z)parameter (nt=20001,nx=101) !nt:時間節(jié)點數(shù);nx:空間節(jié)點數(shù)dimension u(nx,nt),a(nx),b(nx),c(nx-1),d(nx-1),xx(nx)Sx=10!計算長度（m）dt=0.003!時間步長（s）dx=0.1!空間步長（m）alfa=1.0!擴散系數(shù)!-!open(1,file='全隱格式沿程變化.txt')open(2,file='全隱格式隨時間變化.txt')r=alfa*dt/dx*2!-邊界條件-!do n=1,nt u(1,n)=0 u(nx,n)=0enddo!-!-初始條件-!

14、do j=1,nx x=(j-1)*dx if (x.le.sx/2) then u(j,1)=2*x/sx else u(j,1)=2*(1-x/sx) endifenddo!-!-三對角方程組系數(shù)矩陣-!a(1)=1do i=1,nx-1 d(i)=-r c(i)=-r a(i+1)=1+2*renddoc(1)=0!-!-三對角方程組常數(shù)項-!do i=1,nx b(i)=u(i,1)enddo!-!-求解差分方程-!do n=2,nt call systri(nx,d,a,c,b,xx) do j=1,nx-1 b(j)=xx(j) u(j,n)=xx(j) enddoenddo!-!

15、-每隔15s輸出沿程變化-!do n=1,nt t=(n-1)*dt if (mod(int(t*1000),15000).eq.0) then write(1,*) 't=',t,'s' do j=1,nx write(1,*) (j-1)*dx,u(j,n) enddo endifenddo!-!-每隔2.5m輸出隨時間變化-!do j=1,nx x=dx*(j-1) if (mod(int(x*10),25).eq.0) then write(2,*) 'x=',x,'m' do n=1,nt,200 write(2,*)

16、(n-1)*dt,u(j,n) enddo endifenddo!-!end subroutine systri(n,d,a,c,b,x)! d：下對角線! c：上對角線! a：主對角線! b：方程右邊常數(shù)項! x：方程的解! implicit double precision (a-h,o-z) dimension a(n),d(n-1),c(n-1) dimension p(n),q(n-1),y(n),x(n),b(n) p(1)=a(1) do i=1,n-1 q(i)=c(i)/p(i) p(i+1)=a(i+1)-d(i)*q(i) enddo y(1)=b(1)/p(1) do

17、i=2,n y(i)=(b(i)-d(i-1)*y(i-1)/p(i) enddo x(n)=y(n) do i=n-1,1,-1 x(i)=y(i)-q(i)*x(i+1) enddo return end6 計算結(jié)果與分析6.1顯格式計算結(jié)果與分析 6.1.1隨時間變化 經(jīng)計算，沿程斷面u隨時間變化如下圖5-圖7所示，從圖中可以看出，各斷面u值隨著時間逐漸遞減，在t=60s時，u值衰減至接近于0，曲線的斜率隨著時間逐漸趨緩，可見u的衰減速率逐漸減小。由初始條件可知，初始時刻，u分布關(guān)于x=5m對稱，對比x=2.5m及x=7.5m的u值隨時間變化可知，這兩個斷面的u隨時間變化關(guān)系一致，說明對

18、稱分布的初始條件隨著時間的變化，其分布仍具有對稱性。圖5 x=2.5m處u隨時間變化圖6 x=5m處u隨時間變化圖7 x=7.5m處u隨時間變化 6.1.2沿程分布圖8-圖11給出了t=15s、t=30s、t=45s及t=60s時u值的分布情況，從圖中可以看出，u值的分布呈拋物線分布，在x=5m處達到最大值，同時u分布關(guān)于x=5m對稱?？梢?，u分布仍保持初始時刻的對稱關(guān)系，只是分布曲線趨于平緩。對比各時刻u的峰值，可知，u的峰值隨著時間逐漸減小，至t=60s時，其峰值為0.002左右，此時最大值與最小值之差在0.002左右，可見，在t=60s時，擴散趨于穩(wěn)定，即整個區(qū)域u值基本相同。圖8 t=

19、15s時u沿程分布圖9 t=30s時u沿程分布圖10 t=45s時u沿程分布圖11 t=60s時u沿程分布 6.2全隱格式計算結(jié)果與分析 6.2.1隨時間變化 經(jīng)計算，沿程斷面u隨時間變化如下圖12-圖14所示，從圖中可以看出，各斷面u值隨著時間逐漸遞減，在t=60s時，u值衰減至接近于0，曲線的斜率隨著時間逐漸趨緩，可見u的衰減速率逐漸減小。由初始條件可知，初始時刻，u分布關(guān)于x=5m對稱，對比x=2.5m及x=7.5m的u值隨時間變化可知，這兩個斷面的u隨時間變化關(guān)系一致，說明對稱分布的初始條件隨著時間的變化，其分布仍具有對稱性。圖12 x=2.5m處u隨時間變化圖13 x=5.0m處u隨

20、時間變化圖14 x=7.5m處u隨時間變化 6.2.2沿程分布圖15-圖18給出了t=15s、t=30s、t=45s及t=60s時u值的分布情況，從圖中可以看出，u值的分布近似呈拋物線分布，在x=5m處達到最大值，同時u分布關(guān)于x=5m對稱?？梢?，u分布仍保持初始時刻的對稱關(guān)系，只是分布曲線趨于平緩。對比各時刻u的峰值，可知，u的峰值隨著時間逐漸減小，至t=60s時，其峰值為0.002左右，此時最大值與最小值之差在0.002左右，可見，在t=60s時，擴散趨于穩(wěn)定，即整個區(qū)域u值基本相同。圖15 t=15s時u沿程分布圖16 t=30s時u沿程分布圖17 t=45s時u沿程分布圖18 t=60s時u沿程分布6.3計算結(jié)果對比 6.3.1隨時間變化對比 圖19-圖21給出了x=2.5m，x=5m，x=7.5m處顯格式與全隱格式計算所得u值隨時間變化的對比圖，從圖中可以看出，顯格式與全隱格式計算得到的曲線幾乎重合基本一致，說明兩者計算結(jié)