擴(kuò)散方程的差分解法

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

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

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

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

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

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

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

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

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

10、對(duì)角方程組，可采用追趕法進(jìn)行求解。5.3 計(jì)算程序 5.3.1顯格式!-顯格式求解擴(kuò)散方程-!-參數(shù)設(shè)置-！parameter (nt=20001,nx=101) !nt:時(shí)間節(jié)點(diǎn)數(shù);nx:空間節(jié)點(diǎn)數(shù)dimension u(nx,nt)Sx=10.0!計(jì)算長(zhǎng)度（m）dt=0.003!時(shí)間步長(zhǎng)（s）dx=0.1!空間步長(zhǎng)（m）alfa=1.0!擴(kuò)散系數(shù)!-!open(1,file='顯格式沿程變化.txt')open(2,file='顯格式隨時(shí)間變化.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輸出隨時(shí)間變化-!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全隱格式!-全隱格式求解擴(kuò)散方程-!-參數(shù)設(shè)置-！!implicit double precision (a-

13、h,o-z)parameter (nt=20001,nx=101) !nt:時(shí)間節(jié)點(diǎn)數(shù);nx:空間節(jié)點(diǎn)數(shù)dimension u(nx,nt),a(nx),b(nx),c(nx-1),d(nx-1),xx(nx)Sx=10!計(jì)算長(zhǎng)度（m）dt=0.003!時(shí)間步長(zhǎng)（s）dx=0.1!空間步長(zhǎng)（m）alfa=1.0!擴(kuò)散系數(shù)!-!open(1,file='全隱格式沿程變化.txt')open(2,file='全隱格式隨時(shí)間變化.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!-!-三對(duì)角方程組系數(shù)矩陣-!a(1)=1do i=1,nx-1 d(i)=-r c(i)=-r a(i+1)=1+2*renddoc(1)=0!-!-三對(duì)角方程組常數(shù)項(xiàng)-!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輸出隨時(shí)間變化-!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：下對(duì)角線! c：上對(duì)角線! a：主對(duì)角線! b：方程右邊常數(shù)項(xiàng)! 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 計(jì)算結(jié)果與分析6.1顯格式計(jì)算結(jié)果與分析 6.1.1隨時(shí)間變化 經(jīng)計(jì)算，沿程斷面u隨時(shí)間變化如下圖5-圖7所示，從圖中可以看出，各斷面u值隨著時(shí)間逐漸遞減，在t=60s時(shí)，u值衰減至接近于0，曲線的斜率隨著時(shí)間逐漸趨緩，可見(jiàn)u的衰減速率逐漸減小。由初始條件可知，初始時(shí)刻，u分布關(guān)于x=5m對(duì)稱，對(duì)比x=2.5m及x=7.5m的u值隨時(shí)間變化可知，這兩個(gè)斷面的u隨時(shí)間變化關(guān)系一致，說(shuō)明對(duì)

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

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

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