差分法求解偏微分方程MAAB

1、南京理工大學(xué)課程考核論文課程名稱：高等數(shù)值分析論文題目：有限差分法求解偏微分方程姓名：羅晨學(xué)號：成績：任課教師評語:簽名：年月日有限差分法求解偏微分方程一、主要內(nèi)容1. 有限差分法求解偏微分方程，偏微分方程如一般形式的一維拋物線型方程: 具體求解的偏微分方程如下：2. 推導(dǎo)五種差分格式、截斷誤差并分析其穩(wěn)定性；3. 編寫MATLA程序?qū)崿F(xiàn)五種差分格式對偏微分方程的求解及誤差分析；4. 結(jié)論及完成本次實驗報告的感想。二、推導(dǎo)幾種差分格式的過程:有限差分法(finite-differencemethods)是一種數(shù)值方法通過有限個微分方程近似求導(dǎo)從而尋求微分方程的近似解。有限差分法的基本思想是把連

2、續(xù)的定 解區(qū)域用有限個離散點(diǎn)構(gòu)成的網(wǎng)格來代替；把連續(xù)定解區(qū)域上的連續(xù)變量的函數(shù)用在網(wǎng)格上定義的離散變量函數(shù)來近似；把原方程和定解條件中的微商用差商來 近似，積分用積分和來近似，于是原微分方程和定解條件就近似地代之以代數(shù)方 程組，即有限差分方程組，解此方程組就可以得到原問題在離散點(diǎn)上的近似解。推導(dǎo)差分方程的過程中需要用到的泰勒展開公式如下：(xo(x X)n1)f (x0)f (x0)2f(X) f(Xo)1!(X Xo)2f(X Xo)(2-1)求解區(qū)域的網(wǎng)格劃分步長參數(shù)如下：tk 1 tkXk 1Xkh(2-2)2.1古典顯格式2.1.1古典顯格式的推導(dǎo)由泰勒展開公式將U(x,t)對時間展開

3、得u2u(Xi,t) u(x,tk) ()i,k(t tk) o(t tk) )(2-3)當(dāng)t tk 1時有u2u(Xi,tk 1)u(Xi,tk) ( )i,k(tk1 tk) o(tk 1 tk)t(2-4)u(Xi,tk)葉)i,ko( 2)得到對時間的一階偏導(dǎo)數(shù),Uu(Xi,tk 1) u(Xi,tk)()i,k=o( )(2-5)由泰勒展開公式將u(x,t)對位置展開得u(x,tk)u(Xi,tk)u()i,k(x X)X2!2u2)i,k(xXx)23o(x Xi)(2-6)當(dāng)x xi 1和 x xi 1時，代入式(2-6)得x)3) 1 “2u、i,k (xi 1、2o(Xi 1

4、Xi)()Xi)2!X、1 2u、i,k (xi 1、20 1xi)()2!X(2-7) 汀)u& i,tk) u(Xi,tk)u(N i，tk) u(X，tk)FU)i,khX(上)i,khXh2h2o(h3)(2-8)o(h3)得到對位置的二階偏導(dǎo)數(shù)22 丿i,kXu(Xi 1 ,tk) 2u(Xi,tk) u(X i,tk)h1o(h2 )(2-9)將式(2-5)、(2-9)代入一般形式的拋物線型偏微分方程得(2-10為了方便我們可以將式(2-10)寫成k 1kuiuikkku 1 2ui ui 1fik(2-11)h2kui1kkui( 2 ui 1hkk2uiui 1fik (2-1

5、2)最后得到古典顯格式的差分格式為k 1ui(1 2ra)u：kkr ui 1 ui1fik(2-13)其中r-，古典顯格式的差分格式的 截斷誤差是o(h2)h2.1.2古典顯格式穩(wěn)定性分析古典顯格式(2-13)寫成矩陣形式為k 1Uh2ra IraCkUhhk(2-14)上面的C矩陣的特征值是:2cos( j h)2raraC= 1 2racos(j h)2 j h4ra sin 2ra 12 ra cos( j h)(2-15)1,2,N 1u(Xi,tQ u(Xi,tQ U)i,k(Xi i Xu(X i,tk) u(Xi,tk) ()i,k(x-i iX使(H)1，即因為Xk 1 Xk

6、h，代入上式得121結(jié)論：當(dāng)0 ra丄時，所以古典顯格式是穩(wěn)定的22.2古典隱格式2.2.1古典隱格式的推導(dǎo)將t tk i代入式(2-3)得U2u(Xj,tkJu(Xj,tk) (t)j,k(tki tk) o(tk 1 tk) )(2-16)U2U(Xj,tk i) U(Xj,tk) ()j,k0( ) (2-17)得到對時間的一階偏導(dǎo)數(shù)uu(Xj ,tk)-)j,k=u(Xj,tk 1)0( )(2-18)將式(2-9)、(2-18)原方程得到(2-19)為了方便把(2-19)寫成k k 1Uj ujkkuj 1 2ujuj 1h2f；(2-20)k k 1kuj ujL2ukkuj 1最

7、后得到古典隱格式的差分格式為(1 2ra)u： rkkuj 1 uj 1k 1ujfjk (2-21)kfj (2-22)其中r -，古典隱格式的差分格式的 截斷誤差是o(h2)h2.2.2古典隱格式穩(wěn)定性分析將古典隱格式(2-22)寫成矩陣形式如下1 2ra I raC u： 1 u：fhk(r 2) (2-23)h誤差傳播方程1 2ra I raC V 1 v；(2-24)所以誤差方程的系數(shù)矩陣為使(H)1，顯然jH 1恒成立。結(jié)論：對于r 0,即任意網(wǎng)格比下，古典隱格式是絕對穩(wěn)定的2.3Richardson 格式2.3.1Richardso n 格式的推導(dǎo)將t tk i和t tk!，代入

8、式(2-3)得u(Xi,tkJu(x,tQF;)i,k(tk itk)o(tk itk)2)t(2-25)u(Xi,tk i)u(Xi,tk)(u)i,k(tkitk)o(tk 1tk)2)即u(,tk i) u(0tk)( +)i,k0( 2)t (2-26)u(x,tki) u(x,tk) (_U)i,ko( 2)由此得到可得2) (2-27)葉)仆 U(XtkiL U(Xitki) o(將式(2-9)、(2-27)代入原方程得到下式 (2-28)為了方便可以把式(2-28)寫成k 1 k 1UiUi2Uiki2Uikh2kUi ifik (2-29)k 1Uik 1Uikk-7 Ui 1

9、 2Ui hkUi ifik (2-30)最后得至U Richards on顯格式的差分格式為k 1Ui2ruk 1 2ukkUi i2 fik(2-31)1其中r -，古典顯格式的差分格式的 截斷誤差是o( 2 h2) h2.3.2Richards on 穩(wěn)定性分析將Richards on顯格式(2-31)寫成如下矩陣形式k 1Uh2r誤差傳播方程矩陣形式kVhkVhC 2IkUh2 fhk (2-32)2r CkVh2IkVhk 1Vhh (2-33)再將上面的方程組寫成矩陣形式k 1Vk-V2ra(CI21)1 wk (2-34)0系數(shù)矩陣的特征值是8rasin21 0(2-35)解得特

10、征值為8ra sin2 hJ64r2a2 sin4 h 41,2(2-36)Max 1 , 2 =4ra sin2 七+1+16r 2a2 si門4占衛(wèi)1 (恒成立)(2-37)結(jié)論：上式對任意的網(wǎng)比都恒成立，即 Richards on格式是絕對不穩(wěn)定的。4.Cra nk-Nichols on格式3.4.1Cra nk-Nicholso n格式的推導(dǎo)將t tk 1代入式(2-9)得u(Xj,tkJ u(Xj,tk) u)j,k (tk 1 tk) o(tk 1 tk)2)t(2-40)u2u(Xj,tk i) u(Xj,tk+1)()j,k1(tk1 tk+1)o(tk 1tk+1)即u2u(

11、Xj,tk 1)u(Xj,tk) ()j,k 2 o()2t 2(2-41)u(Xj,tk 1)u(Xj,tk+1)(U)j,k 1 2 o( 2)得到如下方程u 2 u(Xj,tk 1) u(Xj,tk)2()j,k=-o()(2-42)u2 u(Xj ,tk 1) u(Xj,tk 1)2(匚)j,k1二2o()所以t tk 1處的一階偏導(dǎo)數(shù)可以用下式表示：-1 uu1 2u(Xj,t)u(Xj,tJ 2u(Xj,tkJ u(Xj,tkJ22()jki (下)j,k 2二二o( 2)u(Xj,tk i) u(Xj,tk)/ 2、0()(2-43)將t tk , X Xj 1和x Xj 1代入

12、式(2-6)可以得到式(2-9);同理ttk 1, XXj 1 和x Xj 1代入式(2-6)可以得到2()j,kXh2u(Xj 1,tk 1) 2u(Xj,tk 1) u(Xj 1,tk 1) 小2、0(h ) (2-44)所以t tk 1 , Xj處的二階偏導(dǎo)數(shù)用式(2-6)、(2-44)表示：(2-45)所以t tk 1， Xj處的函數(shù)值可用下式表示：1f (Xj,tk 1) f(Xj,tk)(2-46)2原方程變?yōu)椋? 2 2 12(十爪1 葉爪 $ 弓)j，k (專)j,k1- f?2h2h2122-f (Xj,tk 1) f(Xj,tk)0(h )(2-48)為了方便寫成：k 1k

13、r k 1 k 1k 1kk k1 k1 kUjUj Uj 12UjUj 1Uj 12UjUj 1- fjfj(2-49)最后得到Crank-Nicholson格式的差分格式為/k 1 r k 1 k1 k r k k1- k 1- k /o 廠c r ujuj 1 uj 1 = 1 r ujuj 1 uj 1 +fjf j(2-50)其中r ，Crank-Nicholson格式的差分格式的 截斷誤差是o( 2 h2)。h3.4.1Cra nk-Nicholso n穩(wěn)定性分析Crank-Nicholson格式寫成矩陣形式如下：f?。52)(2-47)將差分格式代入上式得：u(Xj,tk 1)u

14、(Xj,tk)u(Xj1,tk 1)2u(Xj,tk1)u(Xj1,tk1) u(Xj1,tk)2u(Xj,tk)u(Xj“tQr(1 r )I TCkUh1= (1rr )I TCk1Uh+2hk(2-51)誤差傳播方程是:(1)1k 1Vh = (1 r)I T vk(2-52)可以得到:(11r)1(1 r )I C (2-53)2(1ra cos( j h)h)(1 r ) ra cos( j(2-54)2 j h1 2ra sin -2使(H)1即2rasi n* 1 2U22 j h1 2rasi n2_1 2rasin2h21 2rasin2-21 (2-55)1 2rasi n

15、2U(2-56)22 j h2 j h1 2ra sin1 2ra sin21 2rasin2-222 h (2-57)1 2rasin-22rasin2-2.2 j h ra sin2rasin2-2(2-58)1 ra sin2 2 2上式恒成立。結(jié)論：Crank-Nicholson格式對任意網(wǎng)格比也是絕對穩(wěn)定的。5.DuFortFra nkle格式（Richardson格式的改進(jìn)）將 Uik 1(uik1 Uik1)代入式(2-31)并化簡得到 DuFortFrankle :2 f(2-59)k 1kk 1 k 1 kk 1q 2r q 1 q q * 1*(1 2ra)u 1 2r u

16、-1(1 2ra)u： 12 f(2-60)可以證明DuFortFra nkle格式是絕對穩(wěn)定的。因為此格式是Richards on格式的改進(jìn)格式，因此截斷誤差還是o( 2 h2)。3.6總結(jié)(1) 古典顯格式古典顯格式的差分格式為截斷誤差：o(h2)。1穩(wěn)定性：當(dāng)0 ra 1時，古典顯格式穩(wěn)定。2(2) 古典隱格式古典隱格式的差分格式為截斷誤差：o(h2)。穩(wěn)定性：任意網(wǎng)格比古典隱格式絕對穩(wěn)定。(3) Richardso n 顯格式Richards on顯格式的差分格式為截斷誤差：o( 2 h2)。穩(wěn)定性：任意網(wǎng)格比Richards on格式絕對不穩(wěn)定。(4) Cra nk-Nicholso

17、 n 格式Crank-Nicholson格式的差分格式為截斷誤差：o( 2 h2)。穩(wěn)定性：Crank-Nicholson格式對任意網(wǎng)格比絕對穩(wěn)定。(5) DuFortFra nkle 格式截斷誤差：o( 2 h2)。穩(wěn)定性：DuFortFrankle格式對任意網(wǎng)格比絕對穩(wěn)定。三、MATLA實現(xiàn)五種差分格式對偏微分方程的求解及誤差分析3.1精確數(shù)值解上述偏微分方程的精確解是區(qū)域取值范圍：0 x 1 , 0 t 0.2。用MATLAB寸精確解進(jìn)行編程畫出三維圖像精確解程序如下：closeallclcx,t=meshgrid(0:0.01:1,0:0.001:0.2)u=exp(-pi*pi*t)

18、.*sin(pi*x)mesh(x,t,u);surf(x,t,u);xlabel(x),ylabel(t),zlabel(u);title( 精確數(shù)值解 );shadinginterp圖 3-1 精確數(shù)值解的三維圖(a)精確數(shù)值解X-Y平面(b)精確數(shù)值解X-Z平面圖 3-2 精確數(shù)值解的三個平面圖3.2五種差分格式MATLAB?序3.2.1 古典顯格式 closeall clc T=0.2 X=1.0 M=41 N=11、曰u=zeros(M,N);%構(gòu)造一個M行 N列的矩陣用于存放時間 t 和變量 x尸=(f(M穩(wěn)(定性 fprintf( 穩(wěn)定性 disp(ra);% 顯示網(wǎng) fori=

19、2:N-1 xx=(i-1)*(X/(N-1);% 網(wǎng)格比 =ra 為:n);(c) 精確數(shù)值解 Y-Z 平面und：%即唯時刻賦值，邊界條件 fork=1:M u(k,1)=0; u(k,N)=0;end;%x=0,x=1 處的邊界條件fork=1:M-1% 矩陣是從 y 軸表示行 k， x 軸表示列 i 。由 行開始fori=2:N-11);%此處為古2典顯格廣叭片沁嚇1*叭- end end將區(qū)域劃分成網(wǎng)格對每個disp(u);% 顯示差分法求得的值點(diǎn)賦值再畫畫圖=meshgrid(1:N,1:M);%將區(qū)域劃分成網(wǎng)格對每個surf(x,t,u);xlabel(x),ylabel(t),

20、zlabel(u);顯格title( 古典顯格式 );% 此程序得到的是圖 3-3圖 3-3 古典顯格式程序結(jié)果圖( r=0.5 )圖 3-4 精確數(shù)值解、古典顯格式程序結(jié)果的 Y-Z 平面圖r=0.5 )圖 3-5 古典顯格式在取不同網(wǎng)格比時的誤差傳播結(jié)果圖 圖 3-6 不同時間取值時精確解、 與古典顯格式的值對比圖 （網(wǎng)格比 r=0.5 ）紅線表示精確解、藍(lán)色線表示差分格式的解圖 3-1 、圖 3-3 對比可以看出，精確解和古典顯格式（網(wǎng)格比 r=0.5 ）的圖形是一致的。圖3-4精確數(shù)值解、古典顯格式的Y-Z平面圖結(jié)果可以看出古典顯 格式的邊界值和精確解一樣。圖 3-5 是 r 分別取

21、0.245、 0.5、 0.72、 1.125 時的 誤差傳播圖像，邊界位置網(wǎng)格數(shù)為 5 處的誤差為 0.01 得到的，可以看出 r 小于 等于 0.5 是穩(wěn)定的；但是 r 大于 0.5 出現(xiàn)不穩(wěn)定現(xiàn)象。 很好的驗證了古典顯格式 穩(wěn)定性。3.2.2 古典隱格式closeallclcT=0.2X=1.0M=41N=21航nT/(M-1穩(wěn)定性系數(shù)邈 =%a為:n網(wǎng)格比 disp(ra);% 顯示網(wǎng)格比時刻的賦值，邊界條件end;%x=0,x=1 處的邊界條件A=zeros(N-1,N-1);% 隱格式的矩陣形式中的A(1,1)=1+2*ra;A(N-1,N-1)=1+2*ra;A(1,2)=-ra

22、;A(N-1,N-2)=-ra; form=2:N-2 A(m,m-1)=-ra;A(m,m)=1+2*ra;A(m,m+1)=-ra; end;%以下是追趕法求 u 值d=zeros(N-1,1);% 隱格式右邊初始矩陣 n=length(d);U=zeros(n);L=eye(n);y=zeros(n,1);x=zeros(n,1);fori=1:N-1A 矩陣賦十、自u=zeros(M,N);%構(gòu)造一個M行 N列的矩陣用于存放時間 t和變量x fori=2:N-1 xx=(i-1)*(X/(N-1); u(1,i)=sin(pi*xx);%t=0 end; fork=1:M u(k,1)

23、=0; u(k,N)=0;角線點(diǎn)賦值再畫end%|格式右邊矩陣賦值訓(xùn) 9以下循環(huán)對矩陣A進(jìn)行LU分解U(1,1)=A(1,1)； fori=2:n L(i,i-1)=A(i,i-1)/U(i-1,i-1);U(i-1,i)=A(i-1,i);%U 的上次對角線即為 A 的上次對 U(i,i)=A(i,i)-L(i,i-1)*U(i-1,i);end fork=1:M-1%外層大循環(huán)是求時間網(wǎng)格 2,3 ,M的求%以下是追趕法的追求的解過過程程 y%(1-,-1-)-=-d-(1,1)追; 的過程 fori=2:n y(i,1)=d(i,1)-L(i,i-1)*y(i-1,1);end % 趕的

24、過程 x(n,1)=y(n,1)/U(n,n);fori=n-1:-1:1 心)=叫盤%追趕滋束加； fori=1:n u(k+1,i)=x(i,1) end d=x%每次外循環(huán)更換右邊矩矩陣end fork=1:M u(k,1)=0; end; disp(u);% 顯示差分法求得的值 x,t=meshgrid(1:N,1:M);% 圖 surf(x,t,u);xlabel(x),ylabel(t),zlabel(u);title( 古典隱格式 );% 此程序得到圖是 3-7 圖 3-7 古典隱格式程序結(jié)果圖( r=2 )即 Ly=d 的求解 y即 Ux=y 的求解 x將區(qū)域劃分成網(wǎng)格對每個圖

25、 3-8 精確數(shù)值解、古典隱格式程序結(jié)果的 Y-Z 平面圖（r=2）圖 3-9 古典隱格式在取不同網(wǎng)格比時的誤差傳播結(jié)果圖圖 3-10 不同時間取值時精確解、與古典隱格式的值對比圖（網(wǎng)格比 r=2 ）紅線表示精確解、藍(lán)色線表示差分格式的解圖 3-7 古典隱格式在 r=2 的圖形與精確解是一致的。 圖 3-8 精確數(shù)值解、 古典隱格式的 Y-Z 平面圖結(jié)果可以看出古典隱格式在 t=0.2 處的值的邊界值和精確解還是有誤差的。圖3-5是r分別取0.245、0.5、0.72、1.125時的誤差傳播圖像，邊界位置網(wǎng)格數(shù)為 5 處的誤差為 0.01 得到的，可以看出 r 取不同的值時都 是穩(wěn)定的；即古典

26、隱格式對任意的網(wǎng)格比穩(wěn)定性。從圖 3-10 可以看出隱格式隨著時間的增加， 差分格式計算的結(jié)果和精確結(jié)果越來越大； 隱格式雖然對任意網(wǎng) 格比都是穩(wěn)定的，但是計算的精確度是它的不足。3.2.3Richardson 顯格式closeallclcT=0.2X=1.0000M=41N=11ra=(T/(M-1)/(兇(N-1)八2); fpriAtd穩(wěn)定性系數(shù)S=ra為:n);u=zeros(M,N);%構(gòu)造一個 M行N列的矩陣 fori=2:N-1 xx=(i-1)*(X/(N-1);u(1,i)=sin(pi*xx);% 邊界條件 end;fork=1:Mu(k,1)=0;u(k,N)=0;% 邊

27、界條件 end;%因為 Richadson 格式需要知道前兩行的值，第二行值我們采用 %求下一面行是，隱此格處式不求再解做第注二釋行，和3.2.2 隱格式的程序一樣，只是A=zeros(N-1,N-1);A(1,1)=1+2*ra;A(N-1,N-1)=1+2*ra;A(1,2)=-ra;A(N-1,N-2)=-ra; form=2:N-2 A(m,m-1)=-ra;A(m,m)=1+2*ra;A(m,m+1)=-ra; end;n=N-1; U=zeros(n); L=eye(n); y=zeros(n,1); x=zeros(n,1);U(1,1)=A(1,1); fori=2:n L(i

28、,i-1)=A(i,i-1)/U(i-1,i-1);角線U(i-1,i)=A(i-1,i);%U 的上次對角線即為 A 的上次對U(i,i)=A(i,i)-L(i,i-1)*U(i-1,i); end y(1,1)=u(1,1);fori=2:ny(i,1)=u(1,i)-L(i,i-1)*y(i-1,1);end x(n,1)=y(n,1)/U(n,n); fori=n-1:-1:1 x(i,1)=(y(i,1)-U(i,i+1)*x(i+1,1)/U(i,i);end fuo(r2i=,12):=N0-;1 eu(n2d,i)=x(i,1)%通過隱格式已求得第二行的值u(2,i)，下面是第

29、 f行開始到第M行Richardson 格式的過程 fork=2:M-1% 時間 fori=2:N-1% 位置 u(k+1,i)=2*ra*(u(k,i-1)-2*u(k,i)+u(k,i+1)+u(k- 1,i);%Richardson 格式 end end disp(u);% 顯示求解的值 x,t=meshgrid(1:N,1:M);% 區(qū)域劃分進(jìn)行賦值畫圖 surf(x,t,u); xlabel(x),ylabel(t),zlabel(u);title(Richardson 格式 );% 此程序得到的圖形是圖 3-11圖 3-11Richardson 顯格式程序結(jié)果圖( r=0.5 )圖

30、 3-12 精確數(shù)值解、 Richardson 顯格式程序結(jié)果的 Y-Z平面圖( r=0.5 )圖 3-13Richardson 顯格式在取不同網(wǎng)格比時都 不穩(wěn)定 的結(jié)果圖圖 3-11 是 Richardson 顯格式在 r=0.5 時的程序結(jié)果圖，可以看出不穩(wěn) 定。圖3-12是精確數(shù)值解、Richardson顯格式程序結(jié)果的Y-Z平面圖。從圖3-13 可以看出 Richardson 顯格式在取不同網(wǎng)格比時都出現(xiàn)不穩(wěn)定現(xiàn)象，和理論結(jié)果 相一致。所以說 Richardson 顯格式絕對不穩(wěn)定，這種差分格式不能使用。后面 有改進(jìn)的格式 D-F 格式。3.2.4Crank-Nicholson 格式c

31、loseallclcT=0.2X=1.0000M=41N=21ra=(T/(M-1)/(X/(N-1)八2);fprintf(穩(wěn)定性系數(shù) S=ra 為:n);、 u=zeraS(M,N);%構(gòu)造一個M行 N列的矩陣用于存放時間 t 和變量 x%disp(u);fori=2:N-1xx=(i-1)*(X/(N-1);u(1,i)=sin(pi*xx);% 邊界條件end;fork=1:Mu(k,1)=0; u(k,N)=0;% 邊界條件 end;%Crank-Nicholson“格式需要兩個矩陣，下面的A、B,參照Crank-Nicholson 格式的矩陣形式A=zeros(N-1,N-1);%

32、 下面對 A 進(jìn)行填充賦值 A(1,1)=1+ra;A(N-1,N-1)=1+ra;A(1,2)=-ra/2;A(N-1,N-2)=-ra/2; form=2:N-2A(m,m-1)=-ra/2;A(m,m)=1+ra;A(m,m+1)=-ra/2;end;B=zeros(N-1,N-1);% 下面對B矩陣進(jìn)行填充賦值 B(1,1)=1-ra; B(N-1,N-1)=1-ra;B(1,2)=ra/2;B(N-1,N-2)=ra/2; form=2:N-2 B(m,m-1)=ra/2; B(m,m)=1-ra;B(m,m+1)=ra/2;分解L、Uend; %以下是追趕法求 u 值 d=zero

33、s(N-1,1);% 首先填充右邊向量，然后進(jìn)行 n=length(d);U=zeros(n);L=eye(n); y=zeros(n,1); x=zeros(n,1); U(1,1)=A(1,1);fori=2:nL(i,i-1)=A(i,i-1)/U(i-1,i-1);角線U(i-1,i)=A(i-1,i);%U 的上次對角線即為 A 的上次對U(i,i)=A(i,i)-L(i,i-1)*U(i-1,i);end fori=1:N-1d(i,1)=sin(pi*(i-1)*(1.0/(N-1);endfork=1:M-1% 按層外圍大循環(huán)，即時間步長m=zeros(n,1); m(1,1)

34、=B(1,1)*d(1,1)+B(1,2)*d(2,1););m(N-1,1)=B(N-1,N-2)*d(N-2,1)+B(N-1,N-1)*d(N-1,1 fori=2:N-2*d(i+1,1);end%上是右邊矩陣的填充更新卄亠 口 H+占厶f% 追求解Ly=b,其中b是原來的列向量矩陣乘上B系數(shù)矩陣得到的7y(1,1)=m(1,1);fori=2:nm(i,1)=B(i,i-1)*d(i-1,1)+B(i,i)*d(i,1)+B(i,i+1)y(i,1)=m(i,1)-L(i,i-1)*y(i-1,1); end% 趕 求解 Ux=yx(n,1)=y(n,1)/U(n,n);fori=n

35、-1:-1:1 x(i,1)=(y(i,1)-U(i,i+1)*x(i+1,1)/U(i,i); endfori=1:n u(k+1,i)=x(i,1) endd=zeros(N-1,1);% 右邊向量d=xendfork=1:Mu(k,1)=0;end;disp(u);%u 的值全部求出，以下畫圖 x,t=meshgrid(1:N,1:M);surf(x,t,u); xlabel(x),ylabel(t),zlabel(u);title(Crank-Nicholson格式 );% 此程序得到的圖像是圖 3-14圖 3-14Crank-Nicholson 格式程序結(jié)果圖( r=2 )圖 3-1

36、5 精確數(shù)值解、 Crank-Nicholson 格式程序結(jié)果的 Y-Z 平面圖（ r=2 ）圖 3-16Crank-Nicholson 格式在取不同網(wǎng)格比時的誤差傳播結(jié)果圖圖 3-17 不同時間取值時精確解、與 C-N 格式的解對比圖（網(wǎng)格比 r=2 ）紅線表示精確解、藍(lán)色線表示差分格式的解圖 3-14 是程序運(yùn)行得到的 Crank-Nicholson 格式在網(wǎng)格比取 r=2 的結(jié)果， 和精確解圖像一致。在 t=0.2 時從圖 3-15 的 Y-Z 面可以看出結(jié)果還是有一定的 誤差。理論上 Crank-Nicholson 格式對任意的網(wǎng)格比也是穩(wěn)定的，從圖 3-16 可 以看出在 r 取 0

37、.245 、0.5 、0.72 、1.125 誤差傳播圖像可以看出誤差不會擴(kuò)撒。 圖 3-17 是不同時間取值時精確解、與 C-N 格式 r=2 時的解對比圖，計算還是具 有誤差。3.2.5DuFortFrankle 格式closeallclcT=0.2X=1.0000M=41N=21ra=(T/(M-1)/(X/(N-1)八2);fprintf(穩(wěn)定性系數(shù) S=ra 為:n);u=zeros(M,N);%構(gòu)造一個 M行N列的矩陣 fori=2:N-1xx=(i-1)*(X/(N-1);u(r)=sin(pi*xx);%i表示 x 軸，k 表示 y 軸(即 t)feonrdk ;= 1 : M

38、 %其實初始矩陣已經(jīng)將 i=1 和 i=N 列的初值賦 u(k,1)=0;u(k,N)=0;%第二行用隱格式求得A=zeros(N-1,N-1);% A(1,1)=1+2*ra;A(N-1,N-1)=1+2*ra;零了個賦值角線面對 A 進(jìn)行填充賦值A(chǔ)(1,2)=-ra;% 第一行第二個和最后一行倒數(shù)第二個一A(N-1,N-2)=-ra; form=2:N-2 A(m,m-1)=-ra;A(m,m)=1+2*ra;A(m,m+1)=-ra;end;n=N-1;U=zeros(n);L=eye(n);y=zeros(n,1);x=zeros(n,1); U(1,1)=A(1,1);fori=2:

39、nL(i,i-1)=A(i,i-1)/U(i-1,i-1);U(i-1,i)=A(i-1,i);%U 的上次對角線即為 A 的上次對 U(i,i)=A(i,i)-L(i,i-1)*U(i-1,i);end% 追 y(1,1)=u(1,1);fori=2:n y(i,1)=u(1,i)-L(i,i-1)*y(i-1,1);end% 趕x(n,1)=y(n,1)/U(n,n);fori=n-1:-1:1 x(i,1)=(y(i,1)-U(i,i+1)*x(i+1,1)/U(i,i);endfuo(r2i=,12):=N0-;1end%第二行求出，下面用d-f差分格式fork=2:M-1fori=2

40、:N-1u(k+1,i)=(2*ra*u(k,i-1)+2*ra*u(k,i+1)+(1-2*ra)* u(k-1,i)/(1+2*ra);end end disp(u);x,t=meshgrid(1:N,1:M);surf(x,t,u);xlabel(x),ylabel(t),zlabel(u);title(DuFortFrankle 格式 );% 此程序為 Richardson 格式的改進(jìn)，得到圖 3-18圖 3-18DuFortFrankle 格式程序結(jié)果圖( r=2 )圖 3-19 精確數(shù)值解、 DuFortFrankle 格式程序結(jié)果的 Y-Z 平面圖（ r=2）圖 3-20DuFo

41、rtFrankle 格式在取不同網(wǎng)格比時的誤差傳播結(jié)果圖圖 3-21 不同時間取值時精確解、與 D-F 格式的解對比圖（網(wǎng)格比 r=2 ）紅線表示精確解、藍(lán)色線表示差分格式的解D-F 格式是 Richardson 格式（絕對不穩(wěn)定）的改進(jìn)，從圖 3-18 可以看出 當(dāng)r時D-F格式是穩(wěn)定的；圖3-20表示D-F格式網(wǎng)格比r為0.245、0.5、0.72、 1.125 時誤差傳播圖像， 不同的網(wǎng)格比誤差都不會擴(kuò)撒。 說明這種格式是穩(wěn)定的， 和理論上的結(jié)果相一致。圖3-21是不同時間取值時精確解、與 D-F格式的解對 比，隨時間的增加計算的值和差分得到的值有誤差。 此種格式雖然是絕對穩(wěn)定的， 但是

42、計算的精度還是有待提升。四、結(jié)論及感想從程序得出的結(jié)果與精確解的對比來看和理論是上的結(jié)果基本一致。 比如古 典顯格式網(wǎng)格比r小于等于0.5 （偏微分方程的系數(shù)a取1）才穩(wěn)定，從MATLAB 編程運(yùn)行的結(jié)果也可以看出 r 小于等于 0.5 是穩(wěn)定的， r 大于 0.5 不穩(wěn)定。又如 Richardson 格式理論上是絕對不穩(wěn)定的，從編程的結(jié)果在取不同的網(wǎng)格比 Richardon 格式都是不穩(wěn)定的，理論和結(jié)果一致。理論上對不穩(wěn)定格式進(jìn)行改進(jìn) 使其穩(wěn)定，比如得到的D-F格式，取不同的格式網(wǎng)格比都是穩(wěn)定的， 很好的驗證 了改進(jìn)的格式的穩(wěn)定性。本次報告首先推導(dǎo)了一維拋物線型偏微分方程的一般差分格式。 分

43、別是古典 顯格式、古典隱格式、Richardson顯格式、C-N格式進(jìn)版的Richardson格式D-F 格式。推導(dǎo)中得到各種格式的截斷誤差， 從理論上分析了各種格式的穩(wěn)定性。 對 于不穩(wěn)定的格式進(jìn)行修改得到穩(wěn)定的格式，即 Richardson 格式的修改通過MATLABS程實現(xiàn)了各種格式的程序?qū)崿F(xiàn)。用實驗的方法來驗證理論結(jié) 果，能更好的理解各種差分格式的穩(wěn)定性及操作過程。 報告中的程序都是基本程 序，誤差圖與二維 x-u 的圖像（見附錄） 都是在基本程序的基礎(chǔ)上對參數(shù)的修改 得到的圖像。通過這次報告的完成學(xué)到了很多的東西。首先，對編程有了進(jìn)一步的了解； 尤其是使用MATLA編程。在這個過程中

44、也遇到了很多的問題，通過查閱資料并 利用網(wǎng)絡(luò)資源尋求解決辦法。 其次，對于差分格式在程序上的實現(xiàn)并不是簡單的 書寫，需要步步銜接好； 比如好幾種格式都用到差分格式的矩陣形式， 尤其是隱 式格式不能直接求解， 需要應(yīng)用追趕法進(jìn)行求解； 編寫過程中通過直接求解得出 的結(jié)果都不正確， 通過追趕法才能得到正確的結(jié)果。 最后， 我們專業(yè)是電磁場與 微波技術(shù)且偏計算， 經(jīng)常遇到偏微分方程的求解； 所以這次試驗毫無疑問的對我 們理解有限差分法和求解方程提供了很大的幫助。最后，感謝張老師在課堂上的知識的講解； 同時也感謝在完成報告期間對我 提供幫助的同學(xué)！誤差傳播三維圖（以顯格式（圖closeallclc附錄

45、T=0.2X=1.0M=410.01ssuurbfp(xlo,tt,(u2),;2,1)xlabel(x),ylabel(t),zl行 撚/ 古典顯格式-誤差傳播和c245)；disp（ra）;% 顯示網(wǎng)格比%即,t=0時1 刻賦值，誤差0.01T=0.2X=1.0M=41N=11u=zeros(M,N)ra=(T/(M-1)/(X/(N-1)八2)穩(wěn)定 性 系 數(shù) S=raudi(s1p, 5(r)a=)0; . 0 1 ;% 即 t=0 時刻賦值，誤差 0.01軸 fork=1:Muu(kk,1N)=00;fork=1:M-1 fori=2:N-1endenddisp(u);taitblee(l(u); 古典顯格式 - 誤差傳播 r=0.72);clcT=0.2X=1.0M=41N=16u=zeros(M,N)