104753150634吳麗煥用有限差分格式解拋物型方程

1、用有限差分格式解拋物型方程 姓名：吳麗煥 學(xué)號：104753150634 專業(yè)：應(yīng)用數(shù)學(xué)用有限差分格式解拋物型方程摘要本文以偏微分方程在初值問題的數(shù)值解法為核心研究內(nèi)容，主要研究了向前差分格式，向后差分格式，六點(diǎn)對稱格式以及格式的構(gòu)造，并用不同方法對這幾種不同格式的穩(wěn)定性進(jìn)行了判斷。關(guān)鍵詞：微分方程數(shù)值解 有限差分法 穩(wěn)定性1前言 微分方程是數(shù)學(xué)的一個重要分支，同樣也是微積分、變分法、控制論、復(fù)變函數(shù)、組合數(shù)學(xué)、組合拓?fù)涞葘W(xué)科的基礎(chǔ)。如果一個微分方程中出現(xiàn)的未知數(shù)只含有一個自變量，這個方程叫做常微分方程；如果一個微分方程中出現(xiàn)多元函數(shù)的偏導(dǎo)數(shù)，或者說未知函數(shù)和幾個變量有關(guān)，而且方程中出現(xiàn)未知數(shù)

2、對幾個變量的導(dǎo)數(shù)，那么這種方程就是偏微分方程。求偏微分方程的定解問題可以先求出它的通解，然后再用定解條件確定出函數(shù)。但一般來說，在實(shí)際中通解是不易求解的，用定解條件確定函數(shù)更是比較困難的。應(yīng)該指出偏微分方程的定解雖然各有解法，但是我們不能忽視由于某些原因有許多定解條件是不能嚴(yán)格解出的，只可以用近似方法求出滿足實(shí)際需要的近似程度的近似解。常用的方法有變分法和有限差分法。有限差分法是把定解問題轉(zhuǎn)化成代數(shù)方程，然后用計算機(jī)進(jìn)行計算。2有限差分法簡介 有限元方法最早應(yīng)用于結(jié)構(gòu)力學(xué)，后來隨著計算機(jī)的發(fā)展慢慢用于流體力學(xué)的數(shù)值模擬。在有限元方法中，把計算域離散剖分為有限個互不重疊且相互連接的單元，在每個單

3、元內(nèi)選擇基函數(shù)，用單元基函數(shù)的線性組合來逼近單元中的真解，整個計算域上總體的基函數(shù)可以看為由每個單元基函數(shù)組成的，則整個計算域內(nèi)的解可以看作是有所有單元上的近似解構(gòu)成。 有限差分方法是計算機(jī)數(shù)值模擬最早采用的方法，至今仍被廣泛使用。該方法將求解域劃分為差分網(wǎng)格，用有限個網(wǎng)格節(jié)點(diǎn)代替連續(xù)的求解域。有限差分法以Taylor級數(shù)展開等方法，把控制方程中的導(dǎo)數(shù)用網(wǎng)格節(jié)點(diǎn)上的函數(shù)值的差商代替進(jìn)行離散，從而建立以網(wǎng)格節(jié)點(diǎn)上的值為未知數(shù)的代數(shù)方程組。該方法是一種直接將微分問題變?yōu)榇鷶?shù)問題的近似數(shù)值解法，數(shù)學(xué)概念直觀，表達(dá)簡單，是發(fā)展較早且比較成熟的數(shù)值方法。 對于有限差分格式，從格式的精度來劃分，有一階格式

4、、二階格式和高階格式。 從差分的空間形式來考慮，可分為中心格式和迎風(fēng)格式。 考慮時間因子的影響，差分格式還可以分為顯格式、隱格式、顯隱交替格式等。 目前常見的差分格式，主要是上述幾種形式的組合，不同的組合構(gòu)成不同的差分格式。差分方法主要適用于有結(jié)構(gòu)網(wǎng)格，網(wǎng)格的步長一般根據(jù)實(shí)際地形的情況和柯朗穩(wěn)定條件來決定。 構(gòu)造差分的方法有多種形式，目前主要采用的是泰勒級數(shù)展開方法。其基本的差分表達(dá)式有三種：一階向前差分、一階向后差分、一階中心差分和二階中心差分等，其中前兩種格式為一階計算精度，后兩種格式為二階計算精度。通過對時間和空間這幾種不同差分格式的組合，可以組合成不同的差分計算格式。 有限元方法的基礎(chǔ)

5、是變分原理和加權(quán)余量法，其基本求解思想是把計算域劃分為有限個互不重疊的單元，在每個單元內(nèi)，選擇一些合適的節(jié)點(diǎn)作為求解函數(shù)的插值點(diǎn)，將微分方程中的變量改寫為由各個變量或其導(dǎo)數(shù)的節(jié)值點(diǎn)與所選用的插值函數(shù)組成的線性表達(dá)式，借助于變分原理或加權(quán)余量法，將微分方程離散求解，采用不同的權(quán)函數(shù)和插值函數(shù)形式，便構(gòu)成不同的有限元方法。 在采用數(shù)值計算方法求解偏微分方程時，若將每一處導(dǎo)數(shù)由有限差分近似公式代替，從而把求解偏微分方程的問題轉(zhuǎn)換成求解代數(shù)方程的問題，即所謂的有限差分法。有限差分法求解偏微分方程的步驟如下：1、 區(qū)域離散化，即把所給偏微分方程的求解區(qū)域細(xì)分為由有限個點(diǎn)組成的網(wǎng)格；2、 近似代替，即采用

6、有限差分法代替每一個格點(diǎn)的導(dǎo)數(shù)；3、 逼近求解。換而言之，這一過程可以看作是用一個插值多項式及其微分來替代偏微分方程的解的過程。3差分逼近的基本概念3.1 差分方程 以最簡單的以為對流問題為例，引入用差分方法求偏微分方程數(shù)值解的一些概念?？紤]對流方程的初值問題 （3.1） 區(qū)域的剖分： 網(wǎng)格剖分可以采用兩組平行于x軸和t軸的直線形成的網(wǎng)覆蓋區(qū)域，它們的交點(diǎn)稱為網(wǎng)格節(jié)點(diǎn)節(jié)點(diǎn)記為。間距稱為空間步長，間距稱為時間步長。 Taylor公式：設(shè)在的某個鄰域內(nèi)具有直到階的導(dǎo)數(shù)，則有是余項，且 微商與差商之間的關(guān)系式 （向前差商） （3.2） （向前差商） （3.3） （向后差商） （3.4） （中心差商）

7、 （3.5）由于u是方程（3.1）的解，所以 （3.6）由（3.2）（3.3），得 （3.7）用方程 （3.8）近似代替，其中表示的近似值。將（3.8）改寫成便于計算的形式 （3.9）這里稱為網(wǎng)格比。（3.8），（3.9）稱為（3.1）的（有限）差分方程（差分格式）。問題（3.1）中的初始條件的離散形式是 （3.10）初始問題（3.1）的差分格式 （3.11）對同一微分方程可以建立不同形式的差分格式。在（3.1）中u對t采用向前差商，u對x采用向后差商和中心差商得 （3.12） （3.13）積分插值法：Green公式：設(shè)閉區(qū)域D由分段光滑的曲線L圍成，函數(shù)及在上有一段連續(xù)偏導(dǎo)數(shù)，則有，其中L是

8、D的取正向的邊界曲線。其中為有向曲線弧L上處的切向量的方向角。 在平面上取矩形域為積分區(qū)域，是D的邊界。將方程（3.1）在D上積分，得到。利用格林公式得 （3.14）其中與分別是L的外法向單位向量n沿x方向與沿t方向的兩個分量。把（3.14）左端分成在上的四個積分，得近似方程，即 （3.15）這里是與的長度，是與的長度，是可按不同方式確定的u在上的近似函數(shù)值。在網(wǎng)格中，點(diǎn)E,F,G,H依次為，并取于是，從（3.15）得到 （3.16）稱為蛙跳格式?，F(xiàn)在換一種方式，在網(wǎng)格中點(diǎn)E,F,G,H依次為并取，于是，從（3.15）得到，也可寫成 （3.17）這個格式成為格式。3.2 截斷誤差 對于齊次問題

9、，可以將微分方程和差分方程記為 其中L是微分算子 其中是相應(yīng)的差分算子方程（3.1）微分算子L為格式（3.8）相應(yīng)的差分算子為設(shè)u時所討論的微分方程的充分光滑的解，將算子L和分別作用于u，記兩者在任意的節(jié)點(diǎn)處的差為E，即，差分格式的截斷誤差是指對E的估計。3.3收斂性設(shè)u是微分方程的準(zhǔn)確解，是相應(yīng)差分方程的數(shù)值解，如果當(dāng)步長時，對任何有則稱差分格式是穩(wěn)定的。4差分格式的構(gòu)造4.1 最簡差分格式考慮一維熱傳導(dǎo)方程 （4.1）其中是常數(shù)，是給定的連續(xù)函數(shù)。 可將（4.1）的定解問題分為兩類：第1、 初值問題（也稱問題）：求具有所需次數(shù)偏微商的函數(shù)，滿足方程（4.1）和初始條件 （4.2）第2、 初

10、邊值問題（也稱混合問題）：具有所需次數(shù)偏微商的函數(shù)，滿足方程（4.1），初始條件邊值條件：假定和在相應(yīng)區(qū)域光滑，并且在滿足相容條件，使上述問題有唯一充分光滑的解。剖分：取空間步長和時間步長，其中都是自然數(shù)。用兩族平行直線和將矩形域分割成為矩形網(wǎng)格，網(wǎng)格節(jié)點(diǎn)為。以表示網(wǎng)格內(nèi)點(diǎn)集合，即位于開矩形G的網(wǎng)點(diǎn)集合；表示所有位于閉矩形的網(wǎng)點(diǎn)集合；是網(wǎng)點(diǎn)界點(diǎn)集合。其次，用表示定義在網(wǎng)點(diǎn)的函數(shù)，。用適當(dāng)?shù)牟罘执娣匠蹋?.1）中相應(yīng)的偏微商，便得到以下幾種最簡差分格式。（1） 向前差分格式，即其中.以表示網(wǎng)比。將改寫成便于計算的形式，使得第k層值u（上標(biāo)為k）在等式右邊，第k+1層值在等式左邊，則得 截斷誤差

11、 記 截斷誤差 （4.5）（2） 向后差分格式其中即： 記截斷誤差 （4.7）（3） 六點(diǎn)對稱格式（格式）將向前差分格式和向后差分格式作算術(shù)平均值，即得六點(diǎn)對稱格式： 將改寫為 令截斷誤差于展開，則得 （4） 格式，即 或 截斷誤差4.2 差分格式的穩(wěn)定性與收斂性 前述引進(jìn)的二層格式均可用矩陣和向量的記號表成 （4.2.1）其中和，AB是矩陣，假定A有逆，并令 （4.2.2）則可將（4.2.1）化為 （4.2.3）先按初值穩(wěn)定，此時F=0， （4.2.8）我們說差分格式（4.2.1）按初值穩(wěn)定，如果存在和常數(shù)，使不等式對一切和成立，這里是中的一種范數(shù)，一般取。顯然差分格式（4.2.1）按初值穩(wěn)

12、定，當(dāng)且僅當(dāng) （4.2.10）其次討論按右端穩(wěn)定，此時認(rèn)為初值沒有誤差即。我們說差分格式（4.2.1）按右端穩(wěn)定，如果存在和常數(shù)使不等式對一切和成立，其中是下列方程的解 （4.2.11）反復(fù)利用遞推式（4.2.11），得設(shè)又差分格式按初值穩(wěn)定，即存在常數(shù)使，則，取，即知格式按右端穩(wěn)定。總之，若，則由格式按初值穩(wěn)定可推出它按右端穩(wěn)定，為檢驗格式按初值穩(wěn)定，需要檢驗不等式（4.2.10），即矩陣族 （4.2.12）一致有界。4.2.1 判別穩(wěn)定性的直接法 命題4.2.1（必要條件）以表示矩陣的譜半徑，則差分格式穩(wěn)定的必要條件是存在與無關(guān)的常數(shù)M使 （4.2.13） 命題4.2.2（充分條件）若是正

13、規(guī)矩陣，則（4.2.13）也是差分格式穩(wěn)定的充分條件。 命題4.2.3 若S是對稱矩陣，是矩陣S的實(shí)系數(shù)有理函數(shù)：，則差分格式穩(wěn)定的充要條件是，其中是S的特征值。 對于向前差分格式當(dāng)時穩(wěn)定，當(dāng)時不穩(wěn)定；對于向后差分格式，對于任何穩(wěn)定，即絕對穩(wěn)定；對于六點(diǎn)對稱格式，對任何有穩(wěn)定，因此絕對穩(wěn)定；對于格式，對任意不穩(wěn)定。4.2.2 Fourier方法 考慮線性常系數(shù)一維拋物方程，具有初值和周期（設(shè)周期為）邊值條件，逼近它的二層差分方程的一般形式為 （4.2.14）（只考慮按初值穩(wěn)定，故可設(shè)非齊次等于0），這里在空間處的差分方程，和是包含0及其附近的正負(fù)整數(shù)的有限集合，和不依賴于但和有關(guān)。 由于是周期

14、邊值條件，故可將周期開拓使其對一切有意義，且方程（4.2.14）對所有整數(shù)成立，為了應(yīng)用Fourier方法，再將開拓為上的，為此取半整點(diǎn)，并用如下階梯函數(shù)逼近初始函數(shù)當(dāng)：其中再將（4.2.14）看成再任一成立，則得具連續(xù)變量的差分解。顯然仍是x的周期，且當(dāng)，此外將Fourier方法用于具連續(xù)變量的差分方程 （4.2.15）將展成Fourier級數(shù)： （4.2.16） （4.2.17）由等式 （4.2.18）把（4.2.16）代到（4.2.15）得 （4.2.19）其中 （4.2.20）比較對應(yīng)的系數(shù)，得 （4.2.21）其中 （4.2.22）將（4.2.21）代到（4.2.19），則 （4.2

15、.23）若差分格式穩(wěn)定，則有常數(shù)使由于階梯函數(shù)類于稠密，故由上式可得 （4.2.24）則一致有界。反之，若一致有界，則由（4.2.23）得從而差分格式按初值穩(wěn)定。稱為增長因子。不等式（4.2.23）又等價于 （4.2.25）上式也稱為條件。命題4.2.4 差分格式（4.2.14）穩(wěn)定一致有界成立。命題4.2.5 差分格式，（其中是方陣，一般依賴步長但和h無關(guān)；是s維列向量，其分量為）穩(wěn)定的充要條件是：矩陣族 （4.2.26）一致有界。 命題4.2.6 矩陣族（4.2.26）有界的必要條件是的譜半徑 （4.2.27）即條件成立。 對于向前差分格式用Fourier方法判別穩(wěn)定的條件是；對于六點(diǎn)對稱

16、格式用Fourier方法判別是無條件穩(wěn)定的；對于格式絕對不穩(wěn)定。4.2.3 判別差分格式穩(wěn)定的代數(shù)準(zhǔn)則 命題4.2.7 設(shè)關(guān)于于連續(xù)，則矩陣族（4.2.26）一致有界的充要條件是矩陣族一致有界。 命題4.2.8 （一致對角化）若對任一G，有矩陣H使且H和H關(guān)于和充分小的一致有界，則條件也是穩(wěn)定的充分條件。5數(shù)值例子 用六點(diǎn)對稱格式計算拋物型方程 滿足初始條件 ，和邊界條件，的近似解，并與解析解進(jìn)行比較。 由于熱傳導(dǎo)拋物型方程的計算結(jié)果為數(shù)據(jù)網(wǎng)格，數(shù)據(jù)量較大，計算易出現(xiàn)錯誤，所以下面的差分計算采用matlab進(jìn)行數(shù)值模擬。 由于數(shù)據(jù)較大，只截取t=0.00625時的部分?jǐn)?shù)據(jù)用作對比，下面將解析解

17、與數(shù)值解的結(jié)果做成表格：h精確解 數(shù)值解 誤差（） 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.2942 0.5595 0.7701 0.9054 0.9519 0.9054 0.7701 0.5595 0.2942 0.2941 0.5595 0.7701 0.9053 0.9518 0.9053 0.77010.5595 0.2941 0.2964 0.5638 0.7760 0.9122 0.9592 0.9122 0.7760 0.5638 0.2964 6總結(jié) 本文研究了偏微分方程中拋物型方程的幾種常見的有限差分格式，并用直接法，F(xiàn)ourier方法以及

18、判別穩(wěn)定性的代數(shù)準(zhǔn)則等對這幾種格式的穩(wěn)定性進(jìn)行了判斷，其中六點(diǎn)對稱格式精度較高且是無條件穩(wěn)定的，比較適宜作為科學(xué)計算中最常用的差分格式。參考文獻(xiàn)1 陸金甫、關(guān)治.偏微分方程數(shù)值解法（第二版）M.清華大學(xué)出版社，20042 胡建偉、湯懷民.偏微分方程數(shù)值方法M.科學(xué)出版社，20053 張文生.科學(xué)計算中的偏微分方程有限差分法M.高等教育出版社，20064 余德浩、湯華中.偏微分方程數(shù)值解法M.科學(xué)出版社，20055 周蜀林.偏微分方程M.北京大學(xué)出版社，200715附錄精確解function p=jingque(h,t,L,T)x=0:h:1;n=0:t:0.05;p=zeros(L/h+1,T

19、/t+1);for j=1:length(x) for i=1:length(n) p(i,j)=exp(-(pi*pi*n(j)*sin(pi*x(i); endendp=pX,Y=meshgrid(0:h:1,0:t:0.05);mesh(X,Y,p)title(一維熱傳導(dǎo)方程的精確解的圖像);xlabel(x);ylabel(t);zlabel(一維熱傳導(dǎo)方程的精確解);endp=jingque(0.05,0.00125,1,0.05);追趕法function x=zhuigan(L,D,U,b)n=length(D);m=length(b);n1=length(L);n2=length

20、(U);for i=2:n L(i-1)=L(i-1)/D(i-1); D(i)=D(i)-L(i-1)*U(i-1);endx=zeros(n,1);x(1)=b(1);for i=2:n x(i)=b(i)-L(i-1)*x(i-1);endx(n)=x(n)/D(n);for i=n-1:-1:1 x(i)=(x(i)-U(i)*x(i+1)/D(i);end 數(shù)值解functionu,x,t=PDE(ux,ut,phi,psi1,psi2,M,N)dx=ux/M;dt=ut/N;x=(0:M)*dx;t=(0:N)*dt;r=dt/dx/dx;Diag=zeros(1,M-1);Low=zeros(1,M-2);Up=zeros(1,M-1);for i=1:M-2 Diag(i)=1+r; Low(i)=-r/2; Up(i)=-r/2;endDiag(M-1)=1+r;u=zeros(M+1,N+1);for i=1:M+1 u(i,1)=phi(x(i);endfor j=1:N+1 u(1,j)=psi1(t(j); u(M+1,j)=psi2(t(j);endB=zeros(M-1,M-1);for i=1:M-2 B(i,i)=1-r; B(i,i+