1、南京理工大學課程考核論文課程名稱：高等數值分析論文題目：有限差分法求解偏微分方程姓名：羅晨學號：115104000545成績：任課教師評語:簽名:有限差分法求解偏微分方程一、主要內容1 .有限差分法求解偏微分方程，偏微分方程如一般形式的一維拋物線型方程:-ot-2=f(x,t)(其中a為常數).:t.x具體求解的偏微分方程如下:_.2二u二ur-r2=0二t：x2u(x,0)=sin(二x)u(0,t)=u(1,t)=0t_02 .推導五種差分格式、截斷誤差并分析其穩定性；3 .編寫MATLAB程序實現五種差分格式對偏微分方程的求解及誤差分析;4 .結論及完成本次實驗報告的感想。二、推導幾種差

2、分格式的過程:有限差分法(finite-differencemethods是一種數值方法通過有限個微分方程近似求導從而尋求微分方程的近似解。有限差分法的基本思想是把連續的定解區域用有限個離散點構成的網格來代替；把連續定解區域上的連續變量的函數用在網格上定義的離散變量函數來近似；把原方程和定解條件中的微商用差商來近似，積分用積分和來近似，于是原微分方程和定解條件就近似地代之以代數方程組，即有限差分方程組，解此方程組就可以得到原問題在離散點上的近似解。(n)/一?(x-)no(x-x)n1)n!(2-1)(2-2)推導差分方程的過程中需要用到的泰勒展開公式如下：f(x)=f(x)flx-x。)f2

3、x0)(x-x。)2求解區域的網格劃分步長參數如下：tk1-tk=xk1一人二h2.1.1古典顯格式的推導由泰勒展開公式將u（x,t）對時間展開得u(x,tXu(ixkt4)U(,i)-t(kt+)o(kt)(2-3)；:t當t=tk中時有,、,、；：u,、，、，2、U(x,kiu(iX,kt)一(,ik&tk)o蟲k)(2-4)二u(x,k),平，i)k,o2()Ft得到對時間的一階偏導數：Uu(xi,tk1)-u(xi,tk)7、gG()i,k=o(.)(2-5).t由泰勒展開公式將u(x,t)對位置展開得/,、/，、/二u、，、1-u2.、3、u(x,tk)=u(xi,tk)()i,k(

4、x-xj(-)i,k(x-xi)o(x-x)(2-6)x2!:x當x=為由和x=為工時，代入式（2-6）得u1八：u23、(2-7)u(x+,tk)=u(xi,tk)+()i,k(xi+-xi)+()i,k(xi-xi)十。(為書xi)二x2!二xFu1：2u23u(x.tk)=口(為上)+(丁)蔣(為一為)+(-)i,k(x_1-xi)+o(xx)二x2!二x因為xk書-xk=h,代入上式得u、u(xtk)=u(x,tj()i,k二xu、u(x，tk)-u(xi,tk)-()i,kx2hLdh(22!:x2-2)i,kh2o(h3)2一一h不得)rho(h)(2-8)得到對位置的二階偏導數-

5、2(2-9)/二u、u(xi1,tk)-2u(xi,tk)u(x，tk)/、Di,k=o(h)exh將式(2-5)、(2-9)代入一般形式的拋物線型偏微分方程得u(xi,tk1)-u(x,tk)一u(x1,tk)-2u(xi,tk)u(xi4,tk)2.-二f(xi,tk)o(-h),h(2-10)為了方便我們可以將式（2-10）寫成k1kkkkuiuiUi1-2U,Uj1k一2=fihhk1k,Ji.kkkkUi-Ui2Ui1-2uiUi1=fih最后得到古典顯格式的差分格式為k1kkkkU(1-2ra)5r：工Ui1ut1其中r=J，古典顯格式的差分格式的截斷誤差是o+h2)h2.1.2古

6、典顯格式穩定性分析古典顯格式(2-13)寫成矩陣形式為Uk1=1_2raIraCUkfhk(2-11)(2-12)(2-13)(2-14)其中r，Uk=(Uk,Uk,.,u：n,uN_l)。h一01III01010C=0：101.0a10一(N)MN)上面的C矩陣的特征值是:%=2cos(jnh)j=1,2,N-1H=1-2raIraC,：二1一7ara,C=1r2ra2cjoS()二1-2ra1cojSh)(2-15)=1-4rasi2nj-hj=1,2,N二.，12使P(H)1,即-1-1-4rasin2j-h-120三ra-121結論：當0WraW，時，所以古典顯格式是穩定的。22.2.

7、1古典隱格式的推導將1=17代入式（2-3）得u(Xj,tk)=u(Xj,tk)(4j,k(tk-tk)o(tk-tk)2)ctU(Xj,tk)=U(Xj,tk)-(U)j，ko(2)二t得到對時間的一階偏導數二u、U(Xj,tk)-U(Xj,tkj)()j,k=o()ft(2-16)(2-17)(2-18)將式（2-9）、（2-18）原方程得到Uga。*)一u(Xji,tk)-2u(Xj,tk)u(Xj,tk)1h一2f(Xj,tk)o(h)為了方便把（2-19）寫成kkJUj一山OlTuk書-2心+u；Ik2=fih2jkkJ-.kkkkUj-Uj-斤ujL2UjUp=fj最后得到古典隱格

8、式的差分格式為(12ra)uk-ruk.1u；=uk.fjk(2-19)(2-20)(2-21)(2-22)其中r=：，古典隱格式的差分格式的截斷誤差是。（T+h2）h2.2.2古典隱格式穩定性分析將古典隱格式(2-22)寫成矩陣形式如下_12raI-raCu；1=uk.f；(r=m)(2-23)誤差傳播方程12raI-raCvk1=vk(2-24)A=12raI-raC,B=I所以誤差方程的系數矩陣為.1_一H二A-|J1-2raI-raC使P(H)1,顯然112ra-2racosj二hj=1,2,N-1112ra-2racos(j二h)11+2ra(1cosjh)142j二h14rasin

9、2九H,k=；o()二t,2u(xi,tki)-U(X|&i)2將式(2-9)、(2-27)代入原方程得到下式u(X|1,tk)2u(xi,tk)u(xy,tk)一h2為了方便可以把式(2-28)寫成uk1-Uikj-a2:Ui=-2Uik+UiZ-h2k1kd-kkkkUi-Ui2Ui1-2ui,Ui1=fih最后得到Richardson顯格式的差分格式為u：*=2ra(4書一2u：+u：,廣u：，+2ttk其中r=,，古典顯格式的差分格式的截斷誤差是o(d+h2)h2.3.2 Richardson穩定性分析將Richardson顯格式(2-31)寫成如下矩陣形式u：1=2廠C-2Iu：u：

10、,2.fhk誤差傳播方程矩陣形式v：1=21C-2Iv：v一kkvh=vh再將上面的方程組寫成矩陣形式k+3)；2ra(C-2I)I：wh=-=w1VliIOJ系數矩陣的特征值是4racos(j二h)-4ra1ijil1o22j二h18rasin1=02解得特征值為-8rasin2j-h-士64r2a2sin4j-h+4,221,2一2Ji112jh|224jh八一4、Max%,%=4rasin+J1+16rasin1(包成立)(2-29)(2-30)(2-31)(2-32)(2-33)(2-34)(2-35)(2-36)(2-37)結論：上式對任意的網比都包成立，即Richardson格式是

11、絕對不穩定的4.Crank-Nicholson格式3.4.1Crank-Nicholson格式的推導將1=tf代入式(2-9)得(2-40)產)+爭k(t”k)+o(tk)2)：u2u(Xj,tk另)=u(Xj,tk+1)+(1)j,k書(tk超一tk+1)+o(tkq一tk+1)二t得到如下方程二uu(Xj,tk.i)=u(Xj,tk)(一)r二uj，k2o(2)(2-41)u(Xj,tk.1)=u(Xj,tk+i)-()二t;u(_T)j,k=二t2(u(Xj,tk用)u(Xj,tk)To(2)(2-42)：:u(-7)j,k1二t2(u(Xj,tk號)u(Xj,tk書)o(2)所以t=t

12、f處的一階偏導數可以用下式表示:1,當()j,k1()j,k2卜ttt2口（為上號）u（Xj,tk）2u(Xj,tki)-u(Xj,tk.i)o(2)u(Xj,tk1)-u(Xj,tk),2、o()(2-43)將t=tk,X=Xj書和x=Xj代入式（2-6）可以得到式（2-9）;同理t=tk+,X=Xj由和x=Xj代入式（2-6）可以得到l2(*=(-2)j,k1：Xu(xj1,tk1)2u(xj,tk1)u(xj，tk1).h2o(h2)(2-44)所以1=旁，Xj處的二階偏導數用式（2-6）、（2-44）表示:2o(h2)(2-45)j,k+(-2)j,k卡cX1lu(xj1,tk1)-2

13、U(Xj,tk1)U(Xji,tk1)U(Xj1,tk)2U(Xj,tk)U(Xj,tk)|=一222ILh2h2所以t=tk3,Xj處的函數值可用下式表示：(2-46)1-r2|f(X,t1)f(X,kt)原方程變為:1/：:uu2(力k1(力k-2,二U-Jx：2u)j,k(-r)j,k.1；x1k1k2二(fj書+fj)+o(h)(2-47)將差分格式代入上式得：tk)-2u(Xj,tk)u(Xj,tk)u(Xj,tk1)-u(Xj,tk)：-u(Xj1,tk1)-2u(Xj,tki)u(Xj,tk1)u(Xj1,h22-2ILh=；f(Xj,tkQ+f(Xj,tk)+o(T2+h2)(

14、2-48)為了方便寫成：k1kr-,-；Jk1k+1k-uj-uj21Uj1-2Uj=f：(2-49)最后得到Crank-Nicholson格式的差分格式為k1k1.一kr-：kk1Uj1uj4=1-r-Uj萬Uj1Uj+-fjk1fjk(2-50)其中r=,Crank-Nicholson格式的差分格式的截斷誤差是o(T2+h2)。h23.4.1Crank-Nicholson穩定性分析Crank-Nicholson格式寫成矩陣形式如下:(1r：)Ir:k1Cuk1=(1-r：)ICuk+.2f；1-fhk(2-51)誤差傳播方程是:可以得到:(2-52)(2-53)H=(1+rcc)H(1-r

15、二)racosj(hj=(2-54)(1r:)-racosj(h1-2asi2世=212asfnh2使P(H)1即1-2asi2-112rasin2j2h-1(2-55)22j二h2rh2一1-2asin1r2is4n1ra2-22jih2(2-56)2j二h2j二h1-2rasin12rasin-222j二h2j二h-1-2rasin1-2rasin22(2-57)2j二h2h-2rasin_2rasin222j二h2j二h-rasin-1-rasin22(2-58)結論：Crank-Nicholson格式對任意網格比也是絕對穩定的。5.DuFortFrankle格式(Richardson格

16、式的改進)將Uik=1(Uik+u尸)代入式(2-31)并化簡得到DuFortFrankle：k1kk1k-1kk-1ku=2rui1-u-uiui4ui2fi(2-59)(1+2ra)qk+=29(u：+口：=)+(1-2ra)ur+2”1k(2-60)可以證明DuFortFrankle格式是絕對穩定的。因為此格式是Richardson格式的改進格式，因此截斷誤差還是o(i2+h2)。3.6總結(1)古典顯格式古典顯格式的差分格式為u尸=(12ra)u：+r(u)+u：fj截斷誤差：o(i+h2)o1穩定性：當0raM1時，古典顯格式穩定。2(2)古典隱格式古典隱格式的差分格式為(12ra)

17、u：-:工1u：1u：=u：，.f；截斷誤差：o(i+h2)o穩定性：任意網格比古典隱格式絕對穩定。(3) Richardson顯格式Richardson顯格式的差分格式為u/=2ru(u：書2u：+u：)+qk，+2丁fik截斷誤差：o(T2+h2)。穩定性：任意網格比Richardson格式絕對不穩定。(4) Crank-Nicholson格式Crank-Nicholson格式的差分格式為1 +。)u：*-卷心；+ukJ1)=(1-rJu：+學年書+u：,)+tfjk*+f；)j截斷誤差：o(2h2)。穩定性：Crank-Nicholson格式對任意網格比絕對穩定。DuFortFrankl

18、e格式(1+2ra)uik+=2口(u+uik1)+(1-2ra)uk+2wfik截斷誤差：o(h2)。穩定性：DuFortFrankle格式對任意網格比絕對穩定。三、MATLAB實現五種差分格式對偏微分方程的求解及誤差分析3.1精確數值解一.2二u二u.:t一2-Xt,0u(x,0)=sin(二x)u(0,t)=u(1,t)=0上述偏微分方程的精確解是u(x,t)=eTtsin(二x)(0x1,t_0)區域取值范圍：0ExEl,0t0.20用MATLAB對精確解進行編程畫出三維圖像精確解程序如下:closeallclcx,t=meshgrid(0:0.01:1,0:0.001:0.2)u=e

19、xp(-pi*pi*t).*sin(pi*x)mesh(x,t,u);surf(x,t,u);xlabel(x),ylabel(t),zlabel(u);title(精確數值解);shadinginterp(a)精確數值解X-Y平面(b)精確數值解X-Z平面以前藥倡竄-0.5(c)精確數值解Y-Z平面圖3-2精確數值解的三個平面圖3.2五種差分格式MATLAB程序3.2.1 古典顯格式closeallclcT=0.2X=1.0M=41N=11u=zeros(M,N);%構造一個M行N列的矩陣用于存放時間t和變量xra=(T/(M-1)/(X/(N-1)A2);%網格比fprintf(穩定性系數

20、S=ra為:n);disp(ra);%顯示網格比fori=2:N-1xx=(i-1)*(X/(N-1);u(1,i)=sin(pi*xx);end;%即t=0時刻賦值，邊界條件fork=1:Mu(k,1)=0;u(k,N)=0;end;%x=0,x=1處的邊界條件fork=1:M-1%矩陣是從y軸表示行k,x軸表示列i。由行開始fori=2:N-1u(k+1,i)=(1-2*ra)*u(k,i)+ra*u(k,i+1)+ra*u(k,i-1);%此處為古典顯格式endenddisp(u);%顯示差分法求得的值x,t=meshgrid(1:N,1:M);%將區域劃分成網格對每個點賦值再畫圖sur

21、f(x,t,u);xlabel(x),ylabel(t),zlabel(u);title(古典顯格式);%此程序得到的是圖3-3圖3-3古典顯格式程序結果圖（r=0.5）圖3-5古典顯格式在取不同網格比時的誤差傳播結果圖時，情確顰、古曲顯格式的解曲線對比t=0125,橫礴解.古典顯格式的解曲線對比0.3i03D.2D40.G口田1圖3-6不同時間取值時精確解、與古典顯格式的值對比圖（網格比r=0.5）紅線表示精確解、藍色線表示差分格式的解圖3-1、圖3-3對比可以看出，精確解和古典顯格式（網格比r=0.5）的圖形是一致的。圖3-4精確數值解、古典顯格式的Y-Z平面圖結果可以看出古典顯格式的邊界

22、值和精確解一樣。圖3-5是r分別取0.245、0.5、0.72、1.125時的誤差傳播圖像，邊界位置網格數為5處的誤差為0.01得到的，可以看出r小于等于0.5是穩定的；但是r大于0.5出現不穩定現象。很好的驗證了古典顯格式穩定性。3.2.2古典隱格式closeallclcT=0.2X=1.0M=41N=21ra=(T/(M-1)/(X/(N-1)A2);fprintf(穩定性系數S=radisp(ra);u=zeros(M,N);fori=2:N-1%網格比為：n);%顯示網格比t和變量%構造一個M行N列的矩陣用于存放時間xx=(i-1)*(X/(N-1);%t=0時刻的賦值，邊界條件%x=

23、0,x=1處的邊界條件%隱格式的矩陣形式中的A矩陣賦值u(1,i)=sin(pi*xx);end;fork=1:Mu(k,1)=0;u(k,N)=0;end;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-2A(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,

24、1);fori=1:N-1d(i,1)=sin(pi*(i-1)*(1.0/(N-1);end%隱格式右邊矩陣賦值%以下循環對矩陣A進彳TLU分解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);endfork=1:M-1%外層大循環是求時間網格2,3,M的求解u%以下是追趕法的求解過程%追的過程即Ly=d的求解yy(1,1)=d(1,1);fori=2:ny(i,1)=d(i,1)-L(i,i-1)*y(i-

25、1,1);end%趕的過程即Ux=y的求解xx(n,1)=y(n,1)/U(n,n);fori=n-1:-1:1x(i,1)=(y(i,1)-U(i,i+1)*x(i+1,1)/U(i,i);%顯示差分法求得的值%將區域劃分成網格對每個點賦值再畫圖%此程序得到圖是3-7圖3-7古典隱格式程序結果圖（r=2）endfori=1:nu(k+1,i)=x(i,1)endd=zeros(N-1,1);d=xendfork=1:Mu(k,1)=0;end;disp(u);x,t=meshgrid(1:N,1:M);surf(x,t,u);xlabel(x),ylabel(t),zlabel(u);tit

26、le(古典隱格式);%追趕法結束%更新右邊矩陣%每次外循環更換右邊矩陣古典隱榨式0.01圖3-9古典隱格式在取不同網格比時的誤差傳播結果圖圖3-10不同時間取值時精確解、與古典隱格式的值對比圖（網格比r=2）紅線表示精確解、藍色線表示差分格式的解圖3-7古典隱格式在r=2的圖形與精確解是一致的。圖3-8精確數值解、古典隱格式的Y-Z平面圖結果可以看出古典隱格式在t=0.2處的值的邊界值和精確解還是有誤差的。圖3-5是r分別取0.245、0.5、0.72、1.125時的誤差傳播圖像,邊界位置網格數為5處的誤差為0.01得到的，可以看出r取不同的值時都是穩定的；即古典隱格式對任意的網格比穩定性。從

27、圖3-10可以看出隱格式隨著時間古四圈格式-誤差傳櫓（=0.245古印第格式誤爰傳播用石0005古擔隱格式誤差博棉胃072古隈旗格式-誤差佶播何1250005-nC.0Q5仁口。好時，精璐醉、古韓在楮式的M曲蛙對比t=C05Hl,M說時、古和旌格式的解曲踹對比D.F|t=012SE1k楂裾解-古卷降格式的解日強對比0.M.-,1=。廿5時,福宰翩、古曲踵梏式的解由蛙對比0.2S.-.的增加，差分格式計算的結果和精確結果越來越大；隱格式雖然對任意網格比都是穩定的，但是計算的精確度是它的不足。3.2.3Richardson顯格式closeallclcT=0.2X=1.0000M=41N=11ra=

28、(T/(M-1)/(X/(N-1)A2);fprintf(穩定性系數S=ra為:n);disp(ra);u=zeros(M,N);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;%構造一個M行N列的矩陣%邊界條件%邊界條件%因為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)

29、=-ra;A(N-1,N-2)=-ra;form=2:N-2A(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:nL(i,i-1)=A(i,i-1)/U(i-1,i-1);%U的上次對角線即為A的上次對角線U(i,i)=A(i,i)-L(i,i-1)*U(i-1,i);endy(1,1)=u(1,1);fori=2:ny(i,1)=u(1,i)-L(i,i-1)*y(i-1,1);endx(n,1)=y(n,1)/

30、U(n,n);fori=n-1:-1:1x(i,1)=(y(i,1)-U(i,i+1)*x(i+1,1)/U(i,i);endu(2,1)=0;fori=2:N-1u(2,i)=x(i,1)end%通過隱格式已求得第二行的值u(2,i),下面是Richardson格式的過程fork=2:M-1%時間軸第3行開始到第M行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格式endenddisp(u);%顯示求解的值x,t=meshgrid(1:N,1:M);%區域劃分進行賦值畫圖surf(x

31、,t,u);xlabel(x),ylabel(t),zlabel(u);title(Richardson格式)；此程序得到的圖形是圖3-11Rich曰rdsori格式21x1000圖3-11Richardson顯格式程序結果圖（r=0.5）圖3-13Richardson顯格式在取不同網格比時都不穩定的結果圖圖3-11是Richardson顯格式在r=0.5時的程序結果圖，可以看出不穩定。圖3-12是精確數值解、Richardson顯格式程序結果的Y-Z平面圖。從圖3-13可以看出Richardson顯格式在取不同網格比時都出現不穩定現象，和理論結果相一致。所以說Richardson顯格式絕對不

32、穩定，這種差分格式不能使用。后面有改進的格式D-F格式。3.2.4 Crank-Nicholson格式closeallclcT=0.2X=1.0000M=41N=21ra=(T/(M-1)/(X/(N-1)A2);fprintf(穩定性系數S=ra為:n);disp(ra);u=zeros(M,N);%構造一個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格式需要兩個矩陣，下面

33、的A、B,參照Crank-Nicholson格式的矩陣形式A=zeros(N-1,N-1);%下面對A進行填充賦值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矩陣進行填充賦值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-2B(m,m-1)=ra/2;B(m,m)=1-ra;B(m,m+1

34、)=ra/2;end;%以下是一一追趕法求u值d=zeros(N-1,1);%首先填充右邊向量，然后進行L、U分解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);endfori=1:N-1d(i,1)=sin(pi*(i-1)*(1.0/(N-1);endfork=1:M-1%按層外圍大循環，即

35、時間步長m=zeros(n,1);m(1,1)=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-2m(i,1)=B(i,i-1)*d(i-1,1)+B(i,i)*d(i,1)+B(i,i+1)*d(i+1,1);end%以上是右邊矩陣的填充更新%追求解Ly=b,其中b是原來的列向量矩陣乘上B系數矩陣得到的y(1,1)=m(1,1);fori=2:ny(i,1)=m(i,1)-L(i,i-1)*y(i-1,1);end%趕求解Ux=yx(n,1)=y(n,1)/U(n,n)

36、;fori=n-1:-1:1x(i,1)=(y(i,1)-U(i,i+1)*x(i+1,1)/U(i,i);endfori=1:nu(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格式程序名果圖(r=2)th電第時.，眄薛、c制格式

37、的解曲城對比時撲斛.CH格式的第曲找對比圖3-17不同時間取值時精確解、與C-N格式的解對比圖（網格比r=2）紅線表示精確解、藍色線表示差分格式的解1t聘時,精福解、C用格式的斜曲城對比圖3-14是程序運行得到的Crank-Nicholson格式在網格比取r=2的結果，和精確解圖像一致。在t=0.2時從圖3-15的Y-Z面可以看出結果還是有一定的誤差。理論上Crank-Nicholson格式對任意的網格比也是穩定的，從圖3-16可以看出在r取0.245、0.5、0.72、1.125誤差傳播圖像可以看出誤差不會擴撒。圖3-17是不同時間取值時精確解、與C-N格式r=2時的解對比圖，計算還是具有誤

38、差。3.2.5 DuFortFrankle格式closeallclcT=0.2X=1.0000M=41N=21ra=(T/(M-1)/(X/(N-1)A2);fprintf(穩定性系數S=ra為:n);disp(ra);u=zeros(M,N);%構造一個M行N列的矩陣fori=2:N-1xx=(i-1)*(X/(N-1);u(1,i)=sin(pi*xx);%i表示x軸，k表示y軸(即t)end;fork=1:M%其實初始矩陣已經將i=1和i=N列的初值賦零了u(k,1)=0;u(k,N)=0;end;%第二行用隱格式求得A=zeros(N-1,N-1);%下面對A進行填充賦值A(1,1)=

39、1+2*ra;A(N-1,N-1)=1+2*ra;A(1,2)=-ra;%第一行第二個和最后一行倒數第二個一個賦值A(N-1,N-2)=-ra;form=2:N-2A(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: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%

40、追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:1x(i,1)=(y(i,1)-U(i,i+1)*x(i+1,1)/U(i,i);endu(2,1)=0;fori=2:N-1u(2,i)=x(i,1)end%第二行求出，下面用D-F差分格式fork=2:M-1fori=2: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);endenddisp(u);x,t=meshgrid

41、(1:N,1:M);surf(x,t,u);xlabel(x),ylabel(t),zlabel(u);title(DuFortFrankle格式)；此程序為Richardson格式的改進，得到圖3-18圖3-18DuFortFrankle格式程序Z果圖(r=2)圖3-19精確數彳1解、DuFortFrankle格式程序結果的Y-Z平面圖（r=2）DuForiF幽kle格式誤差悔強（=0245DuForiFankle格式誤差拷Jfi(MJ5DuFmFenHe格式誤差傳JfifO.72DuForiFMk岫格式栽整傳播圖3-20DuFortFrankle格式在取不同網格比時的誤差傳播結果圖圖3-2

42、1不同時間取值時精確解、與D-F格式的解對比圖（網格比r=2）紅線表示精確解、藍色線表示差分格式的解D-F格式是Richardson格式（絕對不穩定）的改進，從圖3-18可以看出當r時D-F格式是穩定的；圖3-20表示D-F格式網格比r為0.245、0.5、0.72、1.125時誤差傳播圖像，不同的網格比誤差都不會擴撒。說明這種格式是穩定的，和理論上的結果相一致。圖3-21是不同時間取值時精確解、與D-F格式的解對比，隨時間的增加計算的值和差分得到的值有誤差。此種格式雖然是絕對穩定的，但是計算的精度還是有待提升。四、結論及感想從程序得出的結果與精確解的對比來看和理論是上的結果基本一致。比如古典

43、顯格式網格比r小于等于0.5（偏微分方程的系數a取1）才穩定，從MATLAB編程運行的結果也可以看出r小于等于0.5是穩定的，r大于0.5不穩定。又如Richardson格式理論上是絕對不穩定的，從編程的結果在取不同的網格比Richardon格式都是不穩定的，理論和結果一致。理論上對不穩定格式進行改進使其穩定，比如得到的D-F格式，取不同的格式網格比都是穩定的，很好的驗證了改進的格式的穩定性。本次報告首先推導了一維拋物線型偏微分方程的一般差分格式。分別是古典顯格式、古典隱格式、Richardson顯格式、C-N格式進版的Richardson格式D-F格式。推導中得到各種格式的截斷誤差，從理論上

44、分析了各種格式的穩定性。對于不穩定的格式進行修改得到穩定的格式，即Richardson格式的修改。通過MATLAB編程實現了各種格式的程序實現。用實驗的方法來驗證理論結果，能更好的理解各種差分格式的穩定性及操作過程。報告中的程序都是基本程序，誤差圖與二維x-u的圖像（見附錄）都是在基本程序的基礎上對參數的修改得到的圖像。通過這次報告的完成學到了很多的東西。首先，對編程有了進一步的了解；尤其是使用MATLAB編程。在這個過程中也遇到了很多的問題，通過查閱資料并利用網絡資源尋求解決辦法。其次，對于差分格式在程序上的實現并不是簡單的書寫，需要步步銜接好；比如好幾種格式都用到差分格式的矩陣形式，尤其是

45、隱式格式不能直接求解，需要應用追趕法進行求解；編寫過程中通過直接求解得出的結果都不正確，通過追趕法才能得到正確的結果。最后，我們專業是電磁場與微波技術且偏計算，經常遇到偏微分方程的求解；所以這次試驗毫無疑問的對我們理解有限差分法和求解方程提供了很大的幫助。最后，感謝張老師在課堂上的知識的講解；同時也感謝在完成報告期間對我提供幫助的同學！誤差傳播三維圖(以顯格式(圖3-5)closeallclcT=0.2X=1.0M=41N=8u=zeros(M,N);%構造一個M行N歹U的矩陣用于存放時間t和變量xra=(T/(M-1)/(X/(N-1)A2);%網格比fprintf(穩定性系數S=ra為:n

46、);disp(ra);%顯示網格比u(1,5)=0.01;%即t=0時刻賦值，誤差0.01fork=1:Mu(k,1)=0;u(k,N)=0;end;%x=0,x=1處的邊界條件fork=1:M-1%矩陣是從y軸表本行k,x軸表本列io由行開始fori=2:N-1u(k+1,i)=(1-2*ra)*u(k,i)+ra*u(k,i+1)+ra*u(k,i-1);%此處為古典顯格式endenddisp(u);%顯示差分法求得的值x,t=meshgrid(1:N,1:M);%將區域劃分應網格對每個點賦值再畫圖subplot(2,2,1)surf(x,t,u);xlabel(x),ylabel(t),

47、zlabel(u);title(古典顯格式-誤差傳播r=0.245,);clcT=0.2X=1.0M=41N=11u=zeros(M,N);ra=(T/(M-1)/(X/(N-1)A2);fprintf(穩定，住系數S=ra為:n);disp(ra);u(1,5)=0.01;%即t=0時刻賦值，誤差0.01fork=1:Mu(k,1)=0;u(k,N)=0;end;fork=1:M-1fori=2:N-1u(k+1,i)=(1-2*ra)*u(k,i)+ra*u(k,i+1)+ra*u(k,i-1);endenddisp(u);x,t=meshgrid(1:N,1:M);subplot(2,2,2)surf(x,t,u);附錄為例):xlabel(x),ylabel(t),zlabel(u);title(古典顯格式-誤